ggplot GLM встроенная кривая без взаимодействия

Question

Я хочу добавить установленную функцию из GLM на ggplot. По умолчанию он автоматически создает график со взаимодействием. Я задаюсь вопросом, могу ли я построить функцию, созданную из модели, без взаимодействия. Например,

dta <- read.csv("http://www.ats.ucla.edu/stat/data/poisson_sim.csv") 
dta <- within(dta, { 
    prog <- factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational")) 
    id <- factor(id) 
}) 

plt <- ggplot(dta, aes(math, num_awards, col = prog)) + 
    geom_point(size = 2) + 
    geom_smooth(method = "glm", , se = F, 
     method.args = list(family = "poisson")) 

print(plt) 

дает сюжет с взаимодействием,

Однако я хочу участок от модели,

`num_awards` = ß0 + ß1*`math` + ß2*`prog` + error 

Я пытался получить это таким образом,

mod <- glm(num_awards ~ math + prog, data = dta, family = "poisson") 

fun.gen <- function(awd) exp(mod$coef[1] + mod$coef[2] * awd) 
fun.acd <- function(awd) exp(mod$coef[1] + mod$coef[2] * awd + mod$coef[3]) 
fun.voc <- function(awd) exp(mod$coef[1] + mod$coef[2] * awd + mod$coef[4]) 

ggplot(dta, aes(math, num_awards, col = prog)) + 
    geom_point() + 
    stat_function(fun = fun.gen, col = "red") + 
    stat_function(fun = fun.acd, col = "green") + 
    stat_function(fun = fun.voc, col = "blue") + 
    geom_smooth(method = "glm", se = F, 
     method.args = list(family = "poisson"), linetype = "dashed") 

Выходной график

Есть ли простой способ в ggplot, чтобы сделать это эффективно?

Answer 1

В этом нет ничего плохого, чтобы обмануть geom_smooth(), но вы можете сделать немного лучше, чем вы сделали. Вы все еще должны подобрать модель себя и добавить строки, но вы можете использовать метод predict() для генерации предсказания и загружать их в кадр данных, с той же структурой, что и исходные данные ...

mod <- glm(num_awards ~ math + prog, data = dta, family = "poisson") 
## generate prediction frame 
pframe <- with(dta, 
      expand.grid(math=seq(min(math),max(math),length=51), 
         prog=levels(prog))) 
## add predicted values (on response scale) to prediction frame 
pframe$num_awards <- predict(mod,newdata=pframe,type="response") 

ggplot(dta, aes(math, num_awards, col = prog)) + 
    geom_point() + 
    geom_smooth(method = "glm", se = FALSE, 
     method.args = list(family = "poisson"), linetype = "dashed")+ 
    geom_line(data=pframe) ## use prediction data here 
      ## (inherits aesthetics etc. from main ggplot call) 

(единственная разница здесь в том, что так, как я это сделал, прогнозы охватывают весь горизонтальный диапазон для всех групп, как если бы вы указали fullrange=TRUE в geom_smooth()).

В принципе кажется, что sjPlot package должен иметь возможность обрабатывать подобные вещи, но похоже, что соответствующий бит кода для этого типа графика жестко запрограммирован, чтобы принять биномиальный GLM ... oh well ,

Answer 2

Я не уверен, но вы написали «без взаимодействия» - может быть, вы ищете участки эффектов? (Если нет, извините, что я предполагаю что-то совершенно неправильное ...)

Для этого, например, можно использовать effects package.

dta <- read.csv("http://www.ats.ucla.edu/stat/data/poisson_sim.csv") 
dta <- within(dta, { 
    prog <- factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational")) 
    id <- factor(id) 
}) 

mod <- glm(num_awards ~ math + prog, data = dta, family = "poisson") 

library(effects) 
plot(allEffects(mod)) 

Другим вариантом был бы пакет sjPlot, так как Бен предложил - однако, текущая версия на CRAN поддерживает только модель логистической регрессии должным образом для эффекта участков. Но в текущей версии разработки на GitHub я добавил поддержку различных семейств моделей и функций ссылок, поэтому, если вам нравится, вы можете загрузить этот снимок.Пакет sjPlot использует ggplot вместо решетки (который используется пакетом эффектов, я думаю):

sjp.glm(mod, type = "eff", show.ci = T) 

Или не-граненый образом:

sjp.glm(mod, type = "eff", facet.grid = F, show.ci = T) 

Answer 3

Идея Бена о прогнозируемом значении ответа для конкретных условий модели вдохновила меня на улучшение функции type = "y.pc" функции sjp.glm. Новое обновление - on GitHub, с номером версии 1.9.4-3.

Теперь вы можете построить предсказанные значения для конкретных условий, один, который используется вдоль оси х, а второй один используется в качестве группировки фактора:

sjp.glm(mod, type = "y.pc", vars = c("math", "prog")) 

, который дает вам следующий сюжет:

Аргумент vars необходим в случае, если ваша модель имеет более двух терминов, чтобы указать термин для оси X-оси и термин для группировки.

Вы также можете фасет группы:

sjp.glm(mod, type = "y.pc", vars = c("math", "prog"), show.ci = T, facet.grid = T)