Почему я получаю значение 1 для Pr (> | z |) в выводе моего glm? - PullRequest
Новая
       18

Почему я получаю значение 1 для Pr (> | z |) в выводе моего glm?

0 голосов
/ 05 июня 2019

Я пытаюсь приспособить обобщенную линейную модель к моим данным. Я собрал пыльцу из воздуха в течение двух (неравных) временных интервалов: дневной (9 часов) и ночной (15 часов) и определил пыльцу как вишня, другие покрытосеменные и другие голосеменные, исходя из морфологических характеристик.

Когда я смотрю на необработанные данные, наблюдаются некоторые интересные тенденции количества пыльцы вишни в воздухе ночью по сравнению с днем, и эти различия довольно сильно различаются в двух садах, которые я пробовал. Я хочу измерить, как временной интервал влияет на мою модель.

Мои проблемы (неполный список: p)

1) Я не понимаю, почему термин Интервал продолжает давать мне значение 1 для Pr (> | z |). 2) Я не знаю, как объяснить неравные промежутки времени.

Я пробовал этот glm с двумя распределениями ошибок: пуассоновским и биномиальным, в зависимости от формата переменной ответа.

Вот подвыборка моих данных: 

Orchard <- c("CSO", "CSO", "CSO", "HBA", "HBA", "HBA")
Interval <- c("AM", "AM", "AM", "PM", "PM", "PM")
Interval.Duration <- c(9,9,9,15,15,15)
PollenType <- c("Cherry", "Other angiosperm", "Other gymnosperm")
Count <- c(0,2,11,245,124,5,0,2,19,80,38,0,1,0,3,200,150,1)
TotalCount <- c(13,13,13,374,374,374,21,21,21,118,118,118,4,4,4,351,351,351)

df <- data.frame(Orchard, Interval, Interval.Duration, PollenType, Count, TotalCount)
df

# Poisson error distribution model
mod <- glm(Count ~ PollenType + Interval + offset(log(TotalCount)), data = df, family = poisson)
summary(mod) 

Call:
glm(formula = Count ~ PollenType + Interval + offset(log(TotalCount)), 
    family = poisson, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.0076  -3.0242  -0.9525   1.3507   8.8612  

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                -5.158e-01  1.646e-01  -3.134  0.00172 ** 
PollenTypeOther angiosperm -5.096e-01  7.117e-02  -7.159  8.1e-13 ***
PollenTypeOther gymnosperm -2.602e+00  1.660e-01 -15.677  < 2e-16 ***
IntervalPM                 -1.532e-14  1.658e-01   0.000  1.00000    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 742.63  on 17  degrees of freedom
Residual deviance: 240.61  on 14  degrees of freedom
AIC: 313.04

Number of Fisher Scoring iterations: 7


# Binomial error distribution model
mod.bin <- glm(Count/TotalCount ~ PollenType + Interval, data = df, family = binomial)
summary(mod.bin)


Call:
glm(formula = Count/TotalCount ~ PollenType + Interval, family = binomial, 
    data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.04298  -0.88408   0.03032   0.56505   1.02296  

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)
(Intercept)                -5.805e-01  9.912e-01  -0.586    0.558
PollenTypeOther angiosperm -6.754e-01  1.300e+00  -0.519    0.603
PollenTypeOther gymnosperm  2.558e-01  1.187e+00   0.216    0.829
IntervalPM                  7.973e-14  1.016e+00   0.000    1.000

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9.7044  on 17  degrees of freedom
Residual deviance: 9.1325  on 14  degrees of freedom
AIC: 28.299

Number of Fisher Scoring iterations: 4

Вывод 1 для Interval говорит мне, что я сделал что-то не так или что я не задаю свой вопрос должным образом в glm.

Я подозреваю, что проблема может быть связана с форматированием моего набора данных, но я действительно не знаю, с чего начать.

~~ РЕДАКТИРОВАТЬ ~~

Спасибо @ Eric-Scott за подробный отзыв. Я переработал фиктивные данные, чтобы они отражали только количество вишен, и сделал подмножества для каждого сада, потому что не думаю, что имеет смысл включать данные из обоих садов в одну модель. Я использовал «Interval.Duration» в качестве смещения, как вы рекомендовали. 

Orchard <- c("HBA", "CSO")
Interval <- c("AM", "AM", "PM","PM")
Interval.Duration <- c(9,15)
Count <- c(13,4,245,0,80,2,98,1,200,1,530,1,196,2,311,1)
TotalCount <- c(38,7,374,21,118,15,144,4,351,10,884,12,338,34,490,15)

df <- data.frame(Orchard, Interval, Interval.Duration, Count, TotalCount)
df
df.H <- subset(df, Orchard == "HBA")
df.C <- subset(df, Orchard == "CSO")

#HBA data
modH <- glm(Count ~ Interval + offset(log(Interval.Duration)),  data = df.H, family = quasipoisson)

Call:
glm(formula = Count ~ Interval + offset(log(Interval.Duration)), 
    family = quasipoisson, data = df.H)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-13.392   -6.225   -1.096    6.204   12.226  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.6088     0.4263   6.120 0.000869 ***
IntervalPM    0.8843     0.5067   1.745 0.131574    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 88.85893)

    Null deviance: 892.38  on 7  degrees of freedom
Residual deviance: 594.73  on 6  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

#CSO data
modC <- glm(Count ~ Interval + offset(log(Interval.Duration)), data = df.C, family = quasipoisson)

Call:
glm(formula = Count ~ Interval + offset(log(Interval.Duration)), 
    family = quasipoisson, data = df.C)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.22474  -0.36170   0.05231   0.27453   1.05020  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8971     0.2400  -7.904 0.000218 ***
IntervalPM   -1.0986     0.4801  -2.289 0.062070 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 0.5185185)

    Null deviance: 6.9044  on 7  degrees of freedom
Residual deviance: 3.7649  on 6  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Используя фиктивные данные, я вижу избыточную дисперсию и недостаточную дисперсию в моделях «HBA» и «CSO» соответственно (хотя, когда я запускаю модель на полном наборе данных, в обоих случаях наблюдается избыточная дисперсия).

Во-первых, я подгоняю модель с отрицательным биномиальным распределением, так как знаю, что модель с нулевым раздувом не поможет с избыточной дисперсией в наборе данных HBA (нулей нет). 

library(MASS)
mod.NB <- glm.nb(formula = Count ~ Interval + offset(log(Interval.Duration)), data = df)
summary(mod.NB)

Call:
glm.nb(formula = Count ~ Interval + offset(log(Interval.Duration)), 
    data = df, init.theta = 0.2822584211, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9557  -1.4501  -0.8920   0.3307   0.8569  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   1.9258     0.6667   2.889  0.00387 **
IntervalPM    0.8753     0.9423   0.929  0.35294   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2823) family taken to be 1)

    Null deviance: 21.300  on 15  degrees of freedom
Residual deviance: 20.463  on 14  degrees of freedom
AIC: 166.75

Number of Fisher Scoring iterations: 1


              Theta:  0.2823 
          Std. Err.:  0.0836 

 2 x log-likelihood:  -160.7510

Следуя советам возводить в степень коэффициенты из https://rpubs.com/kaz_yos/pscl-2 Я "возводил в степень" (это слово?) Коэффициенты из моей модели: 

exp(coef(mod.NB))

(Intercept)  IntervalPM 
   6.860548    2.399691

Если я правильно интерпретирую этот вывод, среднее число пыльцевых зерен вишни составляет 6,86, а ночью пыльца вишневого дерева в 2,4 раза больше. Статистически, однако, увеличение пыльцы в ночное время не является статистически большим, чем в течение дня. (NB: потерпите меня, большая часть моего замешательства связана с интерпретацией результатов).

Я запустил модель с нулевым раздуванием, чтобы справиться с беспорядком нулей в наборе данных "CSO", но мне также не ясно, как интерпретировать результаты здесь. 

library(pscl)
#HBA data set
modzH <- zeroinfl(formula = Count ~ Interval, data = df.H, dist = "negbin")

В наборе данных HBA нет нулей, поэтому появляется следующая ошибка (отсюда и вышеупомянутый NegBin glm). 

Error in zeroinfl(formula = Count ~ Interval, data = df.H, dist = "negbin") :

недопустимая зависимая переменная, минимальное число не равно нулю 

#CSO data
modzC <- zeroinfl(formula = Count ~ Interval, data = df.C, dist = "negbin")
summary(modzC)

Call:
zeroinfl(formula = Count ~ Interval + offset(log(Interval.Duration)), data = df.C, 
    dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.8660 -0.3333  0.0610  0.2887  1.1666 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.8971     0.3333  -5.691 1.26e-08 ***
IntervalPM   -1.0986     0.6667  -1.648   0.0994 .  
Log(theta)   13.6976   510.5736   0.027   0.9786    

Zero-inflation model coefficients (binomial with logit link):
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    -28.49  198393.43       0        1
IntervalPM      12.25  198394.44       0        1
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 888750.384 
Number of iterations in BFGS optimization: 23 
Log-likelihood: -10.13 on 5 Df

По совету Ferran Paüls Vergés Я снова возводил в степень коэффициенты из модели. (Это потому, что выходные данные модели дают коэффициенты в лог-коэффициентах?) 

## Exponentiated coefficients
expCoef <- exp(coef((modzC)))
expCoef

# count_(Intercept)  count_IntervalPM  zero_(Intercept)   zero_IntervalPM 
# 1.500015e-01      3.333271e-01      4.234422e-13      2.092359e+05

И (если я понимаю, как их применять) базовое число пыльцы вишни в образце, в котором содержится пыльца, составляет 1,500015e-01, а шансы базовой линии образца, имеющего ноль зерен пыльцы вишни, составляют 4,234422e-13? Я действительно не уверен в этом предложении, оно не имеет смысла.

1 Ответ

0 голосов
/ 05 июня 2019

Ваша переменная Count определенно перераспределена, что означает, что пуассоновский GLM может давать неверные результаты.Я могу видеть избыточную дисперсию несколькими способами: во-первых, остаточное отклонение намного больше, чем степени свободы, во-вторых, если вы запускаете модель как квазипуассон (https://www.theanalysisfactor.com/glm-r-overdispersion-count-regression/) параметр дисперсии равен 46, что являетсянамного больше 1, что предполагается для распределения Пуассона, и построение гистограммы дает распределение с длинным хвостом. У вас много нулей, и вы можете извлечь выгоду из использования модели с нулевым завышением. Модели с нулевым раздувом предполагают, что вы обнуляетеОн получен из двух источников: нули обнаружения и «настоящие» нули, то есть, возможно, в этом первом наблюдении действительно не было никакой пыльцы вишни, или, возможно, она была там, и вы просто ее не обнаружили.

Я думаю, что вы на правильном пути с первой моделью. Смещения обычно используются для учета различий в усилиях по сэмплированию. Для меня наиболее естественным кандидатом на смещение является interval.length, но я могу понять, почему вы могли быиспользовать TotalCount.

Глядя на коробочный сюжет, я думаю, что кажется разумным, что нет мВлияние интервала (р = 1), но если вы включите взаимодействие между интервалом и типом пыльцы, взаимодействие является значительным.Если вы только интересуетесь пыльцой Cherry, и вы просто используете другие типы пыльцы, чтобы получить TotalPollen для смещения, тогда отфильтруйте df так, чтобы это было только cherry и drop PollenType из модели, и яСпорим, вы увидите эффект интервала. 

mod <- glm(Count ~ PollenType + Interval, offset = log(TotalCount), data = df, family = poisson)
car::Anova(mod) 

mod2 <- glm(Count ~ PollenType*Interval, offset = log(TotalCount), data = df, family = quasipoisson)
car::Anova(mod2)

Но первое, что я бы сделал, это попытайтесь справиться с избыточной дисперсией.

Ответ на редактирование

Во-первых, вы сказали, что будете делать отдельные модели для каждого сада, но не сделали этого (data = df).Если вы используете data = df, тогда вы можете включить фруктовый сад в качестве фактора, и вы получите значительный интервал: взаимодействие фруктового сада, что может быть вам интересно?(попробуйте посмотреть ggplot(df, aes(x = Interval, y = log(Count), color = Orchard)) + geom_boxplot()).

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

Я бы постарался упростить ситуацию, применив одну и ту же модель к обоимфруктовые сады.Тот факт, что в одном из ваших садов нет нулей, заставляет меня думать, что в другом саду на самом деле нет двух источников нулей.То есть нули, которые у вас есть, являются «реальными» и указывают на отсутствие пыльцы, а не просто на неспособность обнаружить пыльцу.Так может не ноль раздутый в конце концов?Вы должны будете судить сами, если у вас больше нулей, чем вы ожидаете при отрицательном биномиальном распределении.Кроме того, как вы заметили, интерпретация моделей с нулевой раздувкой общеизвестно трудна.

Что касается интерпретации выходных данных отрицательного бинома glm, я бы лично избегал слишком много смотреть на summary().Вы можете использовать его для получения коэффициентов (или вы можете использовать coef()), но для проверки значимости интервала используйте либо car::Anova(mod.NB), либо создайте пустую модель (mod.NB.0 <- glm.nb(formula = Count ~ 1 + offset(log(Interval.Duration)), data = df)) и используйте lmtest::lrtest(mod.NB, mod.NB.0).

...