Hi Jarrod,

So, I tried following up on your suggestion to "create predictions for each
posterior sample of the parameters, and then average the predictions"...
and I tried to do this while getting predicted probabilities for every
combination of temp, contact, and judge, based on model MC3.

I've found this surprisingly labour-intensive. The eight lines of code I've
come up with to do this are below... but they seem quite cumbersome, and
I'm not even completely sure I've got this right.

Is there a better way of doing this? Using the posterior means of the
parameters, we could get predictions for a *single* combination of temp,
contact, and judge with just one line of code. Need it be so much more
complicated to get predictions for 36 combinations?

Many thanks in advance for any further suggestions you can provide.

- Malcolm


step1 <- expand.grid(temp=0:1, contact=0:1, judge=levels(wine$judge)) # 36
combinations
nrl <- MC3$Random$nrl
step2 <- cbind(1, step1[,-3], cbind(diag(nrl)[,rep(1:nrl,
MC3$Random$nfl)])[rep(1:nrl, rep(4, nrl)),])
step3 <- cbind(-Inf, 0, MC3$CP, Inf, MC3$Sol %*% t(as.matrix(step2)))
step4 <- t(apply(step3, 1, function(x) rep(x[1:6], each=36) -
x[7:length(x)]))
step5 <- pnorm(step4, mean=0, sd=sqrt(MC3$VCV[,"units"]))
step6 <- matrix(colMeans(step5), nrow=36)
step1 <- data.frame(step1, t(apply(step6, 1, diff)))



On 25 March 2014 08:56, Jarrod Hadfield <j.hadfield@ed.ac.uk> wrote:

> Hi Malcolm,
>
> I've commented in-line ...
>
>
> Quoting Malcolm Fairbrother <M.Fairbrother@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@ed.ac.uk>
>>> To: Shamil Sadigov <shamil@gmail.com>
>>> Cc: r-sig-mixed-models@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@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@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@ed.ac.uk>
>>> To: Jarrod Hadfield <j.hadfield@ed.ac.uk>
>>> Cc: r-sig-mixed-models@r-project.org
>>> Subject: Re: [R-sig-ME] Predicted values in MCMCglmm
>>>         family="threshold"
>>> Message-ID: <20140321142412.47024rvgurxlf12c@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@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@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.
>
>
>

	[[alternative HTML version deleted]]

