I: Data structure issue (two-level nested hierarchy in survival data)
[This analysis was done in R v4.0.1, using the Coxme package for mixed effects Cox regression models]
I have collected data on ants trying to survive in different thermal conditions, and I hypothesize that the surface area of these ants as well as the chemical compounds they possess determine their survivability. These ants were assayed in groups of 20 (called “Samples”) replicated for each treatment, and these samples came from a few different nests. The response variable (survival, as a Surv(Time,Event)
object) is measured for each individual, but each individual’s body size is averaged at the level of “Sample”, and their chemical compound data is averaged at the level of “nest”. This hierarchical structure to the data looks like this:
> head(df, n=10)
## NEST SAMPLE TIME EVENT SA X1675 X1700 X1760
## 1 ABJ ABJ1 16 1 3.374059 0.3677422 2.46877 1.230814
## 2 ABJ ABJ1 16 1 3.374059 0.3677422 2.46877 1.230814
## 3 ABJ ABJ1 16 1 3.374059 0.3677422 2.46877 1.230814
## 4 ABJ ABJ1 20 1 3.374059 0.3677422 2.46877 1.230814
## 5 ABJ ABJ1 20 1 3.374059 0.3677422 2.46877 1.230814
## 6 ABJ ABJ2 16 1 3.184234 0.3677422 2.46877 1.230814
## 7 ABJ ABJ2 16 1 3.184234 0.3677422 2.46877 1.230814
## 8 ABJ ABJ2 16 1 3.184234 0.3677422 2.46877 1.230814
## 9 ABJ ABJ2 16 1 3.184234 0.3677422 2.46877 1.230814
## 10 ABJ ABJ2 20 1 3.184234 0.3677422 2.46877 1.230814
Here we can see one NEST (ABJ) that multiple SAMPLES were taken from (ABJ1,ABJ2,etc.), and while surface area (SA) is unique to each sample (3.37, 3.18, etc.), the chemical data (X1675, X1700, etc.) is unique to each nest. Because of this structure, I don’t think I can just do a usual Cox Proportional Hazard model (like coxph()
).
II: Comparing models, considering random effects
To deal with this nested hierarchy, I tried the following approach of mixed effects model comparison to see if adding my surface area and chemical data as predictor variables would change the AIC scores/log likelihood(loglik) of models allowing each NEST and SAMPLE to have a random intercept. I may not be explaining this well, but here is what that comparison looked like (method based on this YouTube video, and the models are using coxme()
):
> anova(m1,m2,m3)
## Analysis of Deviance Table
## Cox model: response is Surv(Time, Event)
## Model 1: ~1 + (1 | nest/Sample)
## Model 2: ~1 + sa + (1 | nest/Sample)
## Model 3: ~1 + sa + (1 | nest/Sample) + X1675 + X1880 + X2600 + X2700 + X2800 + X2900 + X3100
## loglik Chisq Df P(>|Chi|)
## 1 -9946.9
## 2 -9946.0 1.7814 1 0.182
## 3 -9930.6 30.8123 7 6.733e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model 1 consideres only random effects for each SAMPLE nested within NEST. Model 2 adds on the sample-level metric of surface area, and we see no change in loglik. Model 3 adds on the nest-level chemical data, and we see a significant change to loglik. So, my interpretation of this is that the chemical data, despite being a nest-level metric that is identical across many individuals, still adds meaningful “explanatory power” to the differences we see in survivability. The problem is, I would be satisfied here if this model considered all predictor variables I want to include. There are 7 more chemical compounds I want to include, but coxme()
can’t handle too many predictor variables while considering random intercepts, and I’ve read that, as a general rule of thumb, random effects + many predictor variables = less accurate coefficient estimates. Also, while these model comparisons tell us vaguely what predictor variables explain survivability, the hazard ratios that can be generated from Cox Proportional Hazard models (with no random effects) offer the best metrics for judging how each variable influences survivability.
III: The problem (am I following the rules?)
This is where I am least certain about my model formulation + choice. I figured, if a model containing random effects + chemical data is better than a model of just random effects, then the random effects are not critical to consider. Therefore, we can thus examine a model of just the chemical data, and focus on the hazard ratios to estimate how each chemical compound influences ant survivability. So, my questions, stated finally, would be: is it safe to ignore random effects because of my model comparisons so far, or do I need to stick to a model incorporating random effects, thus limiting the amount of predictor variables I want to include?
I can provide more information where needed. Thanks for the help!