R help - Mar 2011 - Risk differences with survey package

I'm trying to use the survey package to calculate a risk difference with
confidence interval for binge drinking between sexes. Variables are
X_RFBING2 (Yes, No) and SEX. Both are factors. I can get the group
prevalences easily enough with

result <- svyby(~X_RFBING2, ~SEX, la04.svy, svymean, na.rm = TRUE)

and then extract components from the svyby object with SE() and coef() to
do the computations. This gives the correct results, but I'd like to set
this up as a contrast and am having difficulty. What is the best way to do
it?

Thanks!

On Tue, Mar 8, 2011 at 6:50 AM,  <Michael.Laviolette at dhhs.state.nh.us>
wrote:
>
> I'm trying to use the survey package to calculate a risk difference
with
> confidence interval for binge drinking between sexes. Variables are
> X_RFBING2 (Yes, No) and SEX. Both are factors. I can get the group
> prevalences easily enough with
>
> result <- svyby(~X_RFBING2, ~SEX, la04.svy, svymean, na.rm = TRUE)
>
> and then extract components from the svyby object with SE() and coef() to
> do the computations. This gives the correct results, but I'd like to
set
> this up as a contrast and am having difficulty. What is the best way to do
> it?
The easiest way is to use svyglm() -- a risk difference is what you
get in a linear model

For example, using the built-in dclus2 data set
> riskmodel<-svyglm(as.numeric(comp.imp)~sch.wide, design=dclus2)
> riskmodel
2 - level Cluster Sampling design
With (40, 126) clusters.
svydesign(id = ~dnum + snum, fpc = ~fpc1 + fpc2, data = apiclus2)

Call:  svyglm(as.numeric(comp.imp) ~ sch.wide, design = dclus2)

Coefficients:
(Intercept)  sch.wideYes
     1.1068       0.7743

Degrees of Freedom: 125 Total (i.e. Null);  38 Residual
Null Deviance:	    27.02
Residual Deviance: 12.9 	AIC: 124.8
> confint(riskmodel)
                2.5 %    97.5 %
(Intercept) 0.9484637 1.2651862
sch.wideYes 0.6232830 0.9253461


You can't extract the results from the output of svyby(), in general,
because svyby() doesn't know how to compute the covariances between
estimates in the groups -- after all, the function input to svyby()
can be any function including one you wrote.  These estimates won't
typically be independent in a complex design, in contrast to the
situation in a simple random sample.

With a replicate-weights design you can use svyby() with the
covmat=TRUE option to return the covariance matrix, and then compute
contrasts.  This only works if the function returns the replicate
estimates (which all the built-in functions do), allowing the
covariance to be computed.
> groups<-svyby(~comp.imp, ~sch.wide, design=rclus1, svymean, covmat=TRUE)
> svycontrast(groups, quote(`Yes:comp.impYes`-`No:comp.impYes`))
           nlcon     SE
contrast 0.83125 0.0209

   -thomas

-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland