[R] survey package svystat objects from predict()

[Thomas Lumley]

On Tue, Feb 14, 2012 at 11:45 AM, Kieran Healy <kjhealy at gmail.com> wrote:
> Hello,
>
> I'm running R 2.14.1 on OS X (x86_64-apple-darwin9.8.0/x86_64 (64-bit)), with version 3.28 of Thomas Lumley's survey package. I was using predict() from svyglm(). E.g.:
>
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> out <- svyglm(sch.wide~ell+mobility, design=dstrat,
>        family=quasibinomial())
> pred.df <- expand.grid(ell=c(20,50,80), mobility=20)
> out.pred <- predict(out, pred.df)
>
> From the console out.pred looks like this:
>
>> class(out.pred)
> [1] "svystat"
>
>> print(out.pred) # or just 'out.pred'
>    link     SE
> 1 1.8504 0.2414
> 2 1.6814 0.3033
> 3 1.5124 0.5197
>
> From here I wanted to conveniently access the predicted values and SEs. I thought that I might be able to do this using either ftable() or as.data.frame(), as methods for these exist for the objects of class "svystat". But while the predicted values come through fine, the SE only gets calculated for the first element and is then repeated:
>
>> ftable(out.pred)
>          A         B
> 1 1.8504120 0.2413889
> 2 1.6814293 0.2413889
> 3 1.5124466 0.2413889
>
>> as.data.frame(out.pred)
>      link        SE
> 1 1.850412 0.2413889
> 2 1.681429 0.2413889
> 3 1.512447 0.2413889
>
> I think what's happening is that as.data.frame.svystat() method in the survey package ends up calling the wrong function to calculate the standard errors.  From the survey package:
>
> as.data.frame.svystat<-function(x,...){
>  rval<-data.frame(statistic=coef(x),SE=SE(x))
>  names(rval)[1]<-attr(x,"statistic")
>  if (!is.null(attr(x,"deff")))
>    rval<-cbind(rval,deff=deff(x))
>  rval
> }
>
> The relevant SE method seems to be:
>
> SE.svrepstat<-function(object,...){
>  if (is.list(object)){
>    object<-object[[1]]
>  }
>  vv<-as.matrix(attr(object,"var"))
>  if (!is.null(dim(object)) && length(object)==length(vv))
>    sqrt(vv)
>  else
>    sqrt(diag(vv))
> }

Mostly correct, but the relevant SE method is actually SE.default
survey:::SE.default
function (object, ...)
{
    sqrt(diag(vcov(object, ...)))
}

It can't be SE.svrepstat, because the class is wrong.

If you define
SE.svystat<-function(object,...){
   v<-vcov(object)
   if (!is.matrix(v) || NCOL(v)==1) sqrt(v) else sqrt(diag(v))
 }

it should work.

    -thomas

-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland