Я запускаю glitter logit с использованием пакета lme4. Меня интересуют различные два и три способа взаимодействия и их интерпретации. Чтобы упростить, меня интересуют только коэффициенты фиксированных эффектов.
Мне удалось придумать код для расчета и построения этих эффектов по шкале logit, но у меня возникли проблемы с преобразованием их в прогнозируемый масштаб вероятностей. В конце концов, я хотел бы реплицировать вывод пакета effects.
Пример основан на данных UCLA для больных раком.
library(lme4)
library(ggplot2)
library(plyr)
getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
facmin <- function(n) {
min(as.numeric(levels(n)))
}
facmax <- function(x) {
max(as.numeric(levels(x)))
}
hdp <- read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv")
head(hdp)
hdp <- hdp[complete.cases(hdp),]
hdp <- within(hdp, {
Married <- factor(Married, levels = 0:1, labels = c("no", "yes"))
DID <- factor(DID)
HID <- factor(HID)
CancerStage <- revalue(hdp$CancerStage, c("I"="1", "II"="2", "III"="3", "IV"="4"))
})
До сих пор это все управление данными, функции и пакеты, которые мне нужны.
m <- glmer(remission ~ CancerStage*LengthofStay + Experience +
(1 | DID), data = hdp, family = binomial(link="logit"))
summary(m)
Это модель. Требуется минута, и она сходится со следующим предупреждением:
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0417259 (tol = 0.001, component 1)
Несмотря на то, что я не совсем уверен, должен ли я беспокоиться о предупреждении, я использую оценки для расчета средних предельных эффектов для взаимодействия, представляющего интерес. Сначала я готовлю набор данных для подачи в функцию predict, а затем вычисляю маржинальные эффекты, а также доверительные интервалы, используя параметры фиксированных эффектов.
newdat <- expand.grid(
remission = getmode(hdp$remission),
CancerStage = as.factor(seq(facmin(hdp$CancerStage), facmax(hdp$CancerStage),1)),
LengthofStay = seq(min(hdp$LengthofStay, na.rm=T),max(hdp$LengthofStay, na.rm=T),1),
Experience = mean(hdp$Experience, na.rm=T))
mm <- model.matrix(terms(m), newdat)
newdat$remission <- predict(m, newdat, re.form = NA)
pvar1 <- diag(mm %*% tcrossprod(vcov(m), mm))
cmult <- 1.96
## lower and upper CI
newdat <- data.frame(
newdat, plo = newdat$remission - cmult*sqrt(pvar1),
phi = newdat$remission + cmult*sqrt(pvar1))
Я достаточно уверен, что это правильные оценки по шкале logit, но, возможно, я ошибаюсь. Во всяком случае, это сюжет:
plot_remission <- ggplot(newdat, aes(LengthofStay,
fill=factor(CancerStage), color=factor(CancerStage))) +
geom_ribbon(aes(ymin = plo, ymax = phi), colour=NA, alpha=0.2) +
geom_line(aes(y = remission), size=1.2) +
xlab("Length of Stay") + xlim(c(2, 10)) +
ylab("Probability of Remission") + ylim(c(0.0, 0.5)) +
labs(colour="Cancer Stage", fill="Cancer Stage") +
theme_minimal()
plot_remission
Я думаю, что теперь шкала OY измеряется по шкале logit, но чтобы понять это, я бы хотел преобразовать ее в предсказанные вероятности. На основе wikipedia что-то вроде exp(value)/(exp(value)+1) должно сделать трюк, чтобы добраться до прогнозируемых вероятностей. Хотя я мог бы сделать newdat$remission <- exp(newdat$remission)/(exp(newdat$remission)+1), я не уверен, что как это сделать для доверительных интервалов?
В конце концов я хотел бы получить тот же сюжет, который генерирует пакет effects. То есть:
eff.m <- effect("CancerStage*LengthofStay", m, KR=T)
eff.m <- as.data.frame(eff.m)
plot_remission2 <- ggplot(eff.m, aes(LengthofStay,
fill=factor(CancerStage), color=factor(CancerStage))) +
geom_ribbon(aes(ymin = lower, ymax = upper), colour=NA, alpha=0.2) +
geom_line(aes(y = fit), size=1.2) +
xlab("Length of Stay") + xlim(c(2, 10)) +
ylab("Probability of Remission") + ylim(c(0.0, 0.5)) +
labs(colour="Cancer Stage", fill="Cancer Stage") +
theme_minimal()
plot_remission2
Несмотря на то, что я мог просто использовать пакет effects, он, к сожалению, не компилируется с множеством моделей, которые мне приходилось запускать для моей собственной работы:
Error in model.matrix(mod2) %*% mod2$coefficients :
non-conformable arguments
In addition: Warning message:
In vcov.merMod(mod) :
variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX
Фиксация, которая потребует корректировки процедуры оценки, которую я хотел бы избежать в данный момент. плюс, мне также интересно, что effects на самом деле здесь.
Буду благодарен за любой совет о том, как настроить мой первоначальный синтаксис, чтобы получить предсказанные вероятности!
