używam modele glmer logitowe użyciu pakietu lme4. Interesują mnie różne dwu- i trójdrożne efekty interakcji i ich interpretacje. Aby uprościć, interesują mnie tylko stałe współczynniki efektów.glmer logarytmicznej - efekty interakcji w skali prawdopodobieństwa (replikowania `` ze skutkami w predict`)
udało mi się wymyślić kodu obliczyć i wykreślić te skutki w skali logit, ale mam problemy przekształcając je do przewidywanej skali prawdopodobieństwa. Ostatecznie chciałbym powtórzyć wyjście pakietu effects.
Przykład opiera się na UCLA's data on cancer patients.
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"))
})
Do tego czasu wszystko to zarządzanie danymi, funkcje i pakiety, których potrzebuję.
m <- glmer(remission ~ CancerStage*LengthofStay + Experience +
(1 | DID), data = hdp, family = binomial(link="logit"))
summary(m)
To jest model. To trwa chwilę i jest zbieżny z następującym ostrzeżeniem:
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)
Chociaż nie jestem pewien, czy powinienem się martwić o ostrzeżeniu, używam szacunki wykreślić średnie marginalne skutki dla interakcji interesów. Najpierw przygotowuje się zestaw danych do paszy w funkcji predict, a następnie obliczyć efektów brzegowych, jak również przedziały ufności przy użyciu stałych parametrów efektów.
Jestem przekonany, że to poprawne szacunki na skali logitowej, ale może się mylę. Tak czy inaczej, jest to działka:
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
myślę teraz skala OY jest mierzona na skali logitowej ale sens to chciałbym, aby przekształcić go do przewidywanych prawdopodobieństw. Opierając się na wikipedia, coś takiego jak exp(value)/(exp(value)+1) powinno wystarczyć, aby uzyskać przewidywane prawdopodobieństwa. Chociaż mogłem zrobić newdat$remission <- exp(newdat$remission)/(exp(newdat$remission)+1) nie jestem pewien jak powinienem to zrobić dla przedziałów ufności ?.
Ostatecznie chciałbym dostać się do tej samej działki, co generuje pakiet effects. Czyli:
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
Choć może po prostu użyć pakietu effects, to niestety nie skompilować z wielu modeli musiałem biec do własnej pracy:
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
mocowania, które wymagają dostosowania procedury szacowania, które w tej chwili chciałbym uniknąć. plus, jestem także ciekawy, co naprawdę robi tutaj. Byłbym wdzięczny za porady, jak poprawić moją początkową składnię, aby uzyskać przewidywane prawdopodobieństwa!
Myślę, że twoja fabuła będzie łatwiejsza do odczytania, jeśli zrobisz coś takiego: 'ggplot (n ewdat, aes (LengthofStay, fill = factor (CancerStage), color = factor (CancerStage))) + geom_ribbon (aes (ymin = plo, ymax = phi), color = NA, alpha = 0.2) + geom_line (aes (y = remisja), rozmiar = 1,2) + xlab ("długość pobytu") + ylab ("Prawdopodobieństwo remisji") + laboratoria (kolor = "etap raka", wypełnienie = "etap raka") + temat_minimalny () ' – eipi10
Powinieneś zdecydowanie martwić się ostrzeżeniem o zbieżności. –
Naprawdę nie rozumiem, dlaczego tak trudno jest odpowiedzieć na pytanie ... Czy jest coś niejasnego w tym, o co proszę? – eborbath