Re: [R] svyby and svyratio in the Thomas Lumley's survey package

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Fri 24 Feb 2006 - 07:43:07 EST

On Thu, 23 Feb 2006, Smith, Phil wrote:

> Dear R-ers:
>
> I'm using Thomas Lumley's "survey" package.
>
> I'd like to compute survey ratio estimates (numerator=~utd ,
> denominator=~one) for each of serval domains using by=~factor(domain).
>
> I can't quite work out the syntax for the call to the "svyby" function. I try:
>
> svyby( numerator=~ utd, denominator=~one , by=~ factor(domain) ,
> design=nis , svyratio )
>
> and I get an error that says that "domain" is not found. (nis is the
> design object)
>

Hmm. I don't get the same error message (perhaps a question of versions?), but there is a problem. svyby is not really set up to work with svyratio. I will try to fix this, but there are simple workarounds.

If you do
  data(api)
  dstrat<-svydesign(id=~1,strata=~stype, weights=~pw,

        data=apistrat, fpc=~fpc)

  result <- svyby(~api.stu,by=~factor(stype), denominator=~enroll,

        design=dstrat,svyratio)

(ie, don't name the numerator argument) then the computations are done as you want. The result does not print correctly, because svyby() returns some extra attributes that print() doesn't know how to handle.

You can extract the results you need with coef() or just remove the extraneous attributes
> result

Error in unlist(x, recursive, use.names) :

         argument not a list
> coef(result)

   statistic.ratio statistic.var

E       0.8518163  4.945408e-05
H       0.8105702  0.0004193180
M       0.8356958  0.0003307829

> result$statistic.call<-NULL
> result
   factor(stype) statistic.ratio statistic.var         SE
E             E       0.8518163  4.945408e-05 0.00703236
H             H       0.8105702  0.0004193180 0.02047726
M             M       0.8356958  0.0003307829 0.01818744



Finally, the name of the denominator variable looks suspicious. If `one' is a variable whose values are all 1 and you just want domain means you can do the equivalent of

   svyby(~api00,~stype,design=dstrat,svymean) to get the mean of api00 within each value of stype. You don't need svyratio(). In fact, part of the test suite for the survey package is to check that svyratio, svymean, and svyglm give the same answers for domain means.

         -thomas



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Feb 24 08:07:22 2006

This archive was generated by hypermail 2.1.8 : Fri 24 Feb 2006 - 09:09:01 EST