[R-sig-ME] trouble getting data scale predictions from a binomial MCMCglmm model

[Jarrod Hadfield]

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

Cheers,

Jarrod



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
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.