線形モデルの予測関数:predictR()

線形モデルの予測値を導くためには、それなりに長いスクリプトを書く必要があり面倒です(参考)。また、lmerのようにランダム効果を入ってくると、信頼区間を推定するためのスクリプトが一気に煩雑になります。そこで、予測の平均値およびその信頼区間を計算する関数predictRを作りました。第一引数にモデルオブジェクトをいれ、第二引数に予測に用いたい説明変数Xの名前を指定するだけで、Xの最小値から最大値まで動かしたときの予測値と信頼区間を返します。なお、重回帰モデルの場合、第二引数以外の変数は平均(連続値)もしくは一番名前の若いグループ名(因子型)に固定されます。オフセット項についてもoffset=Tとすることでユニット当たりの予測値を返します。

SPONSOR LINK

関数の実装には、以下のスクリプトをRコンソールにコピペします。

predictR <- function(m, f_var, linkinv = T, ci.level=0.95, offset=F,
                     Interaction = F, int_var=NULL, pct = NULL,
                     n=100){
  ## Preparation for a new data set---
  # model class
  CL <- class(m)[1]
  
  # data frame
  dat <- model.frame(m)
    
    # Random effect
    if(CL=="lmerMod"|CL=="glmerMod"){
      REID <- NULL
      for(i in 1:length(names(ranef(m)))){
        tmp <- which(colnames(dat)==names(ranef(m))[i])
        REID <- c(REID, tmp)
      }
      dat <- dat[,-REID]
    }
    
    # offset term
    fixID <- which(colnames(dat)!="(offset)")
    ofsID <- which(colnames(dat)=="(offset)")
    dat <- dat[,fixID]
    if(CL=="lmerMod"|CL=="glmerMod"){
      dat.check <- length(names(dat)[-1]) - length(names(fixef(m))[-1])
    }else{
      dat.check <- length(names(dat)[-1]) - length(names(coefficients(m))[-1])
    }
    
    if(offset==T){
      if(length(ofsID)==0|dat.check>0){ stop("Specify offset term as e.g., 'offset = log(x)' in the model object") }
    }else{
      if(length(ofsID)!=0|dat.check>0){ stop("Should be 'offset=T' and specify as e.g., 'offset = log(x)' in the model object") }
    }
    
  # turn "character" or "logical" into "factor"
  dat <- as.data.frame(lapply(dat,function(x){
    if(is.character(x)|is.logical(x)){
      as.factor(x)
    }
    else{x}
  }))
  
  # choose focal variable
  fID <- which(names(dat)==f_var)
  if(is.factor(dat[,fID])){
    stop("focal variable must be numeric")
  }else{
    X1 <- list(seq(min(dat[,fID]),max(dat[,fID]), length = n))
    names(X1) <- f_var
  }
      
  # Interction terms
  if(Interaction == T){
    iID <- which(names(dat)==int_var) # column ID of interction terms
    if(length(iID)>1){ stop("Only one interaction term is allowed") }
    if(is.null(int_var)==T){ stop("Specify interaction terms") }
    
    if(is.factor(dat[,iID])==1){#if factor
        if(length(iID)==1){ X2 <- list(unique(dat[,iID])); names(X2) <- int_var }
      }else{#if numeric
        if(length(iID)==1){ X2 <- list(quantile(dat[,iID],pct)); names(X2) <- int_var }
        if(is.null(pct)==T){ stop("Specify percentiles of interaction term (ex. pct=c(0.2,0.8))") }
    }
    
    X <- c(X1,X2)
    X <- expand.grid(X)
  }else{
    if(length(int_var)>0){ stop("Should be 'Interaction = T'") }
    X <- X1
  }
  
  # make new data set
  if(Interaction==T){
    dat_r <- as.data.frame(dat[,-c(1,fID,iID),drop=F])
    if(dim(dat_r)[2]==0){
      new_dat <- X
    }else{
      fixed <- rep(lapply(dat_r,function(x) if(is.numeric(x)) mean(x) else factor(levels(x)[1],levels = levels(x))))
      new_dat <- cbind(X,as.data.frame(fixed))
    }
  }else{
    dat_r <- as.data.frame(dat[,-c(1,fID),drop=F])
    if(dim(dat_r)[2]==0){
      new_dat <- X
    }else{
      fixed <- rep(lapply(dat_r,function(x) if(is.numeric(x)) mean(x) else factor(levels(x)[1],levels = levels(x))))
      new_dat <- cbind(X,as.data.frame(fixed))
    }
  }
  
  ## Get predicted values---
  # interval
  cr <- ci.level + (1-ci.level)/2
  CR <- qnorm(cr,0,1)
  
  # inverse link function
  ilink <- ifelse(linkinv==T, family(m)$linkinv, identity)
  
  if(Interaction==T){
    output <- list(NULL)
    Index <- ifelse(is.factor(new_dat[,int_var]), length(unique(new_dat[,int_var])), length(pct))
    for(i in 1:Index){
      st <- I(1+(i-1)*n); en <- i*n
      if(CL == "lm"|CL == "glm"|CL == "negbin"){
        fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
        mm <- model.matrix(fo, new_dat[st:en,])
        us.fit <- mm%*%coefficients(m) 
        se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
        fit <- ilink(us.fit)
        uci <- ilink(us.fit+CR*se)
        lci <- ilink(us.fit-CR*se)
        y.tmp <- data.frame(fit=fit,uci=uci,lci=lci)
      }
      
      if(CL == "lmerMod"| CL == "glmerMod"){
        fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
        mm <- model.matrix(fo, new_dat[st:en,])
        us.fit <- mm%*%fixef(m) # unscaled fitted value
        se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
        fit <- ilink(us.fit)
        uci <- ilink(us.fit + CR*se)
        lci <- ilink(us.fit - CR*se)
        y.tmp <- data.frame(fit=fit,uci=uci,lci=lci)
      }
      output[[i]] <- data.frame(y.tmp, new_dat[st:en,])
    }#i
    u_int_var <- unique(new_dat[,int_var])# unique values or name for interaction terms
    if(is.factor(new_dat[,int_var])==T){
      names(output) <- as.character(u_int_var)
    }else{
      names(output) <- as.factor(paste(pct*100,"%","_",int_var,sep=""))
    }
    return(output)
  }else{
    if(CL == "lm"|CL == "glm"|CL == "negbin"){
      fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
      mm <- model.matrix(fo, new_dat)
      us.fit <- mm%*%coefficients(m) # unscaled fitted value
      se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
      fit <- ilink(us.fit)
      uci <- ilink(us.fit + CR*se)
      lci <- ilink(us.fit - CR*se)
      Y <- data.frame(fit=fit,uci=uci,lci=lci)
    }
    
    if(CL == "lmerMod"| CL == "glmerMod"){
      fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
      mm <- model.matrix(fo, new_dat)
      us.fit <- mm%*%fixef(m) # unscaled fitted value
      se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
      fit <- ilink(us.fit)
      uci <- ilink(us.fit + CR*se)
      lci <- ilink(us.fit - CR*se)
      Y <- data.frame(fit=fit,uci=uci,lci=lci)
    }
    output <- data.frame(Y,new_dat)
    return(output)
  }
}

引数は以下の通りです。
m: モデルオブジェクト(lm, glm, glm.nb, lmer)
f_var: 予測を行いたい説明変数の名前
linkinv: 逆リンク関数による応答変数の変換(デフォルトT)
ci.level: 信頼区間のレベル(デフォルト0.95)
offset: offset項の有無(デフォルトF)
Interaction: f_varとの交互作用項の有無(デフォルトF)
int_var: 交互作用項の名前(ひとつだけ指定可能)
pct: 交互作用項が連続値の場合、予測に用いる分位点(デフォルトNULL)
n: 予測に用いる説明変数の長さ(デフォルト100)

#lm, glm, glm.nb, lmer, glmer
#データを用意する
> D
  y1         x1 x2 x3
1  1 -0.2469068 19  C
2  3  0.7600840 18  T
3  3  0.6211508 16  C
4  3 -0.4965345 18  T
5  4  0.5435650 19  C
6  5  0.5004225 23  T

m <- glm(y1~x1+x2+x3, D, family=poisson)
fit <- predictR(m, f_var="x1")

> head(fit)
       fit      uci      lci        x1    x2 x3
1 2.351697 3.517449 1.185944 -1.610692 19.92  C
2 2.352722 3.501309 1.204135 -1.568597 19.92  C
3 2.353747 3.485274 1.222221 -1.526502 19.92  C
4 2.354772 3.469348 1.240196 -1.484408 19.92  C
5 2.355798 3.453539 1.258057 -1.442313 19.92  C
6 2.356823 3.437849 1.275796 -1.400218 19.92  C
#fit: 予測値の平均
#uci,lci: ci.levelで設定した信頼区間の上限・下限
#x1,x2,x3: 予測値を導く際の説明変数の代入値

#matplotで図示
matplot(fit[,"x1"], fit[,1:3], type="l", col="black")

#交互作用項(int_var)が因子の場合
m <- glm(y1~x1*x3, D, family=poisson)
fit <- predictR(m, f_var="x1", Interaction = T, int_var = "x3")

> head(fit$C,5)
       fit      uci       lci        x1 x3
1 1.295155 2.358462 0.7112373 -1.721629  C
2 1.318219 2.377811 0.7307986 -1.676699  C
3 1.341694 2.397378 0.7508794 -1.631769  C
4 1.365586 2.417169 0.7714919 -1.586839  C
5 1.389905 2.437190 0.7926487 -1.541909  C

> head(fit$T,5)
         fit      uci       lci        x1 x3
101 1.087035 1.653718 0.7145378 -1.721629  T
102 1.119484 1.692855 0.7403141 -1.676699  T
103 1.152902 1.732966 0.7669992 -1.631769  T
104 1.187318 1.774078 0.7946234 -1.586839  T
105 1.222761 1.816220 0.8232178 -1.541909  T

#交互作用項(int_var)が連続値の場合:int_varの分位点を指定
m <- glm(y1~x1*x2, D, family=poisson)
fit <- predictR(m, f_var="x1", Interaction = T, int_var = "x2", pct=c(0.2, 0.8)

> head(fit$`20%_x2`,5)
       fit      uci       lci        x1 x2
1 1.000534 1.678020 0.5965770 -1.721629 16
2 1.026491 1.707934 0.6169346 -1.676699 16
3 1.053121 1.738428 0.6379697 -1.631769 16
4 1.080443 1.769517 0.6597036 -1.586839 16
5 1.108473 1.801214 0.6821577 -1.541909 16
> head(fit$`80%_x2`,5)
         fit      uci       lci        x1 x2
101 1.248705 2.032481 0.7671726 -1.721629 23
102 1.277472 2.063563 0.7908334 -1.676699 23
103 1.306902 2.095166 0.8152064 -1.631769 23
104 1.337010 2.127301 0.8403113 -1.586839 23
105 1.367811 2.159981 0.8661687 -1.541909 23

※交互作用のないモデルでInteraction = Tとすると、int_varで指定した変数の相加効果のみが反映されます。
※式中にlog(x)のような関数が入るとエラーがでます。一度別名変数として格納することで実行可能です。

自作関数になりますので、誤り等がある可能性があります。各自の責任のもとでご利用ください。エラーなどありましたらコメントいただけると幸いです。

線形回帰の95%CI:ポリゴンver.

線形モデルの95%CIをグレイのシェードで描きたいときは、polygon()を使います。ポアソンモデルを例に、グレイシェイドの95%CIを描きます(下図)。

SPONSOR LINK

#1. 統計モデルの構築
model <- glm(y~x, family=poisson)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- exp(pred$fit) #予測値の平均値
up <- exp(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- exp(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogなので、predict()の返値をlogの逆関数exp()で変換する

#4. polygon()のx,yの値を設定する
X <- c(new$x, rev(new$x))
Y <- c(low, rev(up))
# rev()はベクターの並び順を逆転させる

#5. プロット
plot(y~x, type="n", ylim=c(min(Y,y),max(Y,y)),
     ann=F, axes=F)
box(bty="l")
polygon(X,Y,border = F, col=rgb(0,0,0,alpha=0.3))
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish abundance", side = 2, line = 3)

polygon()を後のほうに描いてしまうと、先に描いたポイントやラインも塗りつぶしてしまうので、先に描いておきます。また、プロットのy軸の範囲をylim=c(min(Y,y),max(Y,y))とすることで、データポイントおよび95%CIの最大・最小値をカバーできる描画領域を指定しています。

※polygon()補足
polygon()は、グラフ上のXY座標を指定し、その座標に囲まれた範囲を塗りつぶします。例えば、次のようなデータフレームがあるとします。

#ポリゴンの座標
X <- c(1,2,2,1); Y <- c(1,1,2,2)
data.frame(X,Y)
  X Y
1 1 1
2 2 1
3 2 2
4 1 2

これは、座標を(X, Y)と表記したとき、データ1が(1, 1)を示し、データ4が(1, 2)を示す、という意味なります。これをpolygon()に書かせてみます(右図)。polygon()では、第一引数でX座標を指定し、第二引数でY座標を指定します。

#ポリゴンを描画するプロット領域
plot(0, type="n", xlim=c(0,3),ylim=c(0,3),
     xlab="X", ylab="Y")
polygon(X,Y,col="black")
# colは塗りつぶしの色指定;borderは縁どりの有無を指定

線形回帰の95%CI:二項分布

二項分布は、「複数回の試行のうち、何回成功したか」のようなデータ(上限値が分かっている非負の離散値)に対して適用されることが多い確率分布です(ベルヌーイ分布の一般化)。生態学のデータでは、種の在/不在データに適用されることが多いと思います。二項分布のパラメーターは、ある事象Xが生じる確率pであり、pは0-1の範囲に収まる必要があります。そのため、リンク関数にはlogit関数(logit(p) = log(p/(1-p)))が用いられます。二項分布モデルの場合、predict()の返値をlogitスケールになっていますので、予測値をプロットするためには逆ロジット関数(inverse logit: ilogit(x) = exp(x)/(1+exp(x))で変換する必要があります。その他、注意点などについてはポアソン分布の場合と同様です。サンプルデータは記事末尾にあります。

#1. 統計モデルの構築
model <- glm(y~x, family=binomial)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

#4. inverse logit functionを作る
ilogit <- function(x){ exp(x)/(1+exp(x)) }
# Rのデフォルト関数にilogit関数はない

#5. predict()をつかって予測値を導く
y_pred <- ilogit(pred$fit) #予測値の平均値
up <- ilogit(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- ilogit(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogitなので、predict()の返値をlogitの逆関数ilogit()で変換する

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Probability of fish present", side = 2, line = 3)
ilogit <- function(x){ exp(x)/(1+exp(x)) }
N <- 50
a <- 0.4; b <- 1.5
x <- rnorm(N,0,1)

p <- y <- NULL
for(i in 1:N){
  p[i] <- ilogit(a + b*x[i])
  y[i] <- rbinom(1,1,p[i])
}

線形回帰の95%CI:ポアソン分布

ポアソン分布は、非負の離散値(カウントデータなど)に対して適用されることの多い確率分布です。ポアソン分布を仮定した線形モデルでは、母数である平均ラムダが負の値をとれないため、リンク関数にlogを使います。そのため、predict()の返値は対数スケールの値になっており、予測値曲線をプロットする際にはもとのスケールに戻す必要があります。記事末尾のサンプルデータを使い、ポアソンモデルの予測曲線をプロットします(右図)。

#1. 統計モデルの構築
model <- glm(y~x, family=poisson)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- exp(pred$fit) #予測値の平均値
up <- exp(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- exp(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogなので、predict()の返値をlogの逆関数exp()で変換する
# SEの演算は指数変換を施す前に行う
# ダメな例:exp(pred$fit) + exp(1.96*pred$se.fit)

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish abundance", side = 2, line = 3)

基本的なスクリプトは正規分布の場合と変わりませんが、辺値を変換するステップが必須となります。また、予測値の95%CIを出す際に、exp()で変換後に演算すると負の値が現れたりしますので注意が必要です。

※注 predict(model, new, se.fit=T, type=”response”)とすると、実測値スケールで予測値を返してくれますが、このSEをもとに四則演算を行うと上記の「ダメな例」と同様のことをすることになります。

N <- 50; a <- 1.4; b <- 0.7
x <- rnorm(N,0,1)

lambda <- y <- NULL
for(i in 1:N){
  lambda[i] <- exp(a + b*x[i])
  y[i] <- rpois(1,lambda[i])
}

線形回帰の95%CI:正規分布

lmやglmの予測値をプロットし、それらの信頼区間を合わせて表示したいケースは多いかと思います。正規分布のモデルであればpredict()の引数を調整することで簡単に描けますが、それ以外の確率分布となったときに案外描けなかったりします。以下のスクリプトでは正規分布を例に示していますが、この手順を覚えればほかの確率分布のモデルにも適用できます(多少の変更は必要)。プロットの基本的な調整でつくったサンプルデータを使ってスクリプトを書きます。以下のスクリプトから右図が描けます。

#1. 統計モデルの構築
model <- lm(y~x)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))
#予測用データフレーム内の説明変数の名称はモデルの変数名と一致する必要あり
#seq(min, max, length=サンプル数)となっており、最小値から最大値までの範囲で等間隔のベクターを生成

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- pred$fit #予測値の平均値
up <- pred$fit + 1.96*pred$se.fit #Upper 95%CI:平均予測値に1.96*SEを加える
low <- pred$fit - 1.96*pred$se.fit #Lower 95%CI

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

おまけ:matplot()を使った方法

Y <- cbind(up,y_pred,low)
matplot(new$x,Y, type="l", lty=c(2,1,2), col=rgb(0,0,0), ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

なお、複数の説明変数がある重回帰モデルの場合、注目する説明変数以外は平均値(もしくは中央値)に固定したデータセットを作る必要があります(そうしない限り、ほかの変数の影響を受けた予測値になる)。以下にサンプルを示します。

model <- lm(y~x1+x2+x3)
new_x1 <- seq(min(x1),max(x1),length=100)
new <- data.frame(x1 = new_x1, x2 = median(x2), x3 = median(x3))

半透明のシンボル

半透明のプロットを描く場合、シンボルの中身の色を調整できるシンボルを指定します。ここでは、シンボルの引数pchは21とし(円の中身の色指定ができるシンボル)、bgで中身の色を指定します。pch=21の場合、colはシンボルの縁どりの色を指定します。

rgb()はred,green,blueの三原色の割合を指定することで色調整する関数です。最初の三つの引数がそれぞれRGBに対応しています。すべて0とすると黒、1とすると白になります。alphaが透明度を指定する引数で、0(完全な透明)から1(完全な不透明)の値をとります。

#points()内の引数の値を調整
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x, y, cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3)
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

 

プロットの調整

下のplotを使い、見た目の調整方法についてまとめます。私個人の意見では、一度まっさらなプロット領域を用意し、好みの設定でデータや軸を追加していくのがのちのち微調整しやすいスクリプトになると思います。

調整前
#1. バックグラウンドとなるプロットを作る
plot(y~x, type="n", ann=F, axes=F)
#type="n"でプロット領域に何も書かない
#ann=F で両軸ラベルの消去
#axes=F で両軸の消去

#2. 外枠の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
#bty="l"でL字型の外枠

#3. データポイントの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
#cexでポイントサイズの調整:2とするとデフォルトの二倍の大きさ

#4. 軸の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex = 1.5)
axis(1); axis(2, las = 2)
#axisの第一引数は軸の位置を指定:1=下横軸、2=左縦軸、3=上横軸、4=右横軸
#las=2とすることで軸の値を横表記

#5. 軸ラベルの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
axis(1); axis(2, las = 2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)
#sideはラベル位置の指定:1=下、2=左、3=上、4=右
#lineはラベルの軸からの距離(行数)

※ステップ1~6まで継ぎ足しで書いてます

調整後

今回使ったサンプルデータの生成

N <- 50; a <- 3.4; b <- 0.2
x <- rnorm(N,0,1)
sigma <- 0.3

mu <- y <- NULL
for(i in 1:N){
  mu[i] <- a + b*x[i]
  y[i] <- rnorm(1,mu[i],sigma)
}

混合モデル(Mixing model)

地球上には様々な元素があり、その性質を決めるのは陽子数です。しかし、同じ元素でも(陽子数が同じ)、中性子数が異なるために性質は同じくしながら質量の異なる「安定同位体」が存在します。混合モデルは、安定同位体の比が捕食を通じて濃縮される特性を利用し、捕食者に対する餌生物の相対的貢献度を推定する手法です。炭素(C)や窒素(N)の安定同位体比がもっとも一般的で、それぞれの濃縮係数(栄養段階が一つ上がる際の安定同位体比の濃縮の程度)についての情報が充実しています。

古典的な混合モデル(IsoSourceなど)では平均値しか扱われてこなかったため、種内の安定同位体比のばらつきや濃縮係数の不確実性は無視されてきました。しかし最近では、これらの不確実性を明示的に取り入れたベイズ混合モデルが一般的になっています(SIARなど)。SIARに代表されるベイズ推定を利用した混合モデルでは、以下の点が古典的な混合モデルと異なります。

  • 種内(あるいは任意に設定した栄養ノード内)の安定同位体比のばらつきを考慮
  • 濃縮係数を確率変数とすることで、安定同位体比の濃縮プロセスの不確実性を考慮

SIARはRで実装されており、だれでも簡単に利用できるようになっています。実装の手順は以下の通りです(詳しい解説はSIARのマニュアルを参照)。

#1 パッケージインストール
install.packages("siar")
library(siar)

#2 消費者のデータ
>consumer #should be matrix
      d15N      d13C
1 2.765054 -24.72063
2 9.857108 -24.22871
3 3.562164 -25.20802
4 7.623977 -23.86547
5 2.181553 -22.99689

#3 エサ生物のデータ
>prey
    Source mean_d15N  sd_d15N mean_d13C  sd_d13C
1 detritus      2.23 2.565637    -25.22 0.932000
2 plankton      0.13 3.711336    -22.33 3.461702
3   shrimp      4.43 3.930261    -23.34 2.922652

#4 濃縮係数のデータ
>tef
    Source mean_d15N sd_d15N mean_d13C sd_d13C
1 detritus      3.54    0.74      1.63    0.63
2 plankton      3.54    0.74      1.63    0.63
3   shrimp      3.54    0.74      1.63    0.63

#5 混合モデルの実行
result <- siarmcmcdirichletv4(consumer,prey,tef,concdep=0, iterations = 500000, burnin = 50000)

#6 推定結果の表示
Summary information for the output file ...
         Low 95% hdr High 95% hdr       mode      mean
detritus   0.1470012    0.9550798 0.43382742 0.5338503
plankton   0.0000000    0.5090717 0.03679839 0.2203550
shrimp     0.0000000    0.5362734 0.04320550 0.2457947
SD1        0.0000000    7.2432659 2.25476445 3.1661398
SD2        0.0000000    4.6901279 1.14716487 1.8523370

おまけ:グラフの表示

 siarproportionbygroupplot(model1) 

階層分割(Hierarchical variation partitioning)

階層分割(Hierarchical variation partitioning)は、線形モデルをもとにした分散分割の手法です。ある説明変数セットのもとで、すべての組み合わせのモデルを構築し、それらのフィッティング(決定係数など)を比較することで各説明変数の相対的な説明力を導きます。多重共線性の影響を受けにくいことから、野外データから各説明変数の相対的重要性を示すのに適した手法といえます。ただし、あまりにも説明変数が多い場合、モデル間の適合度指標の差を計算するプロセスが複雑化するため(ただ引いてるだけですが)、あまり信頼のおける結果は得られないように思います(説明変数はせいぜい3~5程度が限界?)。

階層分割はパッケージhier.partを用いることで容易に実装できます。

  • hier.part()
    パッケージ:hier.part
    family指定: 正規分布、ポアソン分布、二項分布
    ランダム効果:不可
    適合度指標:対数尤度(logLik)、決定係数(Rsqu)、Root-mean-square ‘prediction’ error(RMSPE)
    備考:Deviance explainedなどの適合度指標も計算は可能。
#sample
library(hier.part)
hier.part(y, X, family="gaussian", gof = "Rsqu")
# y: response variable
# X: dataframe of explanatory variables

# 結果の表示---------------------------------------------
# GOF for each candidate model
$gfs
[1] 0.00000000 0.31519073 0.01143161 0.32120974

# Independent and joint contributions of each variables
$IJ
             I         J      Total
x1 0.312484426 0.0027063 0.31519073
x2 0.008725311 0.0027063 0.01143161

# Percentage of independent contributions
# Note: the denominator is "total variance explained"
$I.perc
           I
x1 97.283609
x2  2.716391

線形モデル

glm()

パッケージ:デフォルト
family指定:正規分布、ポアソン分布、二項分布
ランダム効果:可
モデル選択:dredge, stepAIC
備考:特になし

# sample
glm(y ~ x, data, family = poisson)

glmer()

パッケージ:lme4
family指定:ポアソン分布、二項分布(正規分布の場合はlmer();R ver 3.Xあたりから仕様が変更?)
ランダム効果:可
モデル選択:dredgestepAIC
備考:
ランダム効果が多くなると推定精度が低下することが指摘されている。Bolkerさんの論文によれば、ランダム効果が3つ以上の場合はMCMCなど別の手法を採用したほうがいいとしている。ただし、これまで実装されていたmcmcsamp()という関数に関しては、分散の推定値がゼロ付近のときに適切な推定を行ってくれない可能性があるらしく、2017年5月現在はパッケージから関数が除外されている模様。なお、推定が収束しなかった場合、control=glmerControl(optCtrl=list(maxfun=20000))などと付け加えることで、試行数の上限を増やすことができる。

# sample 1: varying intercept model
glmer(y ~ x + (1|group), data, family = poisson)

# sample 2: varying intercept and slope model
glmer(y ~ x1 + (1 + x1|group), data, family = poisson)

# sample 3: nested random effect model
glmer(y ~ x1 + (1|region/group), data, family = poisson)
# region is a larger set of groups

# sample 3: nested random effect model
lmer(y ~ x1 + (1|region), data, control=lmerControl(optCtrl=list(maxfun=20000)))
glmer(y ~ x1 + (1|region), data, family=poisson, control=glmerControl(optCtrl=list(maxfun=20000)))
# when did not reach convergence with default settings

glm.nb()

パッケージ:MASS
family:負の二項分布に固定
ランダム効果:不可
モデル選択:dredge, stepAIC
備考:特になし

# sample
glm.nb(y ~ x, data)

glmmadmb()

パッケージ:glmmADMB
family:ポアソン分布、負の二項分布、ベータ分布、ベータ二項分布
ランダム効果:可
モデル選択:不可(関数を自作すれば可)
備考:zeroInflation=TRUEとすることでゼロ過剰モデルも実装可

# sample 1: poisson model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "poisson")

# sample 2: negative binomial model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "nbinom")

# sample 3: beta model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "beta")

# sample 4: zero inflated poisson model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "poisson", zeroInflation=TRUE)