[R-sig-ME] R-sig-mixed-models Digest, Vol 87, Issue 34
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Mar 25 09:56:20 CET 2014
Hi Malcolm,
I've commented in-line ...
Quoting Malcolm Fairbrother <M.Fairbrother at bristol.ac.uk> on Tue, 25
Mar 2014 02:08:39 +0000:
> Dear Jarrod,
>
> Following on from your response to Shamil, and your comments about the new
> features in MCMCglmm version 2.18, can I please ask you to confirm a couple
> details about the extraction of predicted probabilities from fitted
> threshold models?
>
> As I did once before when asking a related question (
> http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/9424), I'll refer
> to the "wine" dataset in the ordinal package. My questions appear below.
>
> library(ordinal)
> library(MCMCglmm)
>
> data(wine)
> str(wine)
> summary(clmod <- clmm2(rating ~ temp + contact, random=judge, data=wine,
> Hess=TRUE, nAGQ=10, link="probit"))
>
> # to get predicted probabilities for each possible response, for one
> combination of covariates:
> diff(pnorm(c(-Inf,clmod$Theta, Inf) - clmod$beta %*% c(0,1), 0,
> sqrt(1+clmod$stDev))) # for tempcold and contactyes
>
> # now fit the same model, using MCMCglmm, two ways:
> prior1 <- list(R = list(V = 1, fix=1), G = list(G1 = list(V = 1, nu =
> 0.002)))
> MC1 <- MCMCglmm(rating ~ temp + contact, random=~judge, data=wine,
> family="ordinal", prior=prior1, nitt=130000, thin=100, verbose=F)
> MC2 <- MCMCglmm(rating ~ temp + contact, random=~judge, data=wine,
> family="threshold", prior=prior1, nitt=130000, thin=100, verbose=F)
>
> # and get predicted probabilities, using the results from MCMCglmm,
> marginalising the random effects:
> diff(pnorm(c(-Inf, 0, colMeans(MC1$CP), Inf)-colMeans(MC1$Sol) %*%
> c(1,0,1), 0, sqrt(1+sum(colMeans(MC1$VCV))))) # for tempcold and contactyes
> diff(pnorm(c(-Inf, 0, colMeans(MC2$CP), Inf)-colMeans(MC2$Sol) %*%
> c(1,0,1), 0, sqrt(1+sum(colMeans(MC2$VCV))))) # for tempcold and contactyes
>
> # Q1: They're not the same... and I'm not sure what precisely would make
> them the same. (The ones from MC1 match the ones from clmod, but those from
> MC2 don't.)
The threshold variance should not include the probit variance of 1.
The probit variance is what is specified as the units variance in
family="threshold":
diff(pnorm(c(-Inf, 0, colMeans(MC2$CP), Inf)-colMeans(MC2$Sol) %*%
c(1,0,1), 0,sqrt(sum(colMeans(MC2$VCV))))) # for tempcold and contactyes
PS your code is a neater way of getting the probabilities than mine,
but perhaps it is easier to understand as:
diff(pnorm(c(-Inf, 0, colMeans(MC2$CP), Inf), colMeans(MC2$Sol) %*%
c(1,0,1),sqrt(sum(colMeans(MC2$VCV))))) # for tempcold and contactyes
>
> # finally, fit a similar model but with a random slope for the binary
> covariate "contact":
> prior2 <- list(R = list(V = 1, fix=1), G = list(G1 = list(V = diag(2), nu =
> 2.002))) # I believe this is an uninformative prior
> MC3 <- MCMCglmm(rating ~ temp + contact, random=~ us(1+contact):judge,
> data=wine, family="threshold", prior=prior2, nitt=130000, thin=100,
> verbose=F, pr=T)
>
> # Q2: How would one get predicted probabilities here, (a) marginalising the
> random effects, and (b) for a particular judge (say #1)?
>
> # for (b), would it be something like this (again for tempcold and
> contactyes)?
> diff(pnorm(c(-Inf, 0, colMeans(MC3$CP), Inf)-colMeans(MC3$Sol)[c(1:4,13)]
> %*% c(1,0,1,1,1), 0, sqrt(1+sum(colMeans(MC1$VCV)))))
>
> # but I'm guessing the "sqrt(1+sum(colMeans(MC1$VCV)))" bit is wrong... I'm
> not sure what you mean by "the sum of the variance components associated
> with that term"
In this model we have to decide when predicting whether we want to
predict for contact=yes or contact=no. We could marginalise this too,
but this would require knowing something about the distribution of
contact in the population, and this isn't modelled in this instance.
Lets say we want to predict for contact=yes as you have attempted. The
prediction for judge 1 is:
diff(pnorm(c(-Inf, 0, colMeans(MC3$CP),
Inf)-colMeans(MC3$Sol)[c(1:4,13)] %*% c(1,0,1,1,1), 0,
sqrt(mean(MC3$VCV[,"units"]))))
because the only random terms left to marginalise are the residuals.
Marginalising the judge effects is a bit trickier in this instance
because the judge effects for contact=yes are the sum of two terms and
we need to know the variance of that sum. It is the sum of the
variances for each effect plus twice the covariance between them.
With the units variance too this is just equal to the sum of all
(co)variances (because MCMCglmm returns matrices so gives the
covariance twice):
diff(pnorm(c(-Inf, 0, colMeans(MC3$CP), Inf)-colMeans(MC3$Sol)[c(1:3)]
%*% c(1,0,1), 0, sqrt(sum(colMeans(MC3$VCV)))))
In some ways it would be easier (and more instructive in this example)
to have fitted us(contact) rather than us(1+contact). With the former
parameterisation the first set of random effects is associated with
contact=no (as before) but the second set of random effects is
associated directly with contact=yes (unlike before). Marginalising
the random effects for contact=yes is then:
diff(pnorm(c(-Inf, 0, colMeans(MC3$CP), Inf)-colMeans(MC3$Sol)[c(1:3)]
%*% c(1,0,1), 0, sqrt(sum(colMeans(MC3$VCV[,c("units",
"yes:yes.judge")])))))
and for contact=no is
diff(pnorm(c(-Inf, 0, colMeans(MC3$CP), Inf)-colMeans(MC3$Sol)[c(1:3)]
%*% c(1,0,1), 0, sqrt(sum(colMeans(MC3$VCV[,c("units",
"no:no.judge")])))))
The covariances are not required under this formulation for these predictions.
Note that these predictions are for based on the posterior mean of the
parameters. It would be more satisfying to create predictions for each
posterior sample of the parameters, and then average the predictions.
This would give the posterior mean predictions rather than the
predictions conditional on the posterior mean of the parameters.
However, in many cases the two things will be close.
Cheers,
Jarrod
>
> As always, any assistance would be much, much appreciated.
>
> - Malcolm
>
>
>
> Date: Fri, 21 Mar 2014 14:14:50 +0000
>> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
>> To: Shamil Sadigov <shamil at gmail.com>
>> Cc: r-sig-mixed-models at r-project.org
>> Subject: Re: [R-sig-ME] Predicted values in MCMCglmm
>> family="threshold"
>>
>> Hi,
>>
>> have
>>
>> cp<-c(-Inf, 0, cp.est, Inf)
>>
>> where cp.est are the estimated cutpoints (if there are any - with 2
>> categories there are none).
>>
>> Have linear predictor nu = xb or nu=xb+zu. If the former (and there
>> are random effects) then have v the sum of the variance components
>> associated with that term, and if the latter have v as the units
>> variance associated with that term.
>>
>> Have obs<-1:k where k is the number of categories (2+the number of
>> estimated cutpoints) and the probability of falling into a category
>> conditional on nu and v is:
>>
>> pnorm(cp[obs+1], nu , sqrt(v)) ? pnorm(cp[obs], nu, sqrt(v))
>>
>> for family=threshold, and
>>
>> pnorm(cp[obs+1], nu , sqrt(v+1)) ? pnorm(cp[obs], nu, sqrt(v+1))
>>
>> for ordinal.
>>
>> For example,
>>
>> cp.est<-1
>> cp<-c(-Inf, 0, cp.est, Inf)
>> k<-2+length(cp.est)
>> obs<-1:k
>> nu<--1
>> v<-2
>> pnorm(cp[obs+1], nu , sqrt(v))-pnorm(cp[obs], nu, sqrt(v))
>>
>> Jarrod
>>
>>
>> Quoting Shamil Sadigov <shamil at gmail.com> on Fri, 21 Mar 2014 14:56:25
>> +0200:
>>
>> > Hi Jarrod,
>> >
>> > I am using the new family="threshold" in MCMCglmm version 2.18 with a
>> > 5-variate ordered response. I would like to obtain the predicted
>> responses
>> > for on the original ordinal scale, but I am not sure how to do so for
>> > either "ordinal" or the "threshold" family.
>> >
>> > 1. For family="threshold" the posterior predicted probabilities are :
>> >
>> > post.pred[, keep] <- pnorm(post.pred[, keep], 0,
>> > sqrt(postvar[, keep]))
>> >
>> > How can I classify these probabilities into the original ordinal scale?
>> >
>> >
>> > 2. I can see that for family="ordinal", cut points (CP) are used in
>> > predict.MCMCglmm():
>> >
>> > for (i in 2:(dim(CP)[2] - 1)) {
>> > q <- q + (pnorm(CP[, i + 1] - post.pred[, keep], 0,
>> > sqrt(postvar[, keep] + 1)) - pnorm(CP[, i] - post.pred[, keep], 0,
>> > sqrt(postvar[, keep] + 1))) * (i - 1)
>> > }
>> > Are the thresholds and the posterior predictive values (using type =
>> > "terms") on the linear (latent variable) scale?
>> >
>> > What would be the interpretation of the predicted values obtained from
>> > using type= "response" with family = "ordinal"? (All 5 ordinal responses
>> > are coded 1-3, and the predicted values from predict.MCMCglmm are real
>> > numbers between 0-6.)
>> >
>> > Regards,
>> > Shamil.
>> >
>> > [[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > 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.
>>
>>
>>
>> ------------------------------
>>
>> Message: 2
>> Date: Fri, 21 Mar 2014 14:24:12 +0000
>> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
>> To: Jarrod Hadfield <j.hadfield at ed.ac.uk>
>> Cc: r-sig-mixed-models at r-project.org
>> Subject: Re: [R-sig-ME] Predicted values in MCMCglmm
>> family="threshold"
>> Message-ID: <20140321142412.47024rvgurxlf12c at www.staffmail.ed.ac.uk>
>> Content-Type: text/plain; charset=ISO-8859-1; DelSp="Yes";
>> format="flowed"
>>
>> Hi,
>>
>> ? should be - (and microsoft should be *)
>>
>> Jarrod
>>
>>
>> Quoting Jarrod Hadfield <j.hadfield at ed.ac.uk> on Fri, 21 Mar 2014
>> 14:14:50 +0000:
>>
>> > Hi,
>> >
>> > have
>> >
>> > cp<-c(-Inf, 0, cp.est, Inf)
>> >
>> > where cp.est are the estimated cutpoints (if there are any - with 2
>> > categories there are none).
>> >
>> > Have linear predictor nu = xb or nu=xb+zu. If the former (and there
>> > are random effects) then have v the sum of the variance components
>> > associated with that term, and if the latter have v as the units
>> > variance associated with that term.
>> >
>> > Have obs<-1:k where k is the number of categories (2+the number of
>> > estimated cutpoints) and the probability of falling into a category
>> > conditional on nu and v is:
>> >
>> > pnorm(cp[obs+1], nu , sqrt(v)) ? pnorm(cp[obs], nu, sqrt(v))
>> >
>> > for family=threshold, and
>> >
>> > pnorm(cp[obs+1], nu , sqrt(v+1)) ? pnorm(cp[obs], nu, sqrt(v+1))
>> >
>> > for ordinal.
>> >
>> > For example,
>> >
>> > cp.est<-1
>> > cp<-c(-Inf, 0, cp.est, Inf)
>> > k<-2+length(cp.est)
>> > obs<-1:k
>> > nu<--1
>> > v<-2
>> > pnorm(cp[obs+1], nu , sqrt(v))-pnorm(cp[obs], nu, sqrt(v))
>> >
>> > Jarrod
>> >
>> >
>> > Quoting Shamil Sadigov <shamil at gmail.com> on Fri, 21 Mar 2014 14:56:25
>> +0200:
>> >
>> >> Hi Jarrod,
>> >>
>> >> I am using the new family="threshold" in MCMCglmm version 2.18 with a
>> >> 5-variate ordered response. I would like to obtain the predicted
>> responses
>> >> for on the original ordinal scale, but I am not sure how to do so for
>> >> either "ordinal" or the "threshold" family.
>> >>
>> >> 1. For family="threshold" the posterior predicted probabilities are :
>> >>
>> >> post.pred[, keep] <- pnorm(post.pred[, keep], 0,
>> >> sqrt(postvar[, keep]))
>> >>
>> >> How can I classify these probabilities into the original ordinal scale?
>> >>
>> >>
>> >> 2. I can see that for family="ordinal", cut points (CP) are used in
>> >> predict.MCMCglmm():
>> >>
>> >> for (i in 2:(dim(CP)[2] - 1)) {
>> >> q <- q + (pnorm(CP[, i + 1] - post.pred[, keep], 0,
>> >> sqrt(postvar[, keep] + 1)) - pnorm(CP[, i] - post.pred[, keep], 0,
>> >> sqrt(postvar[, keep] + 1))) * (i - 1)
>> >> }
>> >> Are the thresholds and the posterior predictive values (using type =
>> >> "terms") on the linear (latent variable) scale?
>> >>
>> >> What would be the interpretation of the predicted values obtained from
>> >> using type= "response" with family = "ordinal"? (All 5 ordinal responses
>> >> are coded 1-3, and the predicted values from predict.MCMCglmm are real
>> >> numbers between 0-6.)
>> >>
>> >> Regards,
>> >> Shamil.
>>
>>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list