[R] survey package svystat objects from predict()

[Kieran Healy]

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))
}

Instead of returning sqrt(vv) on each element, it calculates sort(diag(vv)) instead. At least I think this is what's happening. 

I apologize in advance if all this is the result of some elementary error on my part. 

Thanks,

Kieran
  
-- 
Kieran Healy :: http://kieranhealy.org