[R] Question regarding significance of a covariate in a coxme survival model
Christopher David Desjardins
I am running a Cox Mixed Effects Hazard model using the library coxme. I
am trying to model time to onset (age_sym1) of thought problems (e.g.
hearing voices) (sym1). As I have siblings in my dataset, I have
decided to account for this by including a random effect for family
(famid). My covariate of interest is Mother's diagnosis where a 0 is
bipolar, 1 is control, and 2 is major depression. I am trying to fit
the following model.
thorp1 <- coxme(Surv(age_sym1, sym1) ~ lifedxm + (1 | famid), data =
bip.surv)
Which provides the following output:
-------------------------------------------------
> thorp1
Cox mixed-effects model fit by maximum likelihood
Data: bip.surv
events, n = 99, 189 (3 observations deleted due to missingness)
Iterations= 10 54
NULL Integrated Penalized
Log-likelihood -479.0372 -467.3549 -435.2096
Chisq df p AIC BIC
Integrated loglik 23.36 3.00 3.3897e-05 17.36 9.58
Penalized loglik 87.66 47.27 3.2374e-04 -6.88 -129.54
Model: Surv(age_sym1, sym1) ~ lifedxm + (1 | famid)
Fixed coefficients
coef exp(coef) se(coef) z p
lifedxmCONTROL -1.0576257 0.3472794 0.3872527 -2.73 0.0063
lifedxmMAJOR -0.6329957 0.5309987 0.3460847 -1.83 0.0670
Random effects
Group Variable Std Dev Variance
famid Intercept 0.9877770 0.9757033
--------------------------------------------------
So it appears that there is a difference in hazard rates associated with
Mother's diagnosis between BIPOLAR and CONTROL. Just to be safe, I fit a
model without this covariate present and decided to compare the models
with AIC and see if fit decreased with this covariate not in the model.
----------------------------------------------------------
thorp1.n <- coxme(Surv(age_sym1, sym1) ~ (1 | famid), data = bip.surv)
> thorp1.n
Cox mixed-effects model fit by maximum likelihood
Data: bip.surv
events, n = 99, 189 (3 observations deleted due to missingness)
Iterations= 10 54
NULL Integrated Penalized
Log-likelihood -479.0372 -471.3333 -436.0478
Chisq df p AIC BIC
Integrated loglik 15.41 1.00 0.00008663 13.41 10.81
Penalized loglik 85.98 49.48 0.00099915 -12.97 -141.37
Model: Surv(age_sym1, sym1) ~ (1 | famid)
Random effects
Group Variable Std Dev Variance
famid Intercept 1.050949 1.104495
------------------------------------------------------------
The problem is that the AIC for the model w/o lifedxm is 13.41 and the
model w/ lifedxm is 17.36. So fit actually improved without lifedxm in
the model but really the fit is indistinguishable if I use Burnham &
Anderson, 2002's criteria.
So my conundrum is this - The AIC and the p-values appear to disagree.
The p-value would suggest that lifedxm should be included in the model
and the AIC says it doesn't improve fit. In general, I tend to favor the
AIC (I usually work with large sample sizes and multilevel models) but I
am just beginning with survival models and I am curious if a survival
model guru might suggest whether lifedxm needs to be in the model or not
based on my results here? Would people generally use the AIC in this
situation? Also, I can't run anova() on this models because of the
random effect.
I am happy to provide the data if necessary.
Please cc me on a reply.
Thanks,
Chris
______________________________________________
R-h...@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.
Teresa Iglesias
Hi Chris,
Did you get an answer to why the p-value and AIC score disagreed?
I am getting the same results with my own data. They're not small
disagreements either. The AIC score is telling me the opposite of
what the p-value and the parameter coef are saying.
The p-value and the coef for the predictor variable are in agreement.
I've also noticed in some published papers with tables containing
p-values and AIC scores that the chisq p-value and AIC score
are in opposition too. This makes me think I'm missing something
and that this all actually makes sense.
However given that AIC = − 2 (log L) + 2K (where K is the number of parameters)
and seeing as how the log-likelihood scores below are in the hundreds, shouldn't
the AIC score be in the hundreds as well??
--------------------------------------
model0 (intercept only model)
NULL Integrated Penalized
Log-likelihood -119.8470 -117.9749 -113.1215
Chisq df p AIC BIC
Integrated loglik 3.74 1.00 0.052989 1.74 0.08
Penalized loglik 13.45 7.06 0.063794 -0.68 -12.43
Random effects
Group Variable Std Dev Variance
Site Intercept 0.6399200 0.4094977
_____
model1 (before vs after)
NULL Integrated Penalized
Log-likelihood -119.8470 -112.1598 -108.1663
Chisq df p AIC BIC
Integrated loglik 15.37 2.00 0.00045863 11.37 8.05
Penalized loglik 23.36 7.06 0.00153710 9.25 -2.49
Fixed coefficients
Coef exp(coef) se(coef) z p
Time -1.329997 0.2644779 0.4009229 -3.32 0.00091
Random effects
Group Variable Std Dev Variance
Site Intercept 0.5625206 0.3164294
______
model2 (stim1 vs stim2)
NULL Integrated Penalized
Log-likelihood -119.8470 -117.8677 -113.167
Chisq df p AIC BIC
Integrated loglik 3.96 2.00 0.138160 -0.04 -3.37
Penalized loglik 13.36 7.83 0.093314 -2.31 -15.34
Fixed coefficients
coef exp(coef) se(coef) z p
stim 0.1621367 1.176021 0.3534788 0.46 0.65
Random effects
Group Variable Std Dev Variance
Site Intercept 0.6262600 0.3922016
----------------------------------------------
If you didn't get an answer, hopefully, someone that understands all this will
reply for both of us.
Thanks,
-Teresa
David Winsemius
That is different than my understanding of AIC. I thought that the AIC
and BIC both took as input the difference in -2LL and then adjusted
those differences for the differences in number of degrees of freedom.
That perspective would inform one that the absolute value of the
deviance ( = -2LL) is immaterial and only differences are subject to
interpretation. The p-values are calculated as Wald tests, which are
basically first order approximations and have lower credence than
model level comparisons. The OP also seemed to have ignored the fact
that the penalized estimates of AIC and BIC went in the opposite
direction from what he might have hoped.
Yes. it would be good to get more authoritative opinion. Perhaps my
standing as a non-statistician will prompt a real statistician to
weigh in and correct any misapprehensions I have been laboring under.
--
David Winsemius, MD
West Hartford, CT
Teresa Iglesias
> Hi David,
Maybe I misunderstand your reply about AIC (and BIC) but my understanding is
that the AIC "score" as an absolute value is meaningless but the "weights"
calculated from the models you want to compare is what matters. So the AIC
is used to calculate the relative Akaike weight among all the models you
wish to compare where the weights add to 1.
As far as I can tell, coxme gives the model's AIC "score" but is not
comparing any models yet. Perhaps that is up to the user to calculate
weights given the scores and figure out which model performs best.
The output from coxme above (if I understand correctly) is the AIC score for
each model (evaluated independent of other models) and that the AIC score
(is? should be?) calculated using the LL given as output. If this is true
then I don't see how one could get those AIC scores from those LLs.
Integrated Penalized
Log-likelihood -117.8677 -113.167
AIC BIC
Integrated loglik -0.04 -3.37
Penalized loglik -2.31 -15.34
-2(-117.86)+2(2) ≠ -0.04
--
Teresa Iglesias
Davis, CA
[[alternative HTML version deleted]]
Kaigang Li
[[alternative HTML version deleted]]
Johnson, Cedrick W.
C. Peng
The likelihood ratio test is more reliable when one model is nested in the
other. This true for your case.
AIC/SBC are usually used when two models are in a hiearchical structure.
Please also note that any decision made made based on AIC/SBC scores are
very subjective since no sampling distribution can be used to make a
"rigorous" decision. Regarding the magnitutes between the loglikelihood and
AIC/SBC, I would say the author must used a modified version in coxme()
since several different modified AIC/SBC scores are running in practice.
My suggestion would be to use LR test for your case:
For the integrated likelihhod:
LL.small.model = - 467.3549 (including lifedxm)
LL.large.model = - 471.3333 (excluding lifedxm)
DF.diff = 3 - 1 = 2
LR: -2*(- 471.3333 + 467.3549) = 7.9568
p-value: 1-pchisq(7.9568,2) = 0.01871556
For the penalized likelihhod:
LPL.small.model = -435.2096 (including lifedxm)
LPL.large.model = -436.0478 (excluding lifedxm)
DF.diff = 3 - 1 = 2
PLR: -2*(- 436.0478 + 435.2096 ) = 1.6764
p-value: 1-pchisq(1.6764,2) = 0.4324883
Two different likehood methods produce different results, which one you
should use depends
on which likelihood makes more sense to you (or which likehood is better).
HTH
--
View this message in context: http://r.789695.n4.nabble.com/Question-regarding-significance-of-a-covariate-in-a-coxme-survival-model-tp2313880p2399114.html
Sent from the R help mailing list archive at Nabble.com.
Cheng Peng
The likelihood ratio test is more reliable when one model is nested in the
other. This true for your case.
AIC/SBC are usually used when two models are in a hiearchical structure.
Please also note that any decision made
made based on AIC/SBC scores are very subjective since no sampling
distribution can be used to make a "rigorous" decision.
the author must used a modified version in coxme()
since several different modified AIC/SBC scores are running in practice.
My suggestion would be to use LR test for your case:
For the integrated likelihhod:
LL.small.model = - 467.3549 (including lifedxm)
LL.large.model = - 471.3333 (excluding lifedxm)
DF.diff = 3 - 1 = 2
LR: -2*(- 471.3333 + 467.3549) = 7.9568
p-value: 1-pchisq(7.9568,2) = 0.01871556
For the penalized likelihhod:
LPL.small.model = -435.2096 (including lifedxm)
LPL.large.model = -436.0478 (excluding lifedxm)
DF.diff = 3 - 1 = 2
PLR: -2*(- 436.0478 + 435.2096 ) = 1.6764
p-value: 1-pchisq(1.6764,2) = 0.4324883
Two different likehood methods produce different results, which one you
should use depends
on which likelihood makes more sense to you (or which likehood is better).
HTH
--
View this message in context: http://r.789695.n4.nabble.com/Question-regarding-significance-of-a-covariate-in-a-coxme-survival-model-tp2313880p2399090.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
Cheng Peng
My suggestion:
If compare model 1 and model 2 with model 0 respectively, the (penalized)
likelihood ratio test is valid.
IF you compare model 2 with model 3, the (penalized) likelihood ratio test
is invalid. You may want to use AIC/SBC to make a subjective decision.
--
View this message in context: http://r.789695.n4.nabble.com/Question-regarding-significance-of-a-covariate-in-a-coxme-survival-model-tp2313880p2399095.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
C. Peng
My suggestion for Teresa:
If compare model 1 and model 2 with model 0 respectively, the (penalized)
likelihood ratio test is valid.
IF you compare model 2 with model 3, the (penalized) likelihood ratio test
is invalid. You may want to use AIC/SBC to make a subjective decision.
--
View this message in context: http://r.789695.n4.nabble.com/Question-regarding-significance-of-a-covariate-in-a-coxme-survival-model-tp2313880p2399116.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
Teresa Iglesias
David Winsemius wrote:
>
> That is different than my understanding of AIC. I thought that the AIC
> and BIC both took as input the difference in -2LL and then adjusted
> those differences for the differences in number of degrees of freedom.
>
>
David! Your words make sense to me now. Sorry for the lapse.
A very smart professor took the time out to school me. I misunderstood the
output from coxme. I see now that it's giving the LL for the NULL model and
for the model I have specified and the AIC output is the difference between
the full model and the NULL. So the numbers all make sense to me and in fact
the p-values are "in agreement" in that they decrease as the specified model
is an improvement over the NULL. I am using only the AIC values and akaike
weights in my analyses so the p-values are not my basis for reaching a
conclusion. I was distressed over the seeming disagreement between the AIC,
the p-values and my results using lmer (ignoring that the data were
censored), and bar graphs illustrating a clear difference when coxme was
telling me "otherwise". To spell it out for others like me that need to see
the numbers add up...
Given this output from coxme:
-------------------------------------------------------
NULL Integrated Penalized
Log-likelihood -119.8470 -112.1598 -108.1663
Chisq df p AIC BIC
Integrated loglik 15.37 2.00 0.00045863 11.37 8.05
Penalized loglik 23.36 7.06 0.00153710 9.25 -2.49
--------------------------------------------------------
-2(LL) + (2*df) = AIC
NULL: -2(-119.8470) = 239.694
Integrated: -2(-112.1598) + (2 * 2) = 228.3196
Penalized: -2(-108.1663) + (2 * 7.06) = 230.4526
subtract the integrated model's AIC from the NULL model's AIC, you
get the stated AIC for the integrated model in the output (same for
penalized model).
239.694 - 228.3196 = 11.3744
So the larger (positive range) the AIC, in the coxme output, the better that
model does compared to the NULL model. Incidentally, we see that the p-value
decreases with an increase in the coxme AIC and so there is no
"disagreement".
Thank you very smart professor!
-Teresa Iglesias
Davis, Ca
______________________________________________
R-h...@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.
--
View this message in context: http://r.789695.n4.nabble.com/Question-regarding-significance-of-a-covariate-in-a-coxme-survival-model-tp2313880p2527739.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________