R help - Feb 2012 - survey package svystat objects from predict()

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

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