Меня интересует получение установленных значений в заданных местах из модели clogit. Это включает ответ уровня населения и доверительные интервалы вокруг него. Например, у меня есть данные, что выглядит примерно так:Как получить установленные значения от модели clogit
set.seed(1)
data <- data.frame(Used = rep(c(1,0,0,0),1250),
Open = round(runif(5000,0,50),0),
Activity = rep(sample(runif(24,.5,1.75),1250, replace=T), each=4),
Strata = rep(1:1250,each=4))
В рамках модели Clogit, активность не меняется в пределах слоев, таким образом, нет никакой активности основного эффекта.
mod <- clogit(Used ~ Open + I(Open*Activity) + strata(Strata),data=data)
То, что я хочу сделать, это построить NewData кадр, в котором я могу в конечном итоге сюжет маргинальные подогнанные значения в определенных местах открытого похожих на NewData дизайн в традиционной GLM модели: например,
newdata <- data.frame(Open = seq(0,50,1),
Activity = rep(max(data$Activity),51))
Однако, когда я пытаюсь запустить функцию предсказать на clogit, я получаю следующее сообщение об ошибке:
fit<-predict(mod,newdata=newdata,type = "expected")
ошибка в Surv (Rep (1, 5000L), б): объект 'б' не найден
Я понимаю, это связано с тем, что clogit в r выполняется через Cox.ph, и, таким образом, функция прогнозирования пытается предсказать относительные риски между парами субъектов внутри одной и той же страты (в данном случае = Используется).
Мой вопрос, однако, если есть способ обойти это. Это легко сделать в Stata (с использованием команды полей) и вручную в Excel, однако я бы хотел автоматизировать в R, поскольку там запрограммировано все остальное. Я также создал это вручную в R (пример кода ниже), однако я продолжаю заканчивать то, что кажется неправильным CI в моих реальных данных, поэтому я хотел бы полагаться на функцию прогнозирования, если это возможно. Мой код для ручного предсказания:
coef<-data.frame(coef = summary(mod)$coefficients[,1],
se= summary(mod)$coefficients[,3])
coef$se <-summary(mod)$coefficients[,4]
coef$UpCI <- coef[,1] + (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity
coef$LowCI <-coef[,1] - (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity
fitted<-data.frame(Open= seq(0,50,2),
Activity=rep(max(data$Activity),26))
fitted$Marginal <- exp(coef[1,1]*fitted$Open +
coef[2,1]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,1]*fitted$Open +
coef[2,1]*fitted$Open*fitted$Activity))
fitted$UpCI <- exp(coef[1,3]*fitted$Open +
coef[2,3]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,3]*fitted$Open +
coef[2,3]*fitted$Open*fitted$Activity))
fitted$LowCI <- exp(coef[1,4]*fitted$Open +
coef[2,4]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,4]*fitted$Open +
coef[2,4]*fitted$Open*fitted$Activity))
Мой конечный продукт будет идеально выглядеть примерно так, но продукт функции прогнозирования ....
Example output of fitted values.
Это то, что автор, Т. Ламли, должен был сказать в 2006 году в rhelp, когда его попросили использовать функцию 'preview.clogit':« Я не думаю, что это будет возможно. Пункт условного логистический регрессии, что вероятности зависят от параметров пластовых , которые не могут быть оценены точно. условная вероятность того, удаляет этих параметры, но полученная модель не содержит достаточно информации для оценки вероятности. \t -Томаса» –