このポストでは、Salamandersデータに含まれるサンショウウオの観測カウントデータを題材に、モデルを段階的に拡張しながら当てはまりを比較する。具体的には、まず通常のポアソン混合モデルを当て、次にゼロ過剰を切片のみで導入したモデル、最後にゼロ過剰確率を説明変数で動かすモデルへ進む。AICと残差診断を用いて、ゼロの多さをどこまで説明できるかを検討する。
Salamandersデータ
Salamanders は、glmmTMB に同梱されている「渓流でのサンショウウオの反復カウント」データである。複数の地点（site）を繰り返し調査し、地点の環境要因（site covariates）と、調査時点の条件（sampling covariates）を併せ持つ。23地点をそれぞれ複数回サンプリングした構造を持つと説明されている。
N（観測数）は 644（10変数）である。
変数名と内容（データフレームの列）
site：反復サンプルを取った地点名（ロケーションの識別子）mined：山頂除去型の石炭採掘（mountaintop removal coal mining）の影響を受けた地点かどうかを表す因子cover：渓流内のカバー（隠れ場所になる物体）の量（スケール済み）sample：反復サンプル（何回目の採取か）DOP：降水からの日数（Days since precipitation；スケール済み）Wtemp：水温（スケール済み）DOY：年内日（day of year；スケール済み）spp：種名（省略形）で、場合によっては生活史段階（life stage）も含む可能性があるcount：観測されたサンショウウオの個体数（カウント）
※ヘルプ上は「10変数」となっているが、上の列挙は9項目しか見えない（Formatの記載が10と一致していない）ため、実際にRで names(Salamanders) を見て「もう1列」が何かを確認するのが確実である。
コード
library(glmmTMB) library(DHARMa) # 診断用 library(dplyr) # 任意
データ読み込み（glmmTMBに同梱）
data("Salamanders", package = "glmmTMB")
dat <- Salamanders
## 列名と型を確認
str(dat)
summary(dat)
'data.frame': 644 obs. of 9 variables:
$ site : Ord.factor w/ 23 levels "R-1"<"R-2"<"R-3"<..: 13 14 15 1 2 3 4 5 6 7 ...
$ mined : Factor w/ 2 levels "yes","no": 1 1 1 2 2 2 2 2 2 2 ...
$ cover : num -1.442 0.298 0.398 -0.448 0.597 ...
$ sample: int 1 1 1 1 1 1 1 1 1 1 ...
$ DOP : num -0.596 -0.596 -1.191 0 0.596 ...
$ Wtemp : num -1.2294 0.0848 1.0142 -3.0234 -0.1443 ...
$ DOY : num -1.497 -1.497 -1.294 -2.712 -0.687 ...
$ spp : Factor w/ 7 levels "GP","PR","DM",..: 1 1 1 1 1 1 1 1 1 1 ...
$ count : int 0 0 0 2 2 1 1 2 4 1 ...
site mined cover sample DOP Wtemp DOY spp count
R-1 : 28 yes:308 Min. :-1.59152 Min. :1.00 Min. :-2.1984 Min. :-3.0234 Min. :-2.7122 GP :92 Min. : 0.000
R-2 : 28 no :336 1st Qu.:-0.69629 1st Qu.:1.75 1st Qu.:-0.3018 1st Qu.:-0.6139 1st Qu.:-0.5653 PR :92 1st Qu.: 0.000
R-3 : 28 Median :-0.04974 Median :2.50 Median :-0.0916 Median : 0.0370 Median :-0.0590 DM :92 Median : 0.000
R-4 : 28 Mean : 0.00000 Mean :2.50 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 EC-A :92 Mean : 1.323
R-5 : 28 3rd Qu.: 0.59682 3rd Qu.:3.25 3rd Qu.: 0.0000 3rd Qu.: 0.6032 3rd Qu.: 0.9739 EC-L :92 3rd Qu.: 2.000
R-6 : 28 Max. : 1.88993 Max. :4.00 Max. : 3.1691 Max. : 2.2094 Max. : 1.4600 DES-L:92 Max. :36.000
(Other):476 DF :92
まず0の多さを確認
y <- dat$count
cat("N =", length(y), "\n")
cat("zeros =", sum(y == 0), "\n")
cat("zero proportion =", mean(y == 0), "\n")
cat("mean =", mean(y), " var =", var(y), "\n")
N = 644 zeros = 387 zero proportion = 0.6009317 mean = 1.322981 var = 6.946843
観測644件のうち0が387件で、0の割合が約0.601とかなり高いことを示している。ポアソン回帰の素朴な前提（平均＝分散）に照らすと、平均が約1.323なのに分散が約6.947であり、分散が平均の約5.25倍も大きい。したがってこのデータは「0が多い」だけでなく、一般にいう過分散（overdispersion）も強い。結論として、通常のポアソン単体で押し切るのはかなり無理があり、少なくとも負の二項（NB）か、ゼロ過剰（ZIP/ZINB）を検討すべき状況だと言える。
種別ごとの0比率（sppごとにゼロが偏ることが多い）
if ("spp" %in% names(dat)) {
zp_by_spp <- tapply(dat$count == 0, dat$spp, mean)
print(sort(zp_by_spp, decreasing = TRUE))
}
PR EC-A GP DM EC-L DF DES-L 0.8478261 0.7717391 0.5869565 0.5108696 0.5108696 0.5108696 0.4673913
0が「全体に均一に多い」のではなく、種（あるいは生活史段階を含むカテゴリ）によって0の出やすさがかなり違うことを示している。具体的には PR は0が約0.848と極端に多く、EC-A も0が約0.772と高い一方、DES-L は0が約0.467で相対的に低い。これは、ゼロの発生が単なる偶然ではなく、種や段階に依存する「検出されにくさ／不在が多い／生息条件が厳しい」等の構造的要因を含む可能性を示唆する。モデル化の観点では、spp を単にカウント部（平均λ）の説明変数として入れるだけでなく、ゼロ過剰部（π）にも spp を入れる価値があるタイプの出方である。
モデルの比較
## (A) 普通のポアソン：count ~ mined + cover + DOP + Wtemp + DOY + spp + (1|site) m_pois <- glmmTMB( count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site), family = poisson, data = dat ) ## (B) ZIP：上と同じカウント部 + ゼロ過剰部（まずは切片のみ） m_zip0 <- glmmTMB( count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site), ziformula = ~ 1, # ゼロ過剰確率πは一定（まずはこれが基本） family = poisson, data = dat ) ## (C) ZIP：ゼロ過剰部にも説明変数を入れる例（例：mined と DOY でπを動かす） m_zip1 <- glmmTMB( count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site), ziformula = ~ mined + DOY, family = poisson, data = dat )
結果の要約と比較
summary(m_pois) summary(m_zip0) summary(m_zip1)
Family: poisson ( log )
Formula: count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site)
Data: dat
AIC BIC logLik -2*log(L) df.resid
1961.6 2019.7 -967.8 1935.6 631
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site (Intercept) 0.3095 0.5563
Number of obs: 644, groups: site, 23
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.699746 0.255631 -6.649 2.95e-11 ***
minedno 2.404469 0.330751 7.270 3.60e-13 ***
cover -0.126933 0.163003 -0.779 0.4361
DOP -0.001835 0.043020 -0.043 0.9660
Wtemp -0.045346 0.057511 -0.788 0.4304
DOY 0.114028 0.039513 2.886 0.0039 **
sppPR -1.386280 0.215165 -6.443 1.17e-10 ***
sppDM 0.230521 0.128889 1.789 0.0737 .
sppEC-A -0.770115 0.171054 -4.502 6.73e-06 ***
sppEC-L 0.621177 0.119308 5.207 1.92e-07 ***
sppDES-L 0.679165 0.118127 5.749 8.95e-09 ***
sppDF 0.080045 0.133440 0.600 0.5486
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Family: poisson ( log )
Formula: count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site)
Zero inflation: ~1
Data: dat
AIC BIC logLik -2*log(L) df.resid
1797.9 1860.5 -885.0 1769.9 630
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site (Intercept) 0.2785 0.5278
Number of obs: 644, groups: site, 23
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.33105 0.26116 -5.097 3.46e-07 ***
minedno 2.29217 0.32543 7.043 1.88e-12 ***
cover -0.16607 0.15729 -1.056 0.2911
DOP 0.10768 0.04622 2.330 0.0198 *
Wtemp -0.09543 0.06016 -1.586 0.1127
DOY 0.19485 0.04301 4.530 5.90e-06 ***
sppPR -1.21577 0.24299 -5.003 5.63e-07 ***
sppDM 0.29333 0.13534 2.167 0.0302 *
sppEC-A -0.38907 0.21558 -1.805 0.0711 .
sppEC-L 0.66579 0.12604 5.283 1.27e-07 ***
sppDES-L 0.65660 0.12358 5.313 1.08e-07 ***
sppDF 0.13844 0.14430 0.959 0.3374
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8277 0.1537 -5.386 7.19e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Family: poisson ( log )
Formula: count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site)
Zero inflation: ~mined + DOY
Data: dat
AIC BIC logLik -2*log(L) df.resid
1768.9 1840.4 -868.5 1736.9 628
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site (Intercept) 0.03968 0.1992
Number of obs: 644, groups: site, 23
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.29319 0.24838 -1.180 0.23783
minedno 1.32636 0.22614 5.865 4.49e-09 ***
cover -0.23861 0.08390 -2.844 0.00445 **
DOP 0.10218 0.04391 2.327 0.01995 *
Wtemp -0.12282 0.05730 -2.144 0.03207 *
DOY 0.23787 0.04447 5.349 8.85e-08 ***
sppPR -1.27956 0.24221 -5.283 1.27e-07 ***
sppDM 0.23481 0.13761 1.706 0.08795 .
sppEC-A -0.35685 0.22410 -1.592 0.11130
sppEC-L 0.62263 0.12741 4.887 1.03e-06 ***
sppDES-L 0.61692 0.12516 4.929 8.26e-07 ***
sppDF 0.05768 0.14656 0.394 0.69391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9083 0.2613 3.477 0.000508 ***
minedno -1.9284 0.2974 -6.484 8.93e-11 ***
DOY 0.3354 0.1278 2.625 0.008673 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AICで比較（小さいほど良い：あくまで相対比較）
AIC(m_pois, m_zip0, m_zip1)
df AIC m_pois 13 1961.625 m_zip0 14 1797.948 m_zip1 16 1768.919
AICの結果は、ポアソン単体よりZIPの方が明確に良いことを示している。
まず m_pois（通常のポアソン）のAICが 1961.6 に対し、m_zip0（ゼロ過剰は切片のみ）が 1797.9 で、約164ポイント改善している。AIC差がここまで大きいのは、データの0過剰（あるいはそれに近い構造）をポアソンが全く捉えられておらず、ZIPを入れると尤度が大きく改善したことを意味する。
さらに m_zip1（zi部に共変量を入れる）が 1768.9 で、m_zip0から さらに約29ポイント改善している。自由度（パラメータ数）は増えている（df 14→16）が、それを補って余りある改善なので、「ゼロ過剰確率（π）は一定ではなく、説明変数（ここでは mined と DOY）で動いている」というモデル化が有効だった、という解釈になる。
要するに、この3つの比較だけでも「ゼロ過剰を明示的に入れることが効いている」「しかもゼロ過剰の程度は条件で変動している」という方向性がかなり強く支持されている。
対数尤度とパラメータ数
logLik(m_pois); npar_pois <- attr(logLik(m_pois), "df")
logLik(m_zip0); npar_zip0 <- attr(logLik(m_zip0), "df")
logLik(m_zip1); npar_zip1 <- attr(logLik(m_zip1), "df")
cat("npar: pois=", npar_pois, " zip0=", npar_zip0, " zip1=", npar_zip1, "\n")
'log Lik.' -967.8123 (df=13) 'log Lik.' -884.9738 (df=14) 'log Lik.' -868.4593 (df=16) npar: pois= 13 zip0= 14 zip1= 16
まずポアソン（df=13）の logLik が -967.8123で、ZIP（zi一定；df=14）は -884.9738に上がっている。差は (-884.9738)-(-967.8123)=82.8385) で、尤度が大きく改善している（負の値が0に近づくほど当てはまりが良い）。AICが大きく下がった主因はここである。
さらにZIP（ziに共変量；df=16）は -868.4593で、zi一定のZIPから (-868.4593)-(-884.9738)=16.5145) だけ改善している。パラメータは2つ増えているが（14→16）、それ以上に尤度が伸びているので、AICでも改善として残った、という関係である。
npar: pois=13 zip0=14 zip1=16 は、モデルが順に複雑化していることの確認であり、にもかかわらず logLik が連続的に改善している。結論として、ゼロ過剰成分を入れること自体が強く効いており、さらにゼロ過剰確率を説明変数で動かすことにも追加の価値があった。
LRT（Likelihood Ratio Test; 尤度比検定）
LRT（Likelihood Ratio Test; 尤度比検定）は、「単純なモデル（帰無モデル）」と「それを包含する複雑なモデル（対立モデル）」のどちらがデータに合うかを、尤度（logLik）の差で判定する検定である。ここで重要なのは、2つのモデルがネスト（nested）、つまり「単純モデルが複雑モデルの特殊ケース」として表せる関係にあることである。
anova(m_pois, m_zip0)
Data: dat
Models:
m_pois: count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site), zi=~0, disp=~1
m_zip0: count ~ mined + cover + DOP + Wtemp + DOY + spp + (1 | site), zi=~1, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_pois 13 1961.6 2019.7 -967.81 1935.6
m_zip0 14 1798.0 1860.5 -884.97 1770.0 165.68 1 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
この比較は、通常のポアソン混合モデル（m_pois）と、同じ説明変数・同じランダム効果を保ったまま「ゼロ過剰」を追加したモデル（m_zip0）を尤度比検定で比べた結果である。m_poisはzi=~0でゼロ過剰成分を持たず、m_zip0はzi=~1としてゼロ過剰確率を切片のみで推定する。結果としてAICは1961.6から1798.0へ大きく低下し、logLikも-967.81から-884.97へ改善している。尤度比の統計量はChisq=165.68（自由度1）で、p値は2.2e-16未満と極めて小さい。したがって、ポアソンだけでは説明できない「追加のゼロ」を導入することがモデル適合を大幅に改善しており、データには構造的ゼロが混在している可能性が強く示唆される。ただし、帰無仮説（π=0）が境界にある点には留意し、AICの改善や残差診断も併せて妥当性を確認すべきである。
「ゼロをどれだけ再現できているか」を簡易にチェック
このコードは、ゼロが多いデータに対して各モデルが「ゼロの出やすさ」をどれくらい再現できているかを確認する処理である。まず各モデルから予測値を取り出すが、ZIPでは単に平均の予測値を見るだけではゼロの確率を直接得られないため、モデルが内部で持っている二つの成分を別々に取り出して計算している。具体的には、カウントの発生しやすさを表す部分と、そもそもゼロになりやすい層に入る確率を表す部分をそれぞれ推定値として取得し、それらを元の確率・平均の尺度に戻す。その上で、ZIPの定義に従って「ゼロになる確率」を各観測ごとに算出する。ポアソンモデルについては、平均が分かればゼロ確率が一意に決まるので同様の手順は不要である。最後に、実データで観測されたゼロ比率と、モデルが予測したゼロ確率の平均を並べて、どのモデルがゼロ過剰をより適切に捉えているかを比較している。
mu_pois <- predict(m_pois, type = "response")
mu_zip0 <- predict(m_zip0, type = "response")
mu_zip1 <- predict(m_zip1, type = "response")
# 予測平均から「0確率」を厳密に出すには π と λ の両方が必要なので、
# まずモデルから線形予測子を取って計算する。
## 4-1) ZIP（切片のみ）の予測0確率： P(Y=0)=π+(1-π)*exp(-λ)
eta_cond_zip0 <- predict(m_zip0, type = "link") # log(λ)
lambda_zip0 <- exp(eta_cond_zip0)
eta_zi_zip0 <- predict(m_zip0, type = "zlink") # logit(π)
pi_zip0 <- plogis(eta_zi_zip0)
p0_zip0 <- pi_zip0 + (1 - pi_zip0) * exp(-lambda_zip0)
## 4-2) ZIP（ziに共変量あり）
eta_cond_zip1 <- predict(m_zip1, type = "link")
lambda_zip1 <- exp(eta_cond_zip1)
eta_zi_zip1 <- predict(m_zip1, type = "zlink")
pi_zip1 <- plogis(eta_zi_zip1)
p0_zip1 <- pi_zip1 + (1 - pi_zip1) * exp(-lambda_zip1)
## 4-3) ポアソンの0確率：exp(-λ)
eta_cond_pois <- predict(m_pois, type = "link")
lambda_pois <- exp(eta_cond_pois)
p0_pois <- exp(-lambda_pois)
## 観測0比率と、平均予測0確率（全体の0再現）を比べる
obs_zero <- mean(dat$count == 0)
cat("Observed zero proportion:", obs_zero, "\n")
cat("Pred zero (Poisson): ", mean(p0_pois), "\n")
cat("Pred zero (ZIP zi~1): ", mean(p0_zip0), "\n")
cat("Pred zero (ZIP zi~x): ", mean(p0_zip1), "\n")
Observed zero proportion: 0.6009317 Pred zero (Poisson): 0.4836197 Pred zero (ZIP zi~1): 0.581975 Pred zero (ZIP zi~x): 0.5932345
観測されたゼロ比率は約0.601であり、データの6割がゼロである。これに対して通常のポアソンモデルが予測するゼロ確率は平均で約0.484にとどまり、実際よりゼロが少ない世界を仮定してしまっている。つまりポアソンはゼロの多さを十分に説明できていない。一方、ゼロ過剰ポアソン（切片のみのゼロ過剰）では予測ゼロ確率が約0.582まで上昇し、観測値にかなり近づく。さらにゼロ過剰確率に共変量を入れたモデルでは約0.593となり、観測ゼロ比率にいっそう近い。したがって、ゼロ過剰成分を導入することがデータの特徴に合っており、ゼロの出方が一定ではなく条件によって変わる可能性も示唆される。ただし、全体平均としてのゼロ比率が合うことは最低限のチェックにすぎず、種別や地点別など条件付きの再現性、残差診断、他の分布（負の二項やゼロ過剰負の二項）との比較も併せて判断するのが妥当である。
係数の見方
fixef(m_zip0)
Conditional model:
(Intercept) minedno cover DOP Wtemp DOY sppPR sppDM sppEC-A sppEC-L sppDES-L sppDF
-1.33105 2.29217 -0.16607 0.10768 -0.09543 0.19485 -1.21577 0.29333 -0.38907 0.66579 0.65660 0.13844
Zero-inflation model:
(Intercept)
-0.8277
m_zip0の条件付きモデルの係数は、ゼロ過剰に入らなかった場合に「観測される個体数の大きさ」をどう動かすかを示している。切片が負であることは、基準条件における期待カウントが小さいことを意味する。minednoが大きく正であるため、採掘影響のない地点では個体数が増える方向に強く働く。coverが負なので隠れ場所量が多いほど観測数が減る方向で、これは生態学的には一見逆にも見えるため、観測条件や他変数との関係を疑う余地がある。DOPは正で、降雨から日数が経つほど観測数が増える傾向を示す。Wtempは負で、水温が高いほど減る方向である。DOYは正なので季節が進むほど増える。種効果は基準種との差で、PRは大きく負で特に少なく、EC-LやDES-Lは正で多い。ゼロ過剰側は切片のみで負なので、全体として構造的ゼロに入りやすいとは言いにくいが、種差はすべて条件付き側に吸収されている解釈になる。
fixef(m_zip1)
Conditional model:
(Intercept) minedno cover DOP Wtemp DOY sppPR sppDM sppEC-A sppEC-L sppDES-L sppDF
-0.29319 1.32636 -0.23861 0.10218 -0.12282 0.23787 -1.27956 0.23481 -0.35685 0.62263 0.61692 0.05768
Zero-inflation model:
(Intercept) minedno DOY
0.9083 -1.9284 0.3354
m_zip1では、条件付きモデルの係数は同様に「ゼロ過剰に入らなかった場合の個体数」を動かすが、m_zip0より切片が大きく、minednoの効果は小さくなっている。これは、m_zip0で条件付き側が背負っていた「ゼロの説明」の一部を、m_zip1ではゼロ過剰側が引き受けた結果として理解できる。DOPが正、Wtempが負、DOYが正という方向性は維持され、種差も概ね同様でPRが特に少ない。注目点はゼロ過剰モデルで、切片が正なので基準条件では構造的ゼロに入りやすい状態が想定される一方、minednoが大きく負であるため採掘影響のない地点では構造的ゼロに入りにくくなる。さらにDOYが正なので季節が進むほど構造的ゼロに入りやすくなる。つまりm_zip1は、個体数の増減だけでなく「そもそもゼロが発生する仕組み」自体が地点条件と季節で変わる、というストーリーを与えている。
DHARMaで残差診断
if (requireNamespace("DHARMa", quietly = TRUE)) {
library(DHARMa)
sim_pois <- simulateResiduals(m_pois, n = 500)
sim_zip0 <- simulateResiduals(m_zip0, n = 500)
sim_zip1 <- simulateResiduals(m_zip1, n = 500)
plot(sim_pois); testZeroInflation(sim_pois)
plot(sim_zip0); testZeroInflation(sim_zip0)
plot(sim_zip1); testZeroInflation(sim_zip1)
}
QQプロットのKS検定が有意で、残差分布が理想的な一様性から外れている。外れ値検定も有意で、モデルが一部の観測を強く取りこぼしている。分散検定は有意ではなく、過分散そのものは決定的ではないが、予測値に応じて残差が系統的に偏る兆候があり、ポアソン単体では当てはまりが不十分だと判断できる。
ゼロ過剰検定でp値が0となり、観測されたゼロの数がモデルが想定するゼロの数と大きく食い違っている。赤線（当てはめモデル）がシミュレーション分布の端に位置しており、ポアソンモデルがゼロの多さを説明できていないことを示す。したがってゼロ過剰成分や別分布の導入が必要だと言える。
KS検定は有意で、残差はまだ完全には理想に一致しないが、001より逸脱は弱まっている可能性がある。分散検定は有意でなく、分散の不一致は目立たない。一方で外れ値検定は有意で、極端な観測が残りやすい。予測値との関係でも赤い曲線が出ており、説明変数や分布形の追加余地が残る。
ゼロ過剰検定のp値が高く、観測されたゼロの数はモデルが想定する範囲に入っている。赤線がシミュレーション分布の中心付近に位置し、ゼロの再現という点ではZIP（切片のみのゼロ過剰）が十分に効いていることを示す。ゼロの問題は概ね解消したので、次は残差の形や外れ値の要因を点検すべきである。
KS検定と分散検定がともに強く有意で、残差分布と分散構造の両方に系統的な不一致が残っている。外れ値検定も有意で、極端値が説明しきれていない。ゼロ過剰は後の検定で解消しているため、ここでの問題は主に分散の取り扱いにあり、負の二項やゼロ過剰負の二項などへの拡張が候補になる。
ゼロ過剰検定のp値が高く、観測ゼロの数はモデルの期待範囲と整合している。赤線も分布の中心付近にあり、ゼロの再現性は良好だと解釈できる。ただし005で分散の不一致が強く示されているため、ゼロの問題は解けてもカウントのばらつき自体はまだ説明不足であり、分布選択の見直しが必要である。
まとめ
6つの図は、モデルを段階的に拡張することで「ゼロの多さ」と「ばらつきの大きさ」がそれぞれどこまで説明できたかを示している。最初のポアソンモデルでは、残差分布の逸脱と外れ値が目立ち、特にゼロ過剰検定で観測ゼロが想定より大幅に多いことが明確になる。次にゼロ過剰ポアソンを導入すると、ゼロ過剰検定は非有意となり、ゼロの再現性は大きく改善する一方、残差の形のゆがみや外れ値は完全には解消しない。さらにゼロ過剰確率に共変量を入れても、ゼロの数そのものは整合的なままだが、別の図では分散の不一致が強く検出され、ゼロ以外のカウント部分に過分散が残っている可能性が示唆される。要するに、このデータではゼロ過剰成分は有効だが、それだけでは十分ではなく、カウントのばらつきをより柔軟に扱う分布や構造の検討が必要だと分かる。
ZIPによってゼロ過剰の核心には到達したが、最終的な当てはまりとしては改善の余地が大きい段階である。次の一手としては、ゼロ過剰負の二項（ZINB）への拡張、負の二項＋必要ならゼロ過剰の併用、あるいは種やサイトに関するランダム効果の見直し（ランダムスロープ等）を検討することになる。
