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
______________________________________________
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.
