0
$\begingroup$

In an ecological experiment I have recorded survival of individuals in response to two different food items (treatment). Individuals belong to different families (random term) and are measured at 5 timepoints. From just looking at the raw data, it becomes clear that food as a strong effect on survival across all families. When I report the results, I would like to show just that, the raw data, overlayed by statistical estimates for survival across all families per diet per timepoint.

However, I can't seem to figure out how I can get survial estimates that span across all families and change over time, while at the same time families are a random effect (I want to use the estimates for the random effects in another analysis).

I hope what I want becomes clear with the first example using lmer: I know this is not the right statistical test for survival (is it though??), but I was not able to extract estimates like this from the other tests I tried (coxph and coxme, below).

I don't get why the estimates from coxph are split up by individual - if I specify frailty by family, I should get only family level estimates, no?

The estimates from coxme don't change with time, which cannot be right - is it because of how I specify the random effect? Should individual be nested in family to specifiy the repeatedness?

library(data.table)
library(survival)
library(coxme)
library(lme4)
library(ggplot2)
library(sjPlot)

sub = structure(list(Family = c("fam086", "fam086", "fam086", "fam086", 
"fam086", "fam086", "fam086", "fam086", "fam086", "fam086", "fam086", 
"fam086", "fam086", "fam086", "fam086", "fam086", "fam086", "fam086", 
"fam086", "fam086", "fam086", "fam086", "fam086", "fam086", "fam086", 
"fam086", "fam086", "fam086", "fam086", "fam086", "fam090", "fam090", 
"fam090", "fam090", "fam090", "fam090", "fam090", "fam090", "fam090", 
"fam090", "fam090", "fam090", "fam090", "fam090", "fam090", "fam090", 
"fam090", "fam090", "fam090", "fam090", "fam090", "fam090", "fam090", 
"fam090", "fam090", "fam090", "fam090", "fam090", "fam090", "fam090", 
"fam118", "fam118", "fam118", "fam118", "fam118", "fam118", "fam118", 
"fam118", "fam118", "fam118", "fam118", "fam118", "fam118", "fam118", 
"fam118", "fam118", "fam118", "fam118", "fam118", "fam118", "fam118", 
"fam118", "fam118", "fam118", "fam118", "fam118", "fam118", "fam118", 
"fam118", "fam118"), Individual = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 
6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 
6L, 6L, 6L), Food = c("H", "H", "H", "H", "H", "L", "L", "L", 
"L", "L", "H", "H", "H", "H", "H", "L", "L", "L", "L", "L", "H", 
"H", "H", "H", "H", "L", "L", "L", "L", "L", "H", "H", "H", "H", 
"H", "L", "L", "L", "L", "L", "H", "H", "H", "H", "H", "L", "L", 
"L", "L", "L", "H", "H", "H", "H", "H", "L", "L", "L", "L", "L", 
"H", "H", "H", "H", "H", "L", "L", "L", "L", "L", "H", "H", "H", 
"H", "H", "L", "L", "L", "L", "L", "H", "H", "H", "H", "H", "L", 
"L", "L", "L", "L"), Status = c(FALSE, FALSE, FALSE, TRUE, TRUE, 
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, 
TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, 
TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, 
TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, 
TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, 
TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, 
TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, 
TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, 
TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE), Time = c(1, 2, 
3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 
4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 
5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 
1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 
2, 3, 4, 5)), class = "data.frame", row.names = c(NA, -90L))                                                                                                                  
surv.obj = with(sub, Surv(time=Time, event=Status))

setDT(sub)

# random effects with lme4
sub[,Status_num:=1-as.numeric(Status)]
mod3 = glmer(Status_num ~ Food * Time + (1|Family) , data = sub, family="binomial")
plot_model(mod3, type = "pred", terms = c("Time","Food"))

# random effects with coxph
mod1 = coxph(surv.obj ~ Food + frailty(Family) , data = sub)
sub$pred1 = exp(-predict(mod1, type="expected"))
qplot(Time, pred1, colour=Food, dat=sub)

# random effects with copxme
mod2 = coxme(surv.obj ~ Food + (1|Family) , data = sub)
sub$pred2 = exp(-predict(mod2, type="risk"))
qplot(Time, pred2, colour=Food, dat=sub)
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.