Hi,
I am modelling individual profiles consisting of 25 values each (diff)
evaluated at fractions 1-25 (frac) in two groups of patients (group). I
am fitting a semi-parametric model using a truncated polynomial basis
(Z.spline) as random effect. The model also has a random intercept
(FID).
fit12g <- tryCatch(lme(diff ~ frac * group, random = list(group =
pdIdent(~ Z.spline1-1), FID = pdSymm(~1)), method = "ML", data=g12,
na.action = na.omit, control = list(maxIter = 100, niterEM = 50)), error
= function(e) NA)
I wanted to check similar models without the random intercept:
fit12noIg <- tryCatch(lme(diff ~ frac * group, random = list(group =
pdIdent(~ Z.spline1-1)), method = "ML", data=g12, na.action = na.omit,
control = list(maxIter = 100, niterEM = 50)), error = function(e) NA)
or without the spline basis:
fit12noSg <- tryCatch(lme(diff ~ frac * group, random = list(FID =
pdSymm(~1)), method = "ML", data=g12, na.action = na.omit, control =
list(maxIter = 100, niterEM = 50)), error = function(e) NA)
Then I checked the log likelihood for these models:
logLik(fit12g)[1]
[1] -49818.38
> logLik(fit12noIg)[1]
[1] -51875.06
> logLik(fit12noSg)[1]
[1] -49817.43
I was very surprised to see that the likelihood for the full model
(logLik(fit12g) = -49818.38) was actually more negative than that of the
reduced model (logLik(fit12noSg) = -49817.43), resulting in a negative
log-likelihood ratio test
-2*(logLik(fit12noSg)[1]-logLik(fit12g)[1]) = -1.901062 !
Still, anova produced a positive likelihood test
> anova(fit12g, fit12noSg)
Model df AIC BIC logLik Test L.Ratio p-value
fit12g 1 7 99650.76 99704.55 -49818.38
fit12noSg 2 6 99646.86 99692.96 -49817.43 1 vs 2 1.901062 0.168
The likelihood for the model without the random intercept was higher
than that of the full model so the test was positive
-2*(logLik(fit12noIg)[1]-logLik(fit12g)[1]) = 4113.364
> anova(fit12g, fit12noIg)
Model df AIC BIC logLik Test L.Ratio p-value
fit12g 1 7 99650.76 99704.55 -49818.38
fit12noIg 2 6 103762.13 103808.23 -51875.06 1 vs 2 4113.364 <.0001
I know from the work of Crainesceau et al. that testing such
semi-parametric spline models is not straightforward and that the
distribution of the test statistic is a complicated and sometimes
unknown mixture of chi-square. Still, I would have expected the log
likelihood of the reduced model to be higher that that of the full model
and the test to be positive.
I wonder if I can trust the value of the likelhood of my reduced model
fit12noSg. Could it be the optimizer routine that is having problems
with model fit12noSg? If so, can someone indicate ways to check that?
Best,
Marie-Pierre
[[alternative HTML version deleted]]