рисунок для регрессии
Здравствуйте! У меня случилась маленькая проблемка с визуализацией к регрессии Пуассона.
Модель:
M3 <- glm(Macross ~ location+Sex, data = a, family = "quasipoisson")
Anova(M3)
Analysis of Deviance Table (Type II tests)
Response: Macross
LR Chisq Df Pr(>Chisq)
location 28.165 4 1.155e-05 ***
Sex 15.291 1 9.216e-05 ***
Мне удаётся построить только рисунок одновременно для факторов "Локализация" и "Пол" с помощью кода:
NEWDATA <- with(a, expand.grid(location=unique(location), Sex=unique(Sex)))
pred <- predict(M3,newdata= NEWDATA,se.fit=TRUE, type = "link")
make_ci <- function(pred, data){
fit <- pred$fit
lower <- fit - 1.96*pred$se.fit
upper <- fit + 1.96*pred$se.fit
return(data.frame(exp(fit), exp(lower), exp(upper), data))
}
my_pred <- make_ci(pred, NEWDATA)
limits <- aes(x = factor(location), ymax = my_pred$exp.upper., ymin = my_pred$exp.lower., group = Sex)
ggplot(my_pred, aes(x = factor(location), y = exp.fit., color = Sex))+
geom_errorbar(limits, position = position_dodge(width = 0.75), color = "black")+
geom_point(size = 8, position = position_dodge(width = 0.75), shape = 16)+
geom_point(data = a, aes(x = factor(location), y = Macross, colour = Sex), size = 8, shape = 17, alpha = 0.7,position = position_dodge(width = 0.75))
Вот "грубая" версия рисунка:
значения для доверительных интервалов:
> my_pred
exp.fit. exp.lower. exp.upper. location Sex
1 30.961669 7.6722541 124.946975 Тункинский ПП f
2 78.220600 52.5776978 116.369915 Долганская яма f
3 19.626074 7.8067394 49.339778 Охотничья f
4 8.158858 0.6930274 96.052420 Тигирек f
5 11.424414 4.3469659 30.024904 Барсуковская пещера f
6 9.012777 2.2112737 36.734553 Тункинский ПП m
7 22.769600 12.9430298 40.056672 Долганская яма m
8 5.713046 2.4043790 13.574770 Охотничья m
9 2.375000 0.2196026 25.685599 Тигирек m
10 3.325586 1.1329687 9.761543 Барсуковская пещера m
Но мне нужно построить ОТДЕЛЬНО рисунки для главных эффектов: один для локализации, другой - для пола. Но с этим кодом не очень-то и выходит.
Конечно, я могу построить по модели для каждого главного эффекта и получить рисунки, например:
M4 <- glm(Macross ~ location, data = a, family = "quasipoisson")
но, как мне кажется, доверительные интервалы этой модели не будут соответствовать модели M3.
Возможно есть пути, как обойти эту проблему и построить рисунки для каждого главного эффекта по модели M3? Взаимодействие между факторами статистически не значимо.
Заранее спасибо за помощь!!!
https://dropmefiles.com/C0wL5 - ссылка на файл с данными
Модель:
M3 <- glm(Macross ~ location+Sex, data = a, family = "quasipoisson")
Anova(M3)
Analysis of Deviance Table (Type II tests)
Response: Macross
LR Chisq Df Pr(>Chisq)
location 28.165 4 1.155e-05 ***
Sex 15.291 1 9.216e-05 ***
Мне удаётся построить только рисунок одновременно для факторов "Локализация" и "Пол" с помощью кода:
NEWDATA <- with(a, expand.grid(location=unique(location), Sex=unique(Sex)))
pred <- predict(M3,newdata= NEWDATA,se.fit=TRUE, type = "link")
make_ci <- function(pred, data){
fit <- pred$fit
lower <- fit - 1.96*pred$se.fit
upper <- fit + 1.96*pred$se.fit
return(data.frame(exp(fit), exp(lower), exp(upper), data))
}
my_pred <- make_ci(pred, NEWDATA)
limits <- aes(x = factor(location), ymax = my_pred$exp.upper., ymin = my_pred$exp.lower., group = Sex)
ggplot(my_pred, aes(x = factor(location), y = exp.fit., color = Sex))+
geom_errorbar(limits, position = position_dodge(width = 0.75), color = "black")+
geom_point(size = 8, position = position_dodge(width = 0.75), shape = 16)+
geom_point(data = a, aes(x = factor(location), y = Macross, colour = Sex), size = 8, shape = 17, alpha = 0.7,position = position_dodge(width = 0.75))
Вот "грубая" версия рисунка:
значения для доверительных интервалов:
> my_pred
exp.fit. exp.lower. exp.upper. location Sex
1 30.961669 7.6722541 124.946975 Тункинский ПП f
2 78.220600 52.5776978 116.369915 Долганская яма f
3 19.626074 7.8067394 49.339778 Охотничья f
4 8.158858 0.6930274 96.052420 Тигирек f
5 11.424414 4.3469659 30.024904 Барсуковская пещера f
6 9.012777 2.2112737 36.734553 Тункинский ПП m
7 22.769600 12.9430298 40.056672 Долганская яма m
8 5.713046 2.4043790 13.574770 Охотничья m
9 2.375000 0.2196026 25.685599 Тигирек m
10 3.325586 1.1329687 9.761543 Барсуковская пещера m
Но мне нужно построить ОТДЕЛЬНО рисунки для главных эффектов: один для локализации, другой - для пола. Но с этим кодом не очень-то и выходит.
Конечно, я могу построить по модели для каждого главного эффекта и получить рисунки, например:
M4 <- glm(Macross ~ location, data = a, family = "quasipoisson")
но, как мне кажется, доверительные интервалы этой модели не будут соответствовать модели M3.
Возможно есть пути, как обойти эту проблему и построить рисунки для каждого главного эффекта по модели M3? Взаимодействие между факторами статистически не значимо.
Заранее спасибо за помощь!!!
https://dropmefiles.com/C0wL5 - ссылка на файл с данными