Dear R users, Terry Theneau,

thank you very much for you answer. I'm running R 2.15.1 (32 bits) and
coxme 2.2-3. Here is a small R code which reproduces the problem (I fitted
a model with random effects whereas it is useless), it gives exactly the
same estimations of the variances than on my real data.

library(coxme)
set.seed(1989)
# parameters
N=500
beta1=0.5
beta2=-0.5
beta3=0.5
nb.groups=10
# variables
x=rbinom(N,1,0.5) # first covariate
y=rbinom(N,1,0.5) # second covariate
z=factor(rbinom(N,1,0.5)) # third covariate which will interacts with the
groups
id=factor(sample(1:nb.groups,N,T)) # groups
time=-log(runif(N))/exp(beta1*x+beta2*y+beta3*I(z==1)) # time of event or
censoring
eps=sample(0:1,N,T,c(0.3,0.7)) # event or censoring
data=data.frame(time,eps,x,y,z,id)

# preparing the coxme model
data$id_z=factor(paste(data$id,data$z,sep="-"))
names=paste(rep(levels(data$id),each=2),rep(levels(data$z),nb.groups),sep="-")
mat1=bdsmatrix(rep(c(1,0,0,0),nb.groups),blocksize=rep(2,nb.groups),dimnames=list(names,names))
mat2=bdsmatrix(rep(c(0,0,0,1),nb.groups),blocksize=rep(2,nb.groups),dimnames=list(names,names))
mat3=bdsmatrix(rep(c(0,1,1,0),nb.groups),blocksize=rep(2,nb.groups),dimnames=list(names,names))
varlist=coxmeMlist(list(mat1,mat2,mat3), rescale = F, pdcheck = F,
positive=F)

# models
fit.me=coxme(Surv(time,eps) ~ x + y + z + (1 |
id_z),varlist=varlist,data=data)
fit.me
fit.ph=coxph(Surv(time,eps) ~ x + y + z,data=data)
fit.ph

1-pchisq(-2*(fit.ph$loglik[2]-fit.me$loglik[2]),3)   # 0.9421373
1-pchisq(-2*(fit.ph$loglik[2]-fit.me$loglik[3]),3)   # 0.5929387

If I'm not wrong, the likelihood ratio tests above indicate adding the
random component is not necessary, which fits well with the way I simulated
the data.

Thank you again,

Hugo


2012/10/8 Terry Therneau <therneau@mayo.edu>

> You are right, those look suspicious.  What version of R and of the coxme
> package are you running?  Later version of coxme use multiple starting
> estimates due to precisely this kind of problem.
>   Also, when the true MLE is variance=0 the program purposely never quite
> gets there, in order to avoid log(0).  Compare the log-lik to a fixed
> effects model with those covariates.
>    I can't do more than guess without a reproducable example.
>
> Terry Therneau
>
>
> On 10/08/2012 05:00 AM, r-help-request@r-project.org wrote:
>
>> Dear R users,
>>
>> I'm using the function coxme of the package coxme in order to build Cox
>> models with complex random effects. Unfortunately, I sometimes get
>> surprising estimations of the variances of the random effects.
>>
>> I ran models with different fixed covariates but always with the same 3
>> random effects defined by the argument
>> varlist=coxmeMlist(list(mat1,**mat2,mat3), rescale = F, pdcheck = F,
>> positive=F). I get a few times exactly the same estimations of the
>> parameters of the random effects whereas the fixed effects of the models
>> are different:
>>
>> Random effects
>>   Group Variable Std Dev    Variance
>>   idp   Vmat.1   0.10000000 0.01000000
>>         Vmat.2   0.02236068 0.00050000
>>         Vmat.3   0.02449490 0.00060000
>>
>> The variances are round figures, so I have the feeling that the algorithm
>> didn't succeed in fitting the model.
>>
>> Has anyone ever faced to this problem?
>>
>> Thanks,
>>
>> Hugo
>>
>

	[[alternative HTML version deleted]]


