On Tue, Feb 14, 2012 at 11:45 AM, Kieran Healy <kjhe...@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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
