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 <thern...@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-requ...@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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
