我无法解释我的模型的结果。我已经运行了一个混合效应模型,其中包含两个固定效应、两个使用 frailty() 的随机效应和两个交互项(固定:固定和固定:随机)。我有兴趣重建一个函数来预测个人的生存时间,因此随机效应项的方差估计很重要(因此我不能使用 cluster() 来识别它们)。以下是一些示例代码,可以在非常缩小的范围内生成我的数据集的样子:
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
使用下一节中的 ggplot 库,很明显,治疗对生存有影响,治疗 B 的个体寿命更长(证实为 p < 0.01),而城市没有显着影响(p > > 0.05)或治疗与城市之间的相互作用(p >> 0.05)。我的问题在于理解伽马项的输出。这是在告诉我随机效应的效应大小吗?如果是这样,为什么除了医院之外我看不到医生的输出?为什么有些交互项会输出 NA 和 0?
# 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)
我不确定我的问题是出在我编写模型的方式上还是出在我对统计数据的理解上(可能两者兼而有之),但是看到输出中的错误让我对即使是关于固定效应的结论也犹豫不决。如果有更好的地方可以发布我的问题,请告诉我。