[R] Fitting a kind of Proportional Odds Modell using nlme, polr, lrm or ordgee

[Henning Henke]

Hello,

I'm trying to fit a special kind of proportional odds model from:

 Whitehead et al.  (2001).  Meta-analysis of ordinal outcome using 
individual patient data. Statistics in medicine 20: 2243-2260. (model 2)

The data are as follows:


library(nlme)
library(geepack)
library(Design)
library(MASS)

options(contrasts=c("contr.SAS","contr.poly"))
counts  <- 
c(2,22,54,29,3,4,23,45,22,2,1,22,35,11,3,14,119,180,54,6,7,16,17, 
10,3,13,20,24,10,1,8,24,73,52,13,21,106,175,62,17,2,13,18,7,1,3,14,19,3,0)
category <- c(rep(1:5,10))
treatment <- (rep(c(0,0,0,0,0,1,1,1,1,1,),5))
study <- gl(5,10)
percoftotal <- counts/sum(counts)

cout <- data.frame(study,treatment,category,counts,percoftotal)





The data are from five controlled clinical trials (study). The patients 
were given placebo (treatment=0) or treatment (treatment=1) and had a 
response in one of five categories.

Now I want to fit the model given by Whitehead et al. This assumes 
proportional odds between treatments, but stratifies by study. This 
means that the cut-off points associated with the distribution of the 
underlying latent variable for determining the response category are 
allowed to vary from study to study but are the same for both treatment 
groups within a study.
It is given by

log(Q_{ijk}/(1-Q_{ijk}) =  \alpha_{ik} + \beta*x_{ij} (k=1,...,m-1)

This model can be considered as arising from a latent continuos variable 
. Assume that the response of the j-th subject in study i is truly equal 
to G_{ij} although this latent reponse will never be observerd.
G_{ij} has a logistic distribution with:

 Q_{ijk} = P(G_{ij} <= \alpha_{ik}) = 
1/(1+exp(-(\alpha_{ik}+\beta*x_{ij})      (k=1,...,m-1)

in which Q_{ijk} is the probabilty of having a response for the j-th 
subject in category k or better that means p_{ij1}+....+p_{ijk}= Q_{ijk} 
and Q_{ijm}=1, Q_{ij1}=1.



i   = 1,...,r  (r=5) #number of studies
j   = 1,....,n_{i} #subject j from study i
k  = 1,...,m (m=5)#Category

My problem is how to fit this model in R although I have a SAS code from 
Whitehead which is availabe at  
http://www.rdg.ac.uk/mps/mps_home/misc/publications.htm . There it is 
fit by using Proc NLMIXED.  The result for beta is there:

The problem is  that alpha_{ik} represents the k_th intercept for the 
i_th study. A "normal" Proportional-Odds model
has only alpha_k which represents the k-th intercept without any 
associaton to the study.

I tried to fit it by using the functions of Pinheiro/ Bates nlme, nls, 
gnls;
   polr  by Venables/Ripley ;
   lrm  by Harrell
   ordgee by Yan;

but I came to no conclusion.


I  would be very happy if anyone could help me.

Thank you VERY VERY much in advance.

Kindly regards,

Henning Henke