[R-sig-ME] predictions for MCMCglmm

[Joshua Wiley]

Hi Antonio,

You need to add both levels to the factor; however, both need not be
represented in your particular model matrix.  To do this, construct a
data frame where all the variables math your original data, but have
different values, a lá
data.frame(x = factor(1, levels = c(1, 2), labels = c("a", "b")))
b does not exist in that particular set of the data, but because the
factor has both levels, model matrix will work okay.

For a research project I was on, I developed a set of prediction
methods for MCMCglmm in the package postMCMCglmm that includes new
predictions and marginalizing random effects, with the ability to
handle new data.  That package is limited to the case I was working
with---a single ordinal outcome---but you may find the code and
documentation useful for examples.  Here is the main relevant code:

https://github.com/JWiley/postMCMCglmm/blob/master/R/prediction.R


I also defined extraction methods (e.g., fixef, ranef, etc.) elsewhere
in the package which you can lean on.  Actually, the code _may_ even
work for other families besides ordinal, as long as you stay on the
linear predictor metric---the back transformations definitely are not
implemented generally.

Cheers,

Josh



On Tue, Mar 19, 2013 at 11:35 AM, Antonio P. Ramos
<ramos.grad.student at gmail.com> wrote:
> sorry: the full code here
>> pred.data.pred <- data.frame(maternal_age_c=rep(20,25),wealth=rep("Lowest
> quintile",25),
> +
>  birth_year=1971:1995,birth_order=rep(1,25),
> +
>  sex=rep("Female",25),residence=rep("Rural",25),
> +                                    maternal_educ=rep("Primary",25))
>> # creating predictions "by hand"
>> X <- model.matrix(~ maternal_age_c + I(maternal_age_c^2)  +
> as.factor(birth_year) +
> +                   as.factor(birth_order) + residence + sex + wealth,
> data=pred.data.pred)
> Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
>   contrasts can be applied only to factors with 2 or more levels
>
>
>
>
> On Tue, Mar 19, 2013 at 11:35 AM, Antonio P. Ramos <
> ramos.grad.student at gmail.com> wrote:
>
>> Thanks for you reply Jarrod. The problem is that model.matrix doesn't
>> allow me to created the predication I need as factors cannot have just one
>> level. Any ways around data? thanks a bunch
>>
>>
>> pred.data.pred <- data.frame(maternal_age_c=rep(20,25),wealth=rep("Lowest
>> quintile",25),
>>
>>  birth_year=1971:1995,birth_order=rep(1,25),
>>
>>  sex=rep("Female",25),residence=rep("Rural",25),
>>                                    maternal_educ=rep("Primary",25))
>> # creating predictions "by hand"
>> X <- model.matrix(~ maternal_age_c + I(maternal_age_c^2)  +
>> as.factor(birth_year) +
>>                   as.factor(birth_order) + residence + sex + wealth,
>> data=pred.data.pred)
>>
>>
>>
>> On Tue, Mar 19, 2013 at 2:27 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk>wrote:
>>
>>> Hi Antonio,
>>>
>>> With (simple) random effects marginalised:
>>>
>>> X<-model.matrix(~ maternal_age_c + I(maternal_age_c^2)  +
>>> as.factor(birth_year) + residence + sex + wealth, data=newdata)
>>>
>>> V<-rowSums(glm.MC.2$VCV)
>>>
>>> beta<-glm.MC.2$Sol
>>>
>>> c2 <- (16 * sqrt(3)/(15 * pi))^2
>>>
>>> pred<-t(plogis(t(beta%*%t(X)/**sqrt(1+c2*V))))
>>>
>>> pred[i,j] is the prediction for the jth new data point for the ith MCMC
>>> sample. colSums(pred)  should be equivalent to the output from
>>> predict.MCMCglmm.
>>>
>>> Cheers,
>>>
>>> Jarrod
>>>
>>>
>>>
>>>
>>> Quoting "Antonio P. Ramos" <ramos.grad.student at gmail.com> on Mon, 18 Mar
>>> 2013 20:04:07 -0700:
>>>
>>>  Hi all.
>>>>
>>>> As far as I can tell newdata is still not implemented for this nice
>>>> package. Thus I wonder what would be the best way to get predictions "by
>>>> hand". My model is actually very simple. Still I need to marginalize the
>>>> random effects. Any hints? Thanks in advance, Antonio Pedro.
>>>>
>>>>
>>>> glm.MC.2 <- MCMCglmm(mortality.under.2 ~ maternal_age_c +
>>>> I(maternal_age_c^2)  +
>>>>                        as.factor(birth_year) + residence +
>>>>                        sex + wealth,
>>>>                      nitt=20000, thin=10, burnin=1000,
>>>>                      random= ~CASEID, prior=prior.2,data=egypt2,
>>>> family='categorical')
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> ______________________________**_________________
>>>> R-sig-mixed-models at r-project.**org <R-sig-mixed-models at r-project.org>mailing list
>>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-mixed-models<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.
>>>
>>>
>>>
>>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://joshuawiley.com/
Senior Analyst - Elkhart Group Ltd.
http://elkhartgroup.com