[R-sig-ME] trouble getting data scale predictions from a binomial MCMCglmm model
kari.ruohonen at utu.fi
Wed Mar 23 14:15:07 CET 2011
Hi Jarrod and thanks for a quick reply.
But the example I gave does not have random effects so the sum of
variance components is 1 as used in my calculation:
Iterations = 3001:12991
Thinning interval = 10
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
1 0 0 0
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
1 1 1 1 1
And with marginal=NULL the predictions are the same as without it in
this case. So, I am still not sure I got it. The predict.MCMCglmm seems
to do something differently than my direct model matrix & coefficients
approach in the example below.
On Wed, 2011-03-23 at 12:48 +0000, Jarrod Hadfield wrote:
> Dear Kari,
> The default in the predict function is to obtain the expectation after
> marginalising the random effects, hence you should be multiplying c2
> by the sum of the variance components (if the random structure is a
> simple grouping) if you want to obtain the same answer. Alternatively,
> you could specify marginal=NULL in the call to predict and this should
> be the same as your preds.data.
> See the posts for https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005575.html
> On 23 Mar 2011, at 12:19, Kari Ruohonen wrote:
> > Hi,
> > I have fitted a binomial model (logit link) with MCMCglmm by fixing
> > the
> > unit variance to 1 as discussed in the CourseNotes of MCMCglmm.
> > Now, I would like to get predictions on the data scale. I know there
> > is
> > a predict function that does this but I need to be able to predict
> > outside R. With the usual X%*%b where X is the model matrix and b the
> > vector of coefficients I get predictions on the link scale. The
> > predictions I get match those from the predict.MCMCglmm method with
> > the
> > type="terms" option. However, I cannot figure out how to get the
> > predictions on the data scale. CourseNotes discuss and give an example
> > of scaling the coefficients suitably with c2 <- (16*sqrt(3)/(15*pi))^2
> > and then dividing with (sqrt(1+c2*residual var)) to obtain estimates
> > under residual variance =0 but I cannot figure out how to apply this
> > to
> > the case of predictions.
> > Here is a small example based on one in the CourseNotes:
> > treatment<-gl(2,25)
> > y<-rbinom(50,1,c(.5,.001)[treatment])
> > data.bin<-data.frame(treatment=treatment,y=y)
> > prior.m2c.3<-list(R=list(V=1,fix=1),G=list(G1=list(V=1,nu=.002)))
> > m2c.3<-
> > MCMCglmm
> > (y~treatment,data=data.bin,family="categorical",prior=prior.m2c.
> > 3,verbose=FALSE)
> > # the next predictions correspond to those from predict.MCMCglmm
> > with "
> > # type="terms"
> > preds.link<-m2c.3$X%*%summary(m2c.3$Sol)$statistics[,1]
> > # but these do not correspond those with type="response"
> > preds.data<-plogis(preds.link[,1]/sqrt(1+c2*1))
> > I am obviously missing something and would be grateful to someone
> > helping me out.
> > Many thanks, Kari
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models