So-net無料ブログ作成
  • ブログをはじめる
  • ログイン

GLMMの説明用サンプル(ポアソン分布) [統計]

設定

GLMMの説明用サンプル。以下のようなデータを想定。

  • 目的変数は、ある植物の個体の種子生産数
    • ポアソン分布する。
    • 薬剤Xの投与量に応じて変化する。
      • 無処理の場合の平均値は3。
      • 圃場ごとに平均値に差がある。その標準偏差は0.4。
      • 薬剤X 1単位あたり、種子生産数が1.2倍になる。2単位では、1.2×1.2=1.44倍。
  • 独立変数は、薬剤Xの投与量
  • 実験設定
    • 圃場の6ヶ所(ブロック)に、薬剤Xを投与して種子生産数の変化を調べる
    • 薬剤Xの投与量は、0.0, 0.5, 1.0, 1.5, 2.0, ……, 10.0(単位)

データの生成


以上の設定にあうサンプルデータを生成する。

set.seed(6)

num.block <- 6
num.rep <- 1
treat <- seq(0, 10, 0.5)
inter <- log(3)		# 1.098612
fixef <- log(1.2)	# 0.1823216
ranef <- rnorm(num.block, 0, 0.4)

x <- rep(treat, num.rep * num.block)

p <- c()
y <- c()
for (i in 1:num.block) {
	p <- c(p, rep(i, length(treat) * num.rep))
	r <- ranef[i]
	y <- c(y, rpois(length(treat) * num.rep,
					exp(inter + fixef * treat + r)))
}

data <- data.frame(p = as.factor(p), x, y)
以下のようなデータが生成される。
> head(data)
  p   x y
1 1 0.0 1
2 1 0.5 3
3 1 1.0 0
4 1 1.5 0
5 1 2.0 3
6 1 2.5 0
> tail(data)
    p    x  y
121 6  7.5 31
122 6  8.0 39
123 6  8.5 38
124 6  9.0 39
125 6  9.5 44
126 6 10.0 46
グラフを描く。
xx <- seq(0, 10, 0.1)
plot(y ~ x, data=data, las=1, type="n")
for (i in 1:num.block) {
	points(y ~ x, data=data[data$p == i,], col=i, ps=1.5)
	g <- glm(y ~ x, family=poisson, data=data[data$p == i,])
	lines(xx, exp(g$coefficients[1] + g$coefficients[2] * xx), col=i, lw=2)
	print(g$coefficients)
}

解析

GLMによる解析
> data.glm <- glm(y ~ x, family=poisson, data=data)
> summary(data.glm)

Call:
glm(formula = y ~ x, family = poisson, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.1866  -2.1445  -0.3723   1.4287   4.6584  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.345637   0.068435   19.66   <2e-16 ***
x           0.186274   0.009554   19.50   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1113.93  on 125  degrees of freedom
Residual deviance:  697.82  on 124  degrees of freedom
AIC: 1175.7

Number of Fisher Scoring iterations: 5

GLMによる解析では、過分散(overdispersion)が発生する。
GLMMによる解析
glmmMLを使用。
> library(glmmML)
> data.glmm <- glmmML(y ~ x, cluster=p,
+ 				family=poisson, data=data)
> data.glmm

Call:  glmmML(formula = y ~ x, family = poisson, data = data, cluster = p) 


              coef se(coef)      z Pr(>|z|)
(Intercept) 1.1373 0.278318  4.086 4.38e-05
x           0.1863 0.009554 19.496 0.00e+00

Standard deviation in mixing distribution:  0.6592 
Std. Error:                                 0.1837 

Residual deviance: 181.1 on 123 degrees of freedom 	AIC: 187.1 

nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:学問

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0