R-sig-mixed-models Digest, Vol 87, Issue 34

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)