2

I am having trouble interpreting the results of my model. I have run a mixed effects model with two fixed effects, two random effects using frailty(), and two interaction terms (fixed:fixed and fixed:random). I am interested in reconstructing a function to predict survival time for individuals, so the estimates of variance for the random effects terms are important (Hence why I cannot use cluster() to identify them). Here is some sample code to generate what my dataset looks like on a very reduced scale:

library(survival)
library(ggplot2)
library(ggfortify)

# Generating data

Hospital <- rep(c("Facility 1", "Facility 2", "Facility 3", "Facility 4"), each = 32)
Doctor <- rep(c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32"), each = 4)
Treatment <- rep(c("A", "B"), each = 16, 4)
City <- rep(c("LA", "NY"), each = 64)
Time1 <- c(NA,1,2,3,4,1,2,3,1,4,NA,1,3,4,1,2,4,3,4,3,3,4,2,4,4,2,3,3,4,4,4,3,3,2,1,3,1,1,2,1,NA,1,NA,1,1,NA,2,3,3,4,2,4,2,3,3,4,3,2,2,2,1,3,4,4,1,2,NA,2,1,1,2,1,2,1,1,1,1,2,NA,3,4,3,2,3,4,4,2,2,4,4,4,3,2,3,3,4,1,1,2,NA,NA,1,1,1,2,1,2,3,1,1,2,3,4,3,4,3,3,4,3,4,1,3,4,3,3,4,4,4)
Time2 <- Time1 + 1
data = data.frame(Hospital, Doctor, Treatment, City, Time1, Time2)
data$Time2[is.na(data$Time2)] <- 1
data$Time2[data$Time2 == 5] <- NA

Using the ggplot libraries in the next section, it is obvious that there is an effect on survival due to Treatment, where individuals with Treatment B live longer (confirmed with p < 0.01) and that there is no significant effect due to City (p >> 0.05) or the interaction between Treatment and City (p >> 0.05). My issue lies in understanding the output for the gamma terms. Is this telling me the effect size of the random effects? If so, why do I not see outputs for Doctor in addition to Hospital? Why do some of the interaction terms output NA's and 0's?

# Model fit and results

sobj <- Surv(data$Time1, data$Time2, type = 'interval2')
data$label <- paste(data$City, data$Treatment)
autoplot(survfit(sobj ~ data$label))
fit <- survreg(sobj ~ data$Treatment + data$City + data$Treatment:data$City + 
frailty(data$Hospital) + frailty(data$Doctor) + frailty(data$Hospital):data$City, dist = "weibull")
summary(fit)

I'm not sure if my problem lies in the way I am coding my model or in my understanding of statistics (probably both), but seeing errors in the output makes me hesitate to trust even the conclusions about the fixed effects. Please let me know if there is a better place to post my questions.

M-Student
  • 21
  • 1
  • I suspect you have not looked at the data with a tabulation. The values of -7.04e-16 and -2.30e-16 are each effectively 0. And then you have the NA values that indicate removal due to aliasing. I think the data does not support such a complex model. It's also not clear what sort of "events" this data is supposed to be predicting or summarizing. There's no event vector in the Surv argument. – IRTFM May 25 '17 at 04:38
  • I also get: "Error: Objects of type survfit not supported by autoplot." – IRTFM May 25 '17 at 04:45
  • Voting to migrate to CrossValidated.com. Another venue might be the Mixed Models R SIG list. – IRTFM May 25 '17 at 15:58
  • @42- When using 'interval2' the event is determined based on the entries for T1 and T2 (for left censored: [NA,t2], for interval censored: [t1,t2], etc.). I get the same exact output when I specify 'interval' instead and include the event type vector. It's likely that the model is too complex. I have not been able to replicate that error. Thank you for the suggestion of another venue. – M-Student May 29 '17 at 07:48
  • It's both less clear and less safe to supply "naked vectors" to the formula in `survreg`. Therneau recommends code like this (which give the same results in this case: `fit <- survreg(Surv(Time1, Time2, type = 'interval2') ~ Treatment + City + Treatment:City + frailty(Hospital) + frailty(Doctor) + frailty(Hospital):City, data=data, dist = "weibull")` and the formula interface allows you to collapse those interaction terms to this more compact expression: `fit <- survreg(Surv(Time1, Time2, type = 'interval2') ~ Treatment* City + frailty(Doctor*Hospital), data=data, dist = "weibull")` – IRTFM May 29 '17 at 20:31

0 Answers0