random dispersal

スポンサーサイト

上記の広告は１ヶ月以上更新のないブログに表示されています。新しい記事を書く事で広告が消せます。

# 生態学会で話していて、過分散データの取り扱いが話題に上った。負の二項分布を使えという意見と、GLMMを使えという意見と両方があるが、実際どっちの方がよいというのはあるのだろうか？（負の二項分布でGLMMという意見もあるでしょうが…そんな複雑な推定はMCMCした方がいい気がする…）

# 過分散 overdispersion：ポアソン分布や二項分布のモデルに当てはめる場合、これらの分布型が平均と分散が表裏一体なのだが（ポアソンでは平均=分散)、実際のデータは大抵はそうでないので分散が過剰となることが多い、という問題。
# GLMのsummary出力結果で、下の方に、こんなのがあって…
# Residual deviance: 5094.7 on 71 degrees of freedom
# …ここで、5094.7/71の値が、1.5を超えているならば、何らかの対処をした方がよいだろう（cf. モデル選択ペンギン本、Zuur et al (2009) Mixed effects models and extentions in Ecology with R. Springer。これのp. 225-6に"1.5"が目安と書かれているのを発見）

# 過分散しているデータとは下記のように裾野が片方に広く延びている（左：負の二項分布、右：正規ノイズを加えたポアソン分布）。

par(mfcol=c(1,2)) # グラフ領域をヨコに2分割
plot(density(rnegbin(1000, mu=exp(1), theta=1/4))) # 負の二項分布（平均 exp(1)、θ 1/4）
plot(density(rpois(1000, lambda=exp(1 + rnorm(1000, mean=0, sd=4^0.5)))))
# ポアソン分布（平均λ= exp(1)）に、正規分布（平均 0、分散 4）のノイズを加えた

# では、負の二項分布をしたデータと、ポアソン分布に正規乱数ノイズを加えたデータを人工的に作成し、それぞれGLM (poisson)、GLM（負の二項分布）、GLMM（データ数分の個体差があると見なす）で推定した結果を比べてみよう。

# これを実行するためのプログラムは以下：
# 平均 exp(alpha + beta*x)とし、真の値はalpha = 2、beta = 0.5、あと負の二項分布の歪みのパラメータθ= 1/4、正規乱数ノイズの分散4も設定
# （ただし、θや分散の値については今回はチェックはしない）

library(MASS) # glm.nbの呼び出し
library(lme4) # lmerの呼び出し

OD <- function(alpha, beta, siml, distribution) {
# 切片、傾き、シミュレーション回数、データ分布（1：負の二項分布；２：正規ノイズ）
estim <- rep(0,6) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(4) # 乱数の種を指定（他の値にすると推定エラーが出るかもですglm.nbやlmerの問題）
for (s in 1:siml) {
Noise <- rnorm(100, 0, 4^0.5) # ノイズの分散を4とした
x <- rep(c(0:9), each=10)
y <- rbind(rnegbin(100, mu=exp(alpha + beta*x), theta=1/4), # 負の二項分布
rpois(100, lambda=exp(alpha + beta*x + Noise))) # 普通のポアソン分布にノイズを追加
y <- y[distribution,] # データ分布の選択
ID <- 1:length(y) # データの数だけ個体差（ID）を設定
D <- data.frame(x, y, ID)
pois <- coef(glm(y ~ x, family=poisson, data=D)) # GLM (poisson)
negb <- coef(glm.nb(y ~ x, data=D)) # GLM (負の二項分布)
pois.mm <- fixef(lmer(y ~ x + (1|ID), family=poisson, data=D)) # GLMM
estim2 <- c(cbind(pois, negb, pois.mm))
estim <- rbind(estim, estim2) # simlの数だけ推定値を重ねていく
}
estim <- estim[-1,] # 1行目のイントロン行を除去
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(-2,6), ylim=c(0,1), main="alpha (intercept)")
for (i in 2:3) { lines(density(estim[,(2*i-1)]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,2]), lwd=4, xlim=c(0,1), ylim=c(0,6), main="beta (coefficient)")
for (j in 2:3) { lines(density(estim[,(2*j)]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
}

OD(alpha=2, beta=0.5, siml=100, distribution=1) # データが負の二項分布の時
#（けっこう時間が掛かる、エラーメッセージも沢山出るが無視）
# ランダム効果がN数と同じ数だけあるよ！？というエラーが大量に出てくるが、推定はちゃんとできているので無視
# Number of levels of a grouping factor for the random effects is *equal* to n, the number of observations

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# GLM (poisson)は両方向に大きく広がった。GLM（負の二項分布）はalphaもbetaも真の値に近かった。GLMMはbetaの値は真の値に近かったが分布の裾が広かった、alphaは盛大にズレた

# 次にポアソン分布に正規乱数ノイズを加えたデータを同様に推定してみる
OD(alpha=2, beta=0.5, siml=100, distribution=2)

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# GLM (poisson)は真の値から大きく外れた。GLM（負の二項分布）はalphaは盛大にズレた、betaは真の値に近かったが裾野が広かった。GLMMはalphaもbetaも真の値に近かった

# 結論として、過分散データと言っても一辺倒ではなく、負の二項分布、正規ノイズ、どちらに近いかによって結果が変わってくるという、ごく当然の結果となった。glm.nb、GLMMのどちらがよいかは、結局データ依存と言うことですね。選択基準としては…glm.nb、GLMMでAICのより小さい方を選ぶしかないのだろうか…。
# なお、敢えて言うなら、GLMMの方がいずれのケースでもbetaの推定がマシな感じですね。切片alphaのズレは単にそれぞれの統計モデルが異なるからだということになるのだろうか（こっちのモデルを想定するならばalpha=...、他方ならalpha=...になる、みたいな）。

# もっとも、ポアソン分布の場合は今回のような選択肢があるけれど、二項分布の場合はGLMMを使うしかない（怪しげなQAICの使用を回避するならば）。また、将来的にはMCMC推定が全盛になるだろうと見越して、総じて考えれば、あえてどちらか選ぶならGLMMなのかな。

# P.S. ゼロデータが多い場合にはzero-inflated modelsがいいようだが、私はまだ勉強不足。それについてはまたいずれレポートしたい。

テーマ:研究者の生活 - ジャンル:学問・文化・芸術

R合宿・第五部：現代統計への導入（GLM、GLMMの基礎）

R合宿の統括記事はこちら

# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）

### GLM準備編 ###########
# まずは古典統計をGLM風に解いてみよう
# 回帰の例（yへのxの影響が有意か、ではなく、yへのxの影響の仕方を推定する）
x00 <- c(1:100)
y00 <- rnorm(100, mean=x00) # x00の値を元に正規乱数を発生
m00 <- lm(y00 ~ x00) # 古典統計編でも使ったlmという推定関数
summary(m00)
# Coefficients: # 切片、回帰係数がそれぞれゼロから有意に異なっているかを示す（いわゆる検定ではない）
# Estimate Std. Error t value Pr(>|t|) # 推定値、SE、t値、有意確率
# (Intercept) 0.325584 0.186879 1.742 0.0846 . # 切片
# x00 0.996002 0.003213 310.016 <2e-16 *** # x00の回帰係数（約1）

plot(y00 ~ x00)
curve(coef(m00)[1] + coef(m00)[2]*x, add=T)
# y00 = x00 + α くらいの一次直線になるはず（切片の値は正負にばらつく）
# こうやってデータ（Y）と説明変数（X）の関係を説明するために、係数を求めるのがGLMの基本形

# ANOVAの例
F01 <- factor(rep(c("A","B"), each=20))
y01 <- rnorm(40, mean=rep(c(10,20), each=20)) # A, Bそれぞれ平均10, 20になるように設定する
m01 <- lm(y01 ~ F01)
summary(m01)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 10.2747 0.2316 44.36 <2e-16 *** # 切片=F01のAの推定値（約10）
# F01B 10.0316 0.3276 30.62 <2e-16 *** # 切片にこれを足すとF01のBの推定値（約20）
# F01という要因がAかBかによって、このぐらいyの値が異なるという関係を推定した
#（Aのとき10、Bのとき20という、真の値に近い推定が得られた）
# こういう、推定というプロセスは普通のANOVAではやらないので馴染みがないでしょうね。

plot(y01 ~ F01) # 図示する
points(1, coef(m01)[1], cex=3) # 推定値をプロットしてみる
points(2, coef(m01)[1] + coef(m01)[2], cex=3) # Bの方はこういう足し合わせた値になる

######## 乱数で遊ぼう（GLM準備） ######
# 一般化線形モデルの内部で計算してくれていることを理解するための演習、
# けっこう理解は難しいかも…ここは雰囲気を掴んでもらえるだけでも十分かなと思います
sample(1:10, replace=T) # 1~10から適当に数を選択（replace=T→同じ値を繰り返し取り得る）
rpois(10, 10) # 平均10の整数 x 10個（前の数字は個数、後ろの数字は平均lambda（ラムダ）を表す）
rnorm(10, 0, 1) # 平均0、標準偏差1の実数 x 10個
# 適当に値を変えて、コピペ＆リターンを繰り返してみよう。
# パラメータは同じ値のままでも、毎回の結果は変わります（それが乱数というもの）

### 今度はサイコロ
# まず、各目の出る確率はいずれも 1/6 = 0.1666...
# 1回振って、ある目（"6"とします）が出たら"1"その他なら"0"となる
rbinom(1, 1, 1/6)
# これを10回繰り返し、平均を取ってみる＝出る確率
mean(rbinom(10, 1, 1/6)) # 1/6にはならない
# 繰り返しを増やしてみる（サイコロ10万笑）
mean(rbinom(100000, 1, 1/6)) # ほぼ 0.166 (1/6)
#（繰り返しを増やすと真の値に近くなる：大数の法則、という性質）

### 図示してみよう
hist(rnorm(10000, mean=5, sd=2)) # 正規分布（実数）、平均と分散で分布が決まる分布
hist(rnorm(10000, mean=100, sd=2)) # 平均の大きさにかかわらず、いつも左右対称な性質を持っている
hist(rpois(10000, lambda=5)) # ポアソン分布（正の整数）、lambda: 平均かつ分散（平均=分散、な分布）
hist(rpois(10000, lambda=100)) # 平均値を大きくすると、左右対称に近くなる（正規分布に近くなる）
hist(rbinom(10000, size=10, prob=0.3)) # 二項分布（0~上限の整数）、size：上限数（上限と確率で分布が決まる）
hist(rbinom(10000, size=100, prob=0.3)) # 上限数（試行回数）を大きくすると、やはり左右対称（
正規分布に近くなる)

##### 尤度（そのデータが生じる確率）を計算するツール ########
po <- rpois(10000, 1) # 平均1のポアソン分布に従う乱数を10000個発生させ、poに格納する
length(po[po==1])/10000 # 上記のうち、実際に値1を取った乱数の個数を全個数10000で割る
dpois(1, 1) # 1という値が、平均1のポアソン分布で生じる確率
# 両者の値は、ほぼ同じになるはず。rpois、dpoisの計算していること、分かりますか？？

# もう少し進めていきます
dpois(1:10, 5) # 1~10というデータが、それぞれ平均5のポアソン分布で生じる確率を計算
# dpois(1,5); dpois(2,5); dpois(3,5) ;;;; dpois(9,5); dpois(10,5)# を一括で計算するのと同じ意味
hist(rpois(10000, 5), breaks=c(0:16)) # 平均5のポアソン分布を図示
lines(1:10, dpois(1:10, 5)*10000, type="b") # 確率を重ねて図示（平均である約5で最大となる）
sum(dpois(1:10, 5, log=T)) # 扱いやすいよう対数を取り合計するのが一般的（掛け算のlogは足し算になる）
# これが対数尤度、これを最大化させるようなパラメータ（平均値など）を推定する

###### パラメータ推定をシンプルな例で試してみます ########
# まず、ツールの紹介
curve(-(x^2 - 4*x + 4), xlim=c(-10,10)) # 高校1年くらいで習うあの曲線
abline(v=2); abline(h=0) # 見やすくするよう、縦横のラインを追加
Ys <- function(x) -(x^2 - 4*x + 4) # -(x - 2)^2 = 0 を解くと、x = 2 で最大でしたね？
optimize(Ys, interval=c(-5, 5), maximum=T) # これを計算してくれるツール（区間は指定要）
# X=2のとき、最大値0を取ることが確認できましたか？（二行上の計算と同じ）

# これを用いて推定してみます（尤度計算に即した、もう少し複雑な例へ拡張）
# ちょっと単純化して Y = b*X を考えてみる（切片なしの回帰式）
b <- 5 # bの真の値は5とする
X <- c(1:100) # Xは1~100
Y <- rpois(100, exp(b*X)) # 平均 exp(b*X)、なぜ指数を取るのかは今はスルーします
poi <- function(b) sum(dpois(Y, exp(b*X), log=T))
# bの値を変化させたときの反応を見るため、関数にしてみる（上記のx^2 - 4*x + 4の状態）
# ちなみにlog=TはXの各データの確率の対数を取るという意味
# なぜ？って、同時に各データが生じる確率は掛け算になるわけだけど（Π、パイの大文字）、
# logを取れば足し算に変わるから扱いやすくなるから（対数の性質、大丈夫ですか？）
optimize(poi, interval=c(0,10), maximum=T) # 対数尤度が最大となるbの値を推定
# 便宜上、推定するためのxの範囲は指定しなければならない（0~10にしてみた）
# $maximumの値：poiが最大値となるxの値。ちゃんと真の値5にほぼ一致しているはず
# $objectiveの値：こちらもほぼ0に近い値になっているはず（確率値の対数は0で最大）

############ GLM 実行に用いるサンプルデータを用意する
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 10
X <- rep(c(0:9), each=10)
f1 <- rep(rep(c(0:3), each=5), 5) # 予備操作
F1 <- rep(rep(c("C","T1","T2","T3"), each=5), 5) # f1の0~3をC~T1,2,3に読み替え
Y1 <- rpois(100, lambda=exp(-2 + f1 + 0.5*X)) # 切片-2（=C）, Tは1ずつ増加、傾き0.5、の乱数発生
Y2 <- rbinom(100, N, prob=logistic(-10 + f1 + 2*X)) # logisticの中身を基にして、二項乱数を発生
vol <- 10*(X+50) + abs(rnorm(length(X))) # 下記のoffset編で用います。数字に特に意味はなし
Dg <- data.frame(X, F1, Y1, Y2, vol)

############ GLM（ポアソン、ANCOVAタイプ） #######
curve(exp(3*x+1), lwd=4) # ポアソン分布のGLMで推定される曲線を図示
# exp()内の値を変えていじって、曲線の性質に慣れてみよう

# 上で定義した値を用い、切片-2, treat+1, 傾き0.5、に戻れるか否かを試してみる
m1 <- glm(Y1 ~ X + F1, family=poisson(link="log"), data=Dg)
summary(m1) # 解析結果の全貌

# Call:
# glm(formula = Y1 ~ X + F1, family = poisson(link = "log"), data = Dg)
#
# Deviance Residuals:
# Min 1Q Median 3Q Max
# -2.0815 -0.8012 -0.3862 0.6185 3.0358
#
# Coefficients:
# Estimate Std. Error z value Pr(>|z|) # 切片と回帰係数の推定値、SE、z値（Wald統計量）、P値
# (Intercept) -2.32878 0.16594 -14.034 < 2e-16 *** # 切片の真の値は-2
# X 0.51697 0.01138 45.436 < 2e-16 *** # Xの係数の真の値は0.5
# F1T1 1.18029 0.16677 7.077 1.47e-12 *** # 真の値は1
# F1T2 2.15645 0.15124 14.258 < 2e-16 *** # 真の値は2
# F1T3 3.14441 0.14816 21.222 < 2e-16 *** # 真の値は3…に近い値となりましたか？
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # 各要因が独立の時、ゼロと異なる有意確率
#
# (Dispersion parameter for poisson family taken to be 1) # ポアソン分布なので平均=分散になるべきことを意味
#
# Null deviance: 6470.285 on 99 degrees of freedom
# Residual deviance: 96.548 on 95 degrees of freedom # ここが1対1に近ければ推定はたぶん大丈夫
# AIC: 442.39 # このモデルのAIC値
#
# Number of Fisher Scoring iterations: 5
#
# どうだろう？おおよそ上記の値を推定できているのではないだろうか？
# 上記の値は個別に取り出すことが可能
coef(m1) # 回帰係数。coef(model1)[1]は切片、[]内の数字で抜き出すものを指定可能
logLik(m1) # 最大化対数尤度、さっき上で求めたものとはちょっと違う（詳細は省略）
deviance(m1) # 最大化対数尤度 x -2。小さいほどデータにfitしてることを示す
AIC(m1) # 460.1022。deviance + 2*パラメータ数 ←複雑なモデルほど罰則が追加される
# 乱数を使っているので、実行するたびに値は多少変わるはず

# 比較対照：切片だけのモデル
m0 <- glm(Y1 ~ 1, family=poisson(link="log"), data=Dg)
# ヌルモデル（null model）という。何も要因を入れていない切片のみモデル
AIC(m0) - AIC(m1) # ΔAIC = 6544.604、m1による説明力の向上
# Deviance explained:
1 - deviance(m1)/deviance(m0) # 0.9842949
# "いちおう"データの98.4%がモデルで説明できた！？
# R2決定係数に相当。ただしR2ほど明快な理屈はないのでR2ほど市民権はないようだ

# グラフの図示
plot(Y1 ~ X, cex=2, pch=21, bg=c(1:4)[F1], data=Dg)
curve(exp(coef(m1)[1] + coef(m1)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示（F1のCは切片の値と等しくなっている。T1~3はCからの増分の値）
curve(exp(coef(m1)[1] + coef(m1)[2]*x + coef(m1)[3]), lwd=5, col=2, add=T)
# 二本目、図示（T2）
for (g in 2:4) # 二本目以降を一括で図示
curve(exp(coef(m1)[1] + coef(m1)[2]*x + coef(m1)[g+1]), lwd=5, col=g, add=T)

############# 注：GLMでの、ポアソン分布や二項分布を使用上の注意 #######
# summary(m1)の一番下の方、Residual deviance: 104.55 on 95 degrees of freedom
# この数字が1:1から大きく外れてると推定が怪しいという目安になる（Crawley 2005）＊。
# 左 > 右の場合、overdispersion（過分散）という。
# 過分散はGLMMの利用で軽減できる
# ＊しかし、この判定方法は標本数もかなり多いときにのみ有効
# Crawley MJ (2005) Statistics: an introduction using R. John Wiley & Sons, West Sussex

# cf. 過分散する例
data(CO2)
m11 <- glm(conc ~ uptake + Plant, family=poisson, data=CO2)
# CO2濃度が植物の吸収量によって変化するという温暖化対策になりそうな解析例（苦笑）
summary(m11)
# かなりの過分散になる Residual deviance: 5094.7 on 71 degrees of freedom
# GLMMなどを利用して、対処が必要

########## GLM（ポアソン、offset）#######
# 単位量あたりのy1の解析、密度のようなもの。割り算するべからず。
m2 <- glm(Y1 ~ X + F1 + offset(log(vol)), family=poisson(link="log"), data=Dg)
summary(m2)
plot(Y1/vol ~ X, cex=2, pch=21, bg=c(1:4)[F1], ylab="Density", data=Dg)
curve(exp(coef(m2)[1] + coef(m2)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示
for (g in 2:4) # 二本目以降を一括で図示
curve(exp(coef(m2)[1]+ coef(m2)[2]*x + coef(m2)[g+1]), lwd=5, col=g, add=T)
# ちゃんとプロットを貫くジャストの曲線が描けましたね？

# cf. 成長量（after/before）のようなデータはガンマ分布（正の値の連続量）で同様に解析する。
# やはり割り算すべきではないということ。単に変化率だけでなく、絶対量も評価するため。
# m22 <- glm(after ~ X + F1 + offset(log(before)), family=Gamma(link="log"), data)
# あとのやり方はポアソンと同様

########### GLM（二項分布）ロジスティック回帰型 ########
logistic <- function(xx) 1 / (1 + exp(-xx))
curve(logistic(1*x-5),xlim=c(0,10)) # logisticの数値を変えて図示してみよう

m3 <- glm(cbind(Y2, N-Y2) ~ X + F1, family=binomial(link="logit"), data=Dg)
# cbind()の部分はこう書くしかないようです、ちょっと煩雑ですが
summary(m3)
# Residual deviance: 54.02 on 95 degrees of freedom
# under dispersionのようです（実例ではあまり出ない）。少し推定がよくないはず
# 二項分布ではなるべくGLMMがお勧め（とくにデータ数が二桁程度の少数の場合は）
# cf. Warton & Hui (2011) Ecology 92: 3-11
# N=20以下の過分散サンプルでGLMの推定がひどいことを示す例あり

plot(Y2/10 ~ X, cex=2, pch=21, bg=c(1:4)[F1], ylim=c(0,1), data=Dg)
curve(logistic(coef(m3)[1] + coef(m3)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示
for (g in 2:4) # 二本目以降を一括で図示
curve(logistic(coef(m3)[1] + coef(m3)[2]*x + coef(m3)[g+1]), lwd=5, col=g, add=T)
# 図示すると、当てはまりは悪くないように見える

############ GLM（変数が多い重回帰分析タイプ） ######
# まずは使用するデータの下準備
library(car)
data(Chile) # carパッケージに入ってるChileというデータを使用
head(Chile) # Chileの中身を眺める

### 各連続変数の関係を図示してみる（対戦表のような配置、各変数をx, yとしたときの散布図）
pairs(statusquo ~ income + population + age, panel=panel.smooth, data=Chile)
# 注意すべきは、相関のある説明変数があるか否か（右上がり：incomeとpopulationは怪しい）

### 回帰木によって交互作用の有無を判断（枝分かれの構造によって交互作用が視覚的に見やすい）
library(rpart)
res <- rpart(statusquo ~ income + population + age, method="anova", data=Chile)
# デフォルトは"anova"（正規）。その他、"poisson"（ポアソン）、"class"（二項分布、多項分布）
print(res) # 数値結果
par(xpd=T) # 図をちゃんと収めるための準備
plot(res) # 図示
text(res) # データの図示
# 値いくつ以上のデータがdevianceの説明にどれくらい効いてるかなどを図示する
# 数値結果の階段構造は分岐を表す
# 高さが長い＝deviance explainedが大きい

#### statusquo（満足度？）をincome, population, ageで回帰してみよう
# その前にちょっと説明変数に細工します
# 基準化（scale: 平均ゼロ、分散1）←解析の是非には影響しないが、
# このデータで解析すれば、回帰係数が標準化される、
# すると、絶対値の大きさでその要因による影響の大きさが比較可能になる
Chile2 <- na.omit(Chile)
Chile2[,c(2,4,6)] <- scale(Chile2[,c(2,4,6)]) # income, population, ageだけscale化

### やっとモデリングします
reg2 <- glm(statusquo ~ income + population + age, family=gaussian, data=Chile2)
summary(reg2)

library(MASS) # 雑多な基本パッケージ
stepAIC(reg2) # 簡易的なモデル選択の実行（減少法）
# stepAICは高速で便利だが、下記のlmerには使用不能だし、選択ミスも生じる場合あり

library(MuMIn) # もうひとつのモデル選択の…というより…モデル評価のためのパッケージ
dredge(reg2, rank="AIC") # 鬼のようなモデル選択（総当たり…）
# rank: AIC 値をよりどころにする； , fixed=c("X1", "X2",...): 必ず含める要因
# 変数の数が多いと階乗分、組み合わせが多くなって、かなりの計算時間となります。

# AICの値の順に上から下に向かってモデルを列挙します
# (Int): 切片； X（各要因）：回帰係数； k：パラメータ数
# delta：ΔAIC、ベストのモデルからの AIC の増分
# weigtht：Akaike weight、そのモデルが選択されうる確率
# Akaike weight：便利だが理論的背景が曖昧なので嫌われてもいる（cf. kubolog）
# cf. AICc：N数が少ないときのAICとして好む人がいるが、これは正規分布用なので注意（誤用は沢山）
#（cf.推進してる人（本） Burnham & Anderson 2002）

### 主成分分析の成分に変換すると相関関係（多重共線性）は無くなる
Chile3 <- na.omit(Chile) # NAを含む行を削除
pca <- prcomp(Chile3[,c(2,4,6)], scale=T) # 規格化する項が付いている
summary(pca)
# Importance of components:
# PC1 PC2 PC3
# Standard deviation 1.1070 1.0009 0.8790
# Proportion of Variance 0.4085 0.3339 0.2576
# Cumulative Proportion 0.4085 0.7424 1.0000

# PC2まででも、すでに74.24%の情報量がある→PC3は除いてもよさそうだ
PC1 <- pca$x[,1] # PC1
PC2 <- pca$x[,2] # PC2

plot(pca$rotation[,1], pca$rotation[,2], type="n", xlim=c(-2,2), ylim=c(-2,2)) # limは適宜調整要
text(pca$rotation[,1], pca$rotation[,2], colnames(Chile3[,c(2,4,6)]))
# やはりpopulとincomeは類似→より第一軸、age→より第二軸

reg22 <- glm(statusquo ~ PC1 + PC2, family=gaussian, data=Chile3)
summary(reg22)
AIC(reg2) - AIC(reg22) # -109.2854 このケースではPCAの使用でAICが悪化してしまった
# しかし、reg2が多重共線性を含んでいたことを考えると、一概にreg22がダメとも言えない
# けっきょく説明変数多変量のモデリングは一義的に切れないところがある、と思います

############### GLMM #######
# 一般化線形混合効果モデル（Generalized linear mixed effects model: GLMM）
logistic <- function(xx) 1 / (1 + exp(-xx))

### あらためてサンプルデータを用意します
# noise（平均0、分散4の正規分布の乱数）を追加
# noiseのせいでバラツキがひどくなるが、平均（切片や回帰係数の値）は変わらない
# 過分散データなどもノイズが多い状態と言えるだろう
# こんなノイズの多いデータでもGLMMなら"比較的"うまく推定してくれる

N <- 10
noise <- rnorm(100, 0, 4^0.5) # SD = √4
x <- rep(c(0:9), each=10)
y1 <- rpois(100, lambda=exp(-2 + 0.5*x + noise)) # 普通のポアソン分布にノイズを追加
y2 <- rbinom(100, N, prob=logistic(-10 + 2*x + noise)) # logisticの中身を基にして、二項乱数を発生
ID <- 1:length(y1)
Dmm <- data.frame(N,x,y1,y2,ID)

library(lme4) # GLMMを実行する汎用パッケージの読み込み

######### GLMM データが正の整数（ポアソン分布）の場合 ###
# まずは普通のGLMで推定
plot(y1 ~ x, cex=2, ylab="Frequency", data=Dmm)
pois1 <- glm(y1 ~ x, family=poisson, data=Dmm)
summary(pois1) # 推定値は真の値（-2, 0.5）から大きくずれてしまう
# 案の定ひどい過分散 Residual deviance: 5998.1 on 98 degrees of freedom
# 注：乱数を使用しているので値は実行するたびに多少異なるはず

# 図示してみよう：
curve(exp(coef(pois1)[1] + coef(pois1)[2]*x), lwd=8, lty=3, add=T)
# （太い点線）当てはまってるんだか、どうなんだか、図からはよく分からない

### 次に、GLMMによる推定を試みる
# lmerという関数。IDによるノイズを分散値として推定
# ID毎に推定すると、無数の切片を推定することになる→複雑すぎるモデルなので回避する
# 複雑なデータに対処しつつモデルは単純化できる。これがGLMMの使いどころ。だと思う。
pois2 <- lmer(y1 ~ x + (1|ID), family=poisson, data=Dmm) # 警告は出るが推定はうまくいってる
summary(pois2)

# Generalized linear mixed model fit by the Laplace approximation
# Formula: y1 ~ x + (1 | ID)
# Data: Dmm
# AIC BIC logLik deviance # AIC, deviance以外にもいろいろ出してくれる
# 315.3 323.1 -154.6 309.3
# Random effects: # ランダム要因の推定結果
# Groups Name Variance Std.Dev. # ここでは分散とSDだけ
# ID (Intercept) 4.4614 2.1122 # IDごとに切片を設け、そのバラツキの推定値（真の値4, 2）
# Number of obs: 100, groups: ID, 100 # IDの数
#
# Fixed effects: # 固定要因の切片、回帰係数の推定
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -2.29482 0.53798 -4.266 1.99e-05 ** # 切片の真の値は-2
# x 0.51653 0.09214 5.606 2.07e-08 *** # 回帰係数の真の値は0.5
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Correlation of Fixed Effects: # 固定要因間の相関
# (Intr)
# x -0.888

# 今度はけっこう真の値に近い推定となったはず
#（データDmmに乱数を用いているので、データを定義するたびに値は変わります）

# 再び図示してみよう
curve(exp(fixef(pois2)[1] + fixef(pois2)[2]*x), lwd=8, lty=1, add=T)
# （太い実線）xについての回帰曲線を図示
for (i in 1:nrow(Dmm))
curve(exp(fixef(pois2)[1] + fixef(pois2)[2]*x + ranef(pois2)$ID[[1]][i]), lwd=2, lty=3, add=T)
# さらに、あえてランダム効果の数値（ID毎の無数の切片の値）を取り出し図示してみる
# 普通はこういうことはやらない。だからこそ分散の値だけで割愛している。

# cf. もう一つの過分散対処法：ポアソンよりもより歪んだ整数値を取る分布 → 負の二項分布
library(MASS) # に含まれる glm.nb、歪みの程度を変化させて対処するパラメータを持つ
glmb <- glm.nb(y1 ~ x, data=Dmm) # family指定不要
summary(glmb)
# ノイズではなくデータの本来持つ性質として分布が大きく偏るならば利用したらいいと思う
# GLMMの一種とも言えるが、一般的なGLMMと違い、パラメータ数の節約にはならないのではないか？
# その点で、比較対象と言うよりは別のデータ構造を仮定したモデルだと思う
# （このデータの例ではGLMMの方が推定がよかった）

####### GLMM、二項分布での例 ###
# 図示してみる
plot(y2/N ~ x, cex=1.5, ylab="Probability", data=Dmm)
# まずは普通のGLMで推定を試みる
binom1 <- glm(cbind(y2, N-y2) ~ x, family=binomial, data=Dmm)
summary(binom1) # 真の値（-10, 2）から大きくずれた推定となる
# やはり過分散 Residual deviance: 171.34 on 98 degrees of freedom

# 図示してみる：
curve(logistic(coef(binom1)[1] + coef(binom1)[2]*x), lwd=8, lty=3, add=T)
# （太い点線）今度も、ずれていることが見た目では分かりにくいが…

### GLMM、二項分布の実行 ###
binom2 <- lmer(cbind(y2, N-y2) ~ x + (1|ID), family=binomial, data=Dmm) # 警告は出るが推定はうまくいってる
summary(binom2) # 真の値（-10, 2）により近くなったはず
curve(logistic(fixef(binom2)[1] + fixef(binom2)[2]*x), lwd=8, lty=1, add=T)
for (i in 1:nrow(Dmm))
curve(logistic(fixef(binom2)[1] + fixef(binom2)[2]*x + ranef(binom2)$ID[[1]][i]), lwd=2, lty=3, add=T)

# cf. もうひとつのGLMM関数 glmmML
# メリット：高速なモデル選択関数 stepAICが適用可能
# デメリット：binomialとpoissonにしか使えない
# ランダム項は一つしか指定できない（lmerでもなるべく一つにした方がよい）
# 推定計算自体が実行できずエラーになることもしばしば。

テーマ:研究者の生活 - ジャンル:学問・文化・芸術

R合宿・第四部：古典統計編

R合宿の統括記事はこちら

# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）

###### Rで"古典"統計 ###
# 過去との架け橋として掲載しますが、大半の現実的なデータの解析ではもはや必要悪と私個人は考えています。
# なぜか？→国内の特定分野の学術誌ではなく、英文で書かれた国際的な標準的学術誌を参考にしてください。
# 次第に、効果の推定値を重視した論文や、AIC値など情報量規準を元にした議論が増えてきているはずです。
# 多くの皆様がGLMなど統計モデリングへとステップアップしてくださることを願っています（cf. R合宿のGLM記事、統括記事）

#### 今回はRに内蔵のデータを使用します ###
# データCO2を呼び出す
data(CO2)

############ 正規性の検定（例：Shapiro-Wilk test）######
# 応答変数（= 従属変数、目的変数）が正規性を満たすかをチェック
shapiro.test(CO2$uptake) # CO2のuptake、をドルマークで繋ぐことで表現できる

# Shapiro-Wilk normality test
#
# data: CO2$uptake
# W = 0.941, p-value = 0.0007908
# 正規分布に対し有意差有り→正規性を満たしていないので、変数変換をする

CO2$uptake.L <- log(CO2$uptake + 1) # 対数変換した列をCO2に追加
CO2$uptake.R <- sqrt(CO2$uptake + 0.5) # 平方根変換した列をCO2に追加

# 再度トライする
shapiro.test(CO2$uptake.L)
shapiro.test(CO2$uptake.R)

# 結局、むしろ正規性が悪化したという珍しいケース…
# 「改善を試みたがダメだった」と、次善の策として元のまま行くしかないだろう
# ちなみに、なるべく図示して、どういうデータなのか視覚的にチェックしてみよう：
par(mfrow=c(1,3)) # 複数の図を行数、列数を指定してまとめるコード、でしたね？
plot(CO2$uptake)
plot(CO2$uptake.L)
plot(CO2$uptake.R)
# 片山状に下がっていっている感じ、たしかに正規分布ではない…

############ 等分散性の検定（例：Levene's test）######
library(car) # car というパッケージに入ってるので呼び出す
leveneTest(uptake ~ Plant, data=CO2)
# CO2というデータのuptakeをPlantの群で分割して、という意味

# Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
# group 11 0.5877 0.8328
# 72
# 今度は有意差無し→分散が異なるとは言えなかった→等分散性は否定されなかった

############ 1 way-ANOVA ######
library(car) # ←上で一度呼び出していれば不要なコードです（以下でも同様）
Anova(lm(uptake ~ Plant, data=CO2))

# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# Plant 4862.2 11 6.569 1.842e-07 ***
# Residuals 4844.8 72
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

############ post-hoc test（多重比較）######
pairwise.t.test(CO2$uptake, CO2$Plant, method="holm") # holm（sequential-Bonferroni）
# 群のペア間での差異の有意確率を対戦行列で出力する
# 手計算でも簡単にできる上、無理のない、実用的な多重比較（ノンパラメトリックでもOK）

TukeyHSD(aov(uptake ~ Plant, data=CO2)) # TukeyHSD法
# こちらはペアで有意確率（p adj）を算出、差異の推定と95%信頼区間も出してくれる
# ANOVAとともに使用（パラメトリックな方法）

############ 多元ANOVA ######
library(car)
Anova(lm(uptake ~ Type*Treatment, data=CO2))
# "*"：単効果の項：Type, Treatment；交互作用の項：Type:Treatment、の両方を表せる

# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# Type 3365.5 1 52.5086 2.378e-10 ***
# Treatment 988.1 1 15.4164 0.0001817 ***
# Type:Treatment 225.7 1 3.5218 0.0642128 .
# Residuals 5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# このAnova関数は多元の時、標準ではType II 平方和で検定する
# →単効果を優先する平方和であり、交互作用は特殊ケースであるという考え方
# もし交互作用も等評価したい場合、下記のようにしてType III平方和で検定する
Anova(lm(uptake ~ Type*Treatment, data=CO2), type="III")

# Anova Table (Type III tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# (Intercept) 26217.3 1 409.0389 < 2.2e-16 ***
# Type 924.0 1 14.4165 0.0002839 ***
# Treatment 134.6 1 2.1007 0.1511413
# Type:Treatment 225.7 1 3.5218 0.0642128 .
# Residuals 5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Type II の時とは異なり、Treatmentが有意でなくなった

# cf. 多元ANOVAは多くとも3元で留めましょう。多すぎると摩訶不思議な交互作用が出ます

############ あえて多元ANOVAの後で多重比較をしたいという場合… ######
# 交互作用が有意でなければ、各要因それぞれにつき1wayANOVAと同様に多重比較する
# （今回は交互作用が有意でなかったが、有意だった場合のやり方も示しておこう↓）
# まず、両変数をドッキングさせた新変数を作る
CO2$TyTr <- factor(paste(CO2$Type, CO2$Treatment, sep=""))
#（ただの文字列になってしまったので、要因として認識させる："factor"という関数）
TukeyHSD(aov(uptake ~ TyTr, data=CO2)) # 新変数でTukeyHSD法を実行する

############ 回帰分析 ######
reg1 <- lm(uptake ~ conc, data=CO2) # 上記のAnovaと中身が同じ構造！
summary(reg1)
# 切片や回帰係数（Estimateの部分）、それらがゼロと異なっているか否かの検定結果、
# R^2値、自由度調整R^2値などを出力

# Call:
# lm(formula = uptake ~ conc, data = CO2)
#
# Residuals:
# Min 1Q Median 3Q Max
# -22.831 -7.729 1.483 7.748 16.394
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 19.500290 1.853080 10.523 < 2e-16 ***
# conc 0.017731 0.003529 5.024 2.91e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 9.514 on 82 degrees of freedom
# Multiple R-squared: 0.2354, Adjusted R-squared: 0.2261
# F-statistic: 25.25 on 1 and 82 DF, p-value: 2.906e-06

library(car) # ANOVAと同じように、検定にかけてみる
Anova(reg1)
# 回帰係数の分散分析結果が得られる

# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# conc 2285 1 25.245 2.906e-06 ***
# Residuals 7422 82
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

テーマ:研究者の生活 - ジャンル:学問・文化・芸術

R合宿・第三部：グラフ作成編

R合宿の統括記事はこちら

# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）

# データの準備の仕方からグラフ作成、体裁作りまでを実行します

###### 基本的グラフ関数の紹介 #####################
# まずはシンプルなサンプルデータを使用
x0 <- rep(1:5, each=4) # 1~5を4個ずつ
y0 <- sort(rnorm(20, 10)) # 平均10の乱数を20個、ソートしておく

plot(y0) # 大概のデータは"plot()"に放り込めばとりあえず図示してくれる
boxplot(y0 ~ x0) # 箱ひげ図（~の意味は、y0をx0で分ける感じ → 統計でよく使う）
barplot(y0) # 棒グラフ
hist(y0) # ヒストグラム
plot(y0, type="o") # 折れ線グラフ
plot(x0, y0) # 散布図（先がX軸, 後がY軸）
plot(x0, y0, log="y") # y軸を対数スケール（log="x": x軸を対数、log="xy": 両対数）
plot(x0, y0, cex=c(1:20)) # 散布図でプロットサイズを変更（cex：サイズ指定）

###### Rグラフの取り扱い #################
# "ファイル" → "保存"で取り出せる（PDF形式）

###### 以下、現実的なデータセットと体裁、グラフ形態にしてみよう ###################
## サンプルデータを用意する。カテゴリー変数：２要因、連続変数：１要因
f1 <- rep(c("A", "B", "C","D"), each=10)
f2 <- rep(c("a", "b"), 20)
X <- sort(sample(0:100,40,replace=T)) # 0~100の数字をランダムに抜き出しソート（重複あり）
Y <- runif(40,0,100) # random uniformの意味。0~100の一様乱数を40個発生させる
D <- data.frame(X=X, Y=Y, F1=f1, F2=f2) # データフレーム化する

## 棒グラフ、折れ線グラフ用に平均と標準偏差を用意
MEAN1 <- tapply(D$Y, D$F1, mean) # 各グループの平均（"対象データ", "グループ分け", "関数名"）
SD1 <- tapply(D$Y, D$F1, sd) # 各グループの標準偏差
MEAN2 <- tapply(D$Y, list(D$F2, D$F1), mean) # 各グループの平均
SD2 <- tapply(D$Y, list(D$F2, D$F1), sd) # 各グループの標準偏差

###### グラフの体裁を整える準備 #####################
# par() に入れたパラメータの効果はウインドウを閉じるまで持続する
par(lwd=2, # line widthの意味。線の太さを標準の2倍にした
cex=1.5, # 文字サイズの指定。標準の1.5倍にした
mgp=c(1.75, 0.75, 0), # 軸タイトル、軸ラベル、軸線の間隔を指定
pin=c(5, 5), # 作図領域のサイズをインチ指定
family="Helvetica") # フォントをHelveticaに指定（推奨）。ほか Courier, Times など
# WindowsではHelveticaではなくArialにしてください（残念ながら用意されていないようなので）。
# （cf. ご自身で記述する時はこんなにマメに改行する必要はありません）

###### 箱ひげ図（１要因）の作成例 ######################
plot(Y ~ F1, data=D, # YをF1でグループ分け、データはDを用いる、の意味
lwd=2, ylim=c(0,100), # y-limitの意味。上限と下限を指定（計算式も利用可）
col="gray80", # カラーを指定：80%白の灰色。"colors()" で色一覧を見られる
axes=F) # いったん軸を除去（plotは、目盛りの線太だけは変えられないようなので）
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1)
# 軸の付け直し（二度手間だが、これをすることで目盛りも同じ線太にできる）

### さらに、図のラベルを左上部に追加する（論文に掲載する図を想定）
mtext("A) Coral growth", side=3, # side：3で上部を指定
at=0.8, line=0.5, # at：左からの横位置、line：グラフから上端からの縦位置
cex=1.5, font=2) # font=2：太字（普通でいいならfont項は不要）
# atの値はデータに応じ細かく調整要…正直、ドロー系ソフト上で調整する方が手っ取り早い

###### 棒グラフ（１要因の場合）の例 ########################
bar1g <- (barplot(MEAN1, # 図示したいデータ列
ylim=c(0,100), lwd=2, # 線太を2倍にした
las=1, # Y軸のラベルを水平方向に変更
col="wheat", # カラーを指定
xlab="group", ylab="volume")) # ラベル名を自分で書く
### エラーバー追加（barplot自体が各棒の中心のx座標データを持っているので利用する）
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2) # x1,y1,x2,y2：両端のx,y座標を指定
mtext("B) Coral happiness", side=3, at=1.2, line=0.8, cex=1.5, font=2)

# cf. points(x,y) # x,y座標を指定しプロット（cex=数字などの項を加え装飾可能）

###### 応用：間隔を変えた棒グラフ（control１群に対し、treatment複数群の場合）##############
barVSctrl <- (barplot(MEAN1,
space=c(1,1,0.5,0.5), # space：棒の幅1に対する棒の間隔を数値で指定
xlim=c(0,7), # xlim：x軸の範囲。棒の幅、棒の間隔の合計より少し大きくなるように設定
ylim=c(0,100), lwd=2, las=1,
col=c("white", "lightblue", "lightgreen", "mistyrose"),
xlab="group", ylab="amount"))
segments(barVSctrl, MEAN1+SD1, barVSctrl, MEAN1-SD1, lwd=2)
mtext("C) Fish sadness", side=3, at=1.0, line=0.8, cex=1.5, font=2)

####### 棒グラフ（２要因（２x２グループ)）の例 #############################
bar2g <- (barplot(MEAN2, beside=T, # TをFにすると積み重ねグラフになる
ylim=c(0,100), lwd=2, las=1, col=c("white", "salmon"),
xlab="group", ylab="Intensity"))
segments(bar2g, MEAN2+SD2, bar2g, MEAN2-SD2, lwd=2)
mtext("D) Shrimp madness", side=3, at=3.0, line=0.8, cex=1.5, font=2)

###### ヒストグラム ######################
hist(D$Y, breaks=c(0,20,40,60,80,100), # 区切りの指定（しなくてもOK）
main="", # でかでかタイトルを除去
lwd=2, las=1, col="seashell", xlab="velocity")
mtext("E) Mollusca marathon", side=3, at=25, line=0.8, cex=1.5, font=2)

###### 折れ線グラフ ######################
plot(MEAN1, type="o", ylim=c(0,100), xlab="factor", ylab="richness",
lwd=2, axes=F, cex=3, # プロット（文字サイズcex）を倍にする
pch=23, # プロットのシンボルを指定（21~25では輪郭と塗りつぶしの色が指定可能）
col="black", # 輪郭の色
bg="tomato") # 塗りつぶしの色
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
segments(c(1:4), MEAN1+SD1, c(1:4), MEAN1-SD1, lwd=2) # 今度はx軸は連続変数なので数値で指定
mtext("F) Polychaeta beauty", side=3, at=1.7, line=0.8, cex=1.5, font=2)

###### 散布図 + 回帰直線 ######################
plot(Y ~ X, xlim=c(0,100), ylim=c(0,120), xlab="Factor", ylab="Loneliness",
axes=F, cex=2, pch=21, col="black", # シンボルの形、輪郭の色は全系列で統一した
bg=c("seagreen","turquoise","slateblue","plum")[F1], data=D)
# 塗りつぶし色を系列によって変える（"F1"のA~D、4カテゴリーで色分けした）
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
legend("topright", legend=levels(D$F1), pt.cex=2, pch=21, # 凡例を追加
pt.bg=c("seagreen","turquoise","slateblue","plum"),horiz=T)
curve(0.2*x + 30, lwd=4, col="orangered", add=T) # 回帰直線：xの一次関数の例
mtext("G) Damselfish", side=3, at=5, line=0.8, cex=1.5, font=2)

############## 複数の図の重ね合わせ ######################
###### 散布図の複数系列を重ねる（上記とは別のやり方で描いてみる） ###############
# まず、枠を用意
matplot(c(0,100), c(0,100), xlab="Factor", ylab="Jealousy", axes=F,
type="n") # type="n"：グラフエリア内の余計な文字を消す
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
# １系列ずつ描く
matpoints(D$X, D$Y, cex=2, pch=21, col="black", bg="turquoise") # １系列目
matpoints(D$X*1.1, D$Y*1.1, cex=2, pch=24, col="black", bg="slateblue") # ２系列目
mtext("H) Sea urchin", side=3, at=5, line=0.8, cex=1.5, font=2)

############## 複数の図の並列 ######################
library(lattice) # latticeというパッケージを呼び出す
xyplot(Y ~ X | F1, cex=1.5, data=D) # ぴったり並列して書き出すので見やすい
# ただし、par()が効いていないようだし、plot()とは使い勝手が異なる模様…マニュアル複雑…

###### 異なるタイプのグラフを重ねる ################
# まず、一つ目のグラフ（棒グラフとする）
bar1g <- (barplot(MEAN1, xlim=c(0,5), # x範囲を二つ目のグラフと合致させる
ylim=c(0,100), lwd=2, las=1, col="khaki",
xlab="group", ylab="conflict"))
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2)
par(new=T) # 上書きするコマンド

# 二つ目のグラフ（折れ線グラフとする）
plot(bar1g, # エラーバーと同じ要領でbar1gをx座標に利用
MEAN1*1.5, xlim=c(0,5), # x範囲を一つ目のグラフと合致させる
type="l", # 線だけのグラフタイプを指定
axes=F, xlab="", ylab="", # 枠やラベルを描かないためのコード
lwd=4, col="royalblue")
axis(side=4, lwd=2, las=1) # 右サイドに軸だけを追加
mtext(expression(paste(factor2(m^2))), side=4, line=2, cex=1.5)
# 二つ目のグラフのラベルを右サイドに追加（expression, pasteを組合わせると上付き等も書ける）
mtext("I) Sea cucumber", side=3, at=0.1, line=0.8, cex=1.5, font=2)

####### 画面を分割して複数グラフをまとめる ########################
# parにmfrowを加える
par(mfrow=c(2,2), # multi frame by row の意味。この例では行方向に2x2=4分割
lwd=2, cex=1, mgp=c(1.75, 0.75, 0), pin=c(2.5, 2.5), family="Helvetica")
# 適当な4枚のグラフを実行してみよう
#（ウィンドウをドラグしてサイズを調整するとグラフの間隔を変えることができる）

####### 便利技：全変数を一括で図示してみる ########################
# 散布図形式で、全変数の関係を対戦表に配置する
# （自身vs.自身の位置に変数名が入ります。横方向は自身がY、縦方向は自身がXの場合）
pairs(D, panel=panel.smooth) # panel.smoothで変数間の関係を表す曲線が追加される
pairs(Y ~ X+F1+F2, panel=panel.smooth, data=D) # これでも同じ図になる

###### おまけ：EPSで出力する場合 ######################
# 下記、postscript…を実行してから、通常のグラフ作成コマンドを実行、
postscript("graphEPS.eps", height=4, width=4, # サイズを指定
horizontal=F, onefile=F, paper="special") # ここは暗号…
# ファイル名とheight、width（インチ単位…）をお好みで書き換えてください
# 作成終了したらdev.off()を実行（括弧内は空白）
dev.off() # 作業ディレクトリ内へ保存される

テーマ:研究者の生活 - ジャンル:学問・文化・芸術

R合宿・第二部：R基本操作編

R合宿の統括記事はこちら

# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）

###### R基本操作 ######################
# 今日はRらしく、コードで行ってみます。１行１行、数式を書いて実行して進めていくソフトです。コード打つのを怖がらないでくださいね…

# Excel: 広大な平地の駐車場に停められた車を扱うような感じ
# （配置を覚えておかなければ迷子になる。何をするにも毎回停めてある場所を探し出さなければ始まらない。方向音痴泣かせ…哀）

# R: 立体駐車場の入り口でIDを告げて預けられた車を呼び出すような感じ
#（こちらでIDさえ指定しておけば、配置や順番は問題にならない。ただし、現実の駐車場のような出し入れの煩わしさは無く、一瞬で預けられ呼び出せる！）

# コード打ちは試行錯誤が付きもの。なので、エディタ（ワープロソフト）に書き出し、保存しておくと利用しやすい。

###### 基本の"き" ##############################
# シャープ以降はその行は無視される（メモ書きに便利）
# 以下、これ通りに打って、出力を確認してみよう
# [1] というのは、単に１番目の出力結果であることを示す。そのうち、複数結果の例が出てきます。

2 + 1 # 足し算（+ の前後の半角スペースは付けなくてもいいが、付けた方が分かりやすい。ただし、全角スペースは厳禁）
2*5 # 掛け算
9/3 # 割り算
a <- 3 # 変数を定義する（< と - の組み合わせ。これの前後も半角スペースをいれるべし）
a # 定義した後でaと打つと3を返す、以下同様に、定義したものを出力してみよう
majimun <- 4 # 変数名は自由（数字から始めてはいけない、スペースやハイフンは使用不可）
b <- "A" # 文字は""で囲む
d <- "haisai" # こんなのもOK
d2 <- "夜露死苦" # 全角も行けないことはないが、半角英数だけにしておいた方が無難
d3 <- "(-_-#)/" # こら！そこ、脱線しないように！
e <- c(1,2,3,4,5) # cで囲むと数字列（ベクトルという）
f <- c(1:5) # 上と同じことを表す、コロンで結んでもよい
f + 1 # 全部に1を足す
f <- f + 1 # 矢印はイコールではない。同じ変数を使ってアップデート（再定義）できる。
f + c(5:1) # ベクトル同士の足し算
c(2:4) * c(3:5) # ベクトル同士の掛け算
rep(1,10) # 1を10個
rep(c(1:5), each=3) # 111222333444555
seq(1,10,by=2) # 1~10を2区切り

letters # 小文字のアルファベット一覧（下記の”データの操作”の要領で一部抜き出し可能）
LETTERS # 大文字のアルファベット一覧
month.abb # 月の3文字表記の一覧
month.name # 月のフル表記の一覧

paste("A",1,sep="") # 文字の合成（文字列作りに便利）
substr(“abcdef”, 2, 4) # 文字の抜き出し（この例では2~4文字目を抜き出し、bcdとなる）

matrix(0,5,5) # 行列（0を5行5列で）
matrix(c(1:16),4,4) # 1~16の4×4行列
v <- matrix(c(1:16),4,4,byrow=T) # 行方向に並べる場合（by row equal True）
colnames(v) <- c("sp1","sp2","sp3","sp4") # ラベルを付ける（cf. 行に名前を付けるのはrownames() ）
v # ラベルを付けたあとはこういう表示になる（実行例）

###### 基本的な関数 ##############################
# 平均 mean, 最大 max, 最小 min, 標準偏差 sd, 分散 var, データ個数 length
# 絶対値 abs, 指数 exp, 対数 log, 正弦&符号 sin, 小数点の丸め round, 小数点の切り捨て trunc
mean(f) # 平均値の例
se <- function(x) sqrt(var(x) / length(x)) # 標準誤差はfunction()を用いて自作関数で
# この場合、"x" 部分は他のアルファベットや単語でもOK。半角スペースは一応、空けておきましょう
log(f + 1) # 対数変換をするなら（0データ混じりの場合、1を足しておけばエラーにならない）
sqrt(f + 0.5) # べき乗変換（0データ混じりの場合、0.5を足しておけばエラーにならない）

# help!を呼び出したいとき
?関数名 # とやると、解説が呼び出される（ただし英語。必ずしも懇切丁寧ではないが…）
example(関数名) # とやると、実行例が見られる（必ずしも分かりやすくはないが…）

##### データの操作 ################################
f[2] # ベクトルの一部を取り出し（この場合は2番目を取り出し）
f[c(1,3)] # ベクトルの一部を取り出し（飛び飛びで、1,3番目を取り出し）
f[f>=3] # ベクトルの一部を取り出し（条件で）
f[-1] # ベクトルの一部を削除
v[2,] # 行列の行を取り出し。[行番号, 列番号]を意味する
v[,3] # 行列の列を取り出し
v[2,2] # 行列の 1 点を取り出し
v[1:2, ] # 複数行を取り出し
rbind(v[2,], v[2,3]) # 上下で連結する（cf. 左右で連結は cbind。r: row, c: column）

##### 作業場所（データ階層）の指定（作業ディレクトリの指定）#############################
# Rに読み込みたいデータファイルを置いた場所を作業ディレクトリに指定せよ
# Mac: メニューの"その他"→"作業ディレクトリの変更"。 command + d, command + d, returnで"デスクトップを指定できる。
# Win: メニューの"ファイル"→"ディレクトリの変更"
# cf. 作業ディレクトリの変更 http://cse.naro.affrc.go.jp/takezawa/r-tips/r/06.html
getwd() # 現在のディレクトリを確認する。この場合（）内は空白
setwd("でぃれくとり名") # ディレクトリの指定をする

##### データの準備・読み込み ##############################
# （予め用意したデータのファイルを読み込む場合）
# Excel などで csv（カンマ区切りファイル）保存、一行目を見出しとしたデータベース型を取ること

# データを "test.csv" というファイル名に保存し、そのデータを "test" としてRに読み込む場合のコード：
test <- read.csv("test.csv", header=T) # read.csv: データの読み込み; header=T: 一行目を見出しとして使用

# 今回は test.csv のサンプルデータを用意しておきました
X <- rep(c(1:5),each=40) # rep(何を,何回) 繰り返す
F1 <- rep(rep(c("A","B"),each=20),5) # 文字を繰り返す場合は""で挟む
F2 <- rep(rep(c("A1","A2","B1","B2"),each=10),5)
Y <- rpois(rep(1,200), # rpois(いくつ,平均値) のポアソン乱数（整数値）を発生
rep(seq(20,200,by=40),each=40)*rep(rep(c(1,1.4),each=20),5)*rep(rep(c(1.0,1.2,1.4,1.6),each=10),5))
test <- data.frame(X,F1,F2,Y) # ひとつのデータセットに統合

##### データの表記ルール ##################################
# 空白やスペース禁止（学名使用は注意）。ゼロなら0、欠損ならNAと入力する（NA：欠損、no data）
# 半角英数で記述すること（全角はダメ）。ピリオドはいいがハイフンはダメ、アンダーバーはOK。
# アルファベット順に順序づけられる
#（e.g. January, Februaryと入れると、Februaryの方が先 → 頭に数字を付けると簡単。1January, 2February）

##### 欠損値（NA）の扱い方 ##################################
# Rでは欠損値をNAと表記しておく。しかし、これを修正したり、ゼロに置き換える等の操作はデータフレームでは案外厄介…下記のようにやるとよいでしょう（cf. "!"は"それ以外"を意味します）。
test[!complete.cases(test),] # 欠損値のある行を表示（箇所を見つけて手動で修正する場合に使用）
test[,1][is.na(test[,1])] <- 0 # 欠損値を0に置き換える（この例はtestの1列目のNAを置き換える場合。列単位で作業する必要がある）
test2 <- na.omit(test) # 欠損値を含む行を除く（ベクトルならば、v[!is.na(v)]のようにしてNAを除ける）

##### パッケージのインストール #############################
install.packages("Rcmdr", dependencies=T) # パッケージをインストール（この場合はRコマンダーの例）、リストアップしてメモしておくとRを再インストールしたとき便利（定期的に最新にすべし）
update.packages("Rcmdr") # パッケージのアップデート

##### データの取り出し #######################
write.csv(test, "output.csv") # ファイルへ書き出し（データ構造により書き出せない場合あり）

# cf. Rの出力画面からコピペしたデータを加工するのに便利なサイト
# （単にコピペしただけだと、空白はスペースになってしまう））
# http://filter.webings.net/
# "連続半角スペース（1文字含む）→タブに変換" を選択、実行

fix(test) # （x：取り出したいものを入力）。コピー&ペーストで取り出せる（ベクトル、行列、データフレームなど）Win版では全体選択ができず利用できない？

# 最終手段はワープロアプリにコピペして、連続半角スペースを ,（カンマ）に変換し、拡張子をcsvに書き換えるとExcel上で復活させられる（多少ずれるでしょうけれど）

###### データの集計、ちょっとした計算など ##########################
test$F1 # testのF1だけを取り出し
test[test$F1=="A", ] # F1がAであるデータを取り出し
test[test$F1=="A" & test$F2=="A1", ] # F1がAで、F2がA1であるデータを取り出し

###### ある列を基準にデータの並べ替えをする ######################
# Excelには対応関数無し（メニューから実行するしかない）
test[order(test$F2, test$F1),] # F2, F1の優先順位で全体行を並べ替え
test[order(test$F2, decreasing=T),] # 降順で並べ替え

###### ヒストグラム #####################
# 頻度別に集計しなくていい（Excelではアドインを使った大がかりな作業…）
hist(test$Y)
hist(test$Y, breaks=c(0:5)*100) # 区切りを指定することもできる

###### 種の集計に便利（必須？）#####################
# （Excelではcountifを使った繁雑な作業が必要…）
spp1 <- rep(c("sp1","sp2","sp3","sp4"),c(5,4,3,0)) # サンプルデータ
spp10 <- rep(c("sp1","sp2","sp3","sp4"),c(0,2,3,1)) # サンプルデータ2
table(spp1) # 種の集計：これだけでも行けるが、いない種はインデックスが作られない
species <- factor(c("sp1","sp2","sp3","sp4")) # これを解消するにはインデックスを予め用意する
spp2 <- factor(spp1, levels=species) # インデックスと照らし合わせる
table(spp2) # いない種についてはゼロとなる

unique(spp1) # spp1の要素をリストアップ
union(spp1, spp10) # spp1、spp10に含まれる全種リスト
intersect(spp1, spp10) # spp1、spp10に共通な種のみのリスト

# 各種の多様度指数、多様度の解析
# （以下、community1, 2というデータを仮定し、計算してみる）
community1 <- sample(1:5, 10, replace=T) # 10種の個体数を想定。1~5個体を取るようにしてある
community2 <- t(sapply(1:5, function(s) sample(1:5, 10, replace=T))) # さらに、行列版

library(vegan) # 生物群集解析のためのパッケージ
diversity(community1, index="shannon", base=2) # Shannon-Wiener index
diversity(community1, index="simpson") # Simpson index (D)
diversity(community1, index="invsimpson") # Simpson index (1/D)

sim <- vegdist(community2, method="bray")
# 各行の群集間類似度を対戦表で表示（Bray-Curtis dissimilarity）
# 類似度（Bray-Curtis similarity）にするには全体を 1 から引く（1 - sim）

cruster <- hclust(sim, method="average") # 群平均法を用いたクラスター分析
plot(cruster, hang=-1) # custerを図示、hang=-1 で葉を揃える

# その他、MDS（多次元尺度構成法）なども"vegan"に含まれる
# cf. http://strata.uga.edu/software/pdf/mdsTutorial.pdf

##### 集計やまとめて統計量を算出するのに便利な技 #######################
apply(test[,c(1,4)],2,mean) # 1:各行に対し、2:各列に対し、関数を適用（この例ではmean）
tapply(test$Y,test$F1,mean) # F1の各要素について関数を適用（excelのsumif等の拡張版）
tapply(test$Y,list(test$F1,test$F2),mean) # F1, F2の各要素について（無い部分はNAになる）
# cf. 応用： tapply(test$Y,test$F1, function(x) sqrt(var(x) / length(x))) # SE を定義しつつ利用

sapply(1:5, exp)
sapply(1:5, function(x) x^2 + 2*x + 1)
# 各値を、x^2 + 2*x + 1 みたいに回帰式に当てはめるなど

##### 繰り返し作業の自動化＝プログラミング？の初歩 ######################
sapply(1:5, function(s) sample(1:5, 10, replace=T))
# 関数を定義しつつ、繰り返し計算の結果を行列に結合するなど、荒っぽい処理でも対応しやすい
# ただし帰ってくる値が、ちゃんと目的のものに合っているかどうか、
# 結果が明確なよりわかりやすい例でチェックすることを習慣づけるべきです

# もうひとつの繰り返し作業の例（for ループという）
# 他の言語でもよく使われるが、ちょっと処理速度は遅め
x <- 0 # xの元の値は0
for (i in 1:5) { x <- x + 1 }
# iを1~5まで繰り返し。左辺のxは新しく定義するx、右辺のは直前に定義したx
x # 繰り返し作業後のxを入力すると値が変わっているのが確認できる

# 7行おきに合計数を集計
xx <- matrix(1:56,28,2) # 集計をする元データ
xxx <- matrix(0,4,2) # 集計結果を納める入れ物
for (i in 1:4) {
xxx[i,1] <- sum(xx[((7*(i-1)+1):(7*i)),1] )
xxx[i,2] <- sum(xx[((7*(i-1)+1):(7*i)),2] )
}
xxx

##### RやRの機能を論文に引用する場合 #####################
citation() # 使用しているRを引用する場合のreference情報
citation(“パッケージ名”) # 特定のパッケージのreference情報

citation(“stats”) # 関数 glm を引用する場合
citation(“lme4”) # 関数 lmer を引用する場合

テーマ:研究者の生活 - ジャンル:学問・文化・芸術