[R] Risk differences with survey package

[Thomas Lumley]

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