濃淡のある信用区間のプロット

WinBUGSやJAGSを用いたベイズ推定を使うと、各パラメーターに対する大量のMCMCサンプルが得られます。このMCMCサンプルを使って95%信用区間を計算し、それをもとに議論しますが、こうした任意の信用区間ではなく、事後確率がどのような分布をしているのかのほうを見せたい場合もあります。そんな時に便利なパッケージdenstripを紹介します。

SPONSOR LINK

濃淡のあるパラメーター推定結果の図示

ここでは、行列XにパラメーターのMCMCサンプルの入っているとしてみましょう。各列が各パラメーター、各行がMCMCサンプルに対応します。

X <- matrix(NA,nrow=1000,ncol=3)
for(i in 1:ncol(X)){
  X[,i] <- rnorm(1000,rnorm(1,2,2),runif(1,0.1,1))  
}
colnames(X) <- c("p1","p2","p3")
> head(X)
           p1        p2       p3
[1,] 2.527981 0.9171227 4.432790
[2,] 3.007233 2.9618068 4.689230
[3,] 1.891604 2.8134809 4.616570
[4,] 2.741595 0.8363671 4.741110
[5,] 2.602177 3.3400819 4.502327
[6,] 1.603801 2.5635456 4.356336

このデータを使って、事後分布と比例するグラデーション付きのパラメーター推定値を図示してみます。手順としては、まず空のプロットを呼び出し、そこにdenstrip関数をつかってストライプを図示していくことになります。

library(denstrip) #load package
plot(0,type="n",xlim=range(X),ylim=c(1,3)) # make plot region
for(i in 1:3) denstrip(X[,i],at=i,mticks=median(X[,i]),mcol="white")
# at = i is y-axis values density stribe to be plotted
# mticks is the x-position for median or similar
# mcol is the tick color to be used

出来上がりはこんな感じです。濃淡は、特に指定しなければdensity関数から求められる確率密度カーネルを基に図示されます。

また、以下のようにすれば見栄えも調節できます。

plot(0,type="n",xlim=range(X),ylim=c(1,3), ann=F, axes=F)
for(i in 1:3){
  denstrip(X[,i],at=i,mticks=quantile(X[,i],c(0.025,0.5,0.975)),
            mcol=c("black","white","black"))
}#adding 95% credible intervals
box(bty="l")
axis(1); axis(2,at=c(1,2,3),labels=colnames(X),las=2)

濃淡のある線形回帰の信用区間

パラメーターの推定結果に加え、線形回帰の推定結果についても濃淡をつけて図示することが可能です。基本的な手順は一緒ですが、denstrip関数の代わりにdensregion関数を使います(densregion関数についての詳細はこちらを参照)。まず、説明変数xについて、線形モデルによる予測値のMCMCサンプルが入った行列を用意します。先ほどと同じように、列が予測値と対応しており、行がMCMCサンプルの繰り返しになります。

> head(MCMC[,1:3])
     ypred[1] ypred[2] ypred[3]…
[1,] 0.336092 0.333215 0.330362
[2,] 0.384028 0.379466 0.374958
[3,] 0.469216 0.467234 0.465260
[4,] 0.259841 0.256886 0.253964
[5,] 0.350182 0.347407 0.344654
[6,] 0.597118 0.593395 0.589696
…

このデータをもとに、下記のように整形していきます。

est <- apply(MCMC,2,median) # median estimate
x <- seq(min(P),max(P),length=100) # range of explanatory variable
y <- seq(0,1.6,length=512) # range of y-axis; probability density will be estimated in this region
z <- matrix(nrow=length(x),ncol=length(y)) # empty case for density values
for(i in 1:length(x)) z[i,] <- density(MCMC[,i],from=min(y),to=max(y))$y # insert density values

library(denstrip)
plot(0,type="n",ylim=range(y),xlim=range(x), ann=F, axes=F)
densregion(x,y,z, colmax="red", gamma=0.8,scale=0.75,nlevels=75)
lines(x,est,lwd=2,col=grey(0,0.65))
axis(1);axis(2,las=2);box(bty="l")

こんな感じになります。

行列操作

Rの行列操作になれる

Rの行列指定になれると、いろいろと作業が簡略化できます。備忘録もかねて、整理したいと思います。

SPONSOR LINK

数字による行列指定

Rでは「行列[行番号,列番号]」の形で行列が指定されます。例として、3×3の行列を作り、特定のデータを行列の番号で指定してみます(Xは行列ではなく、データフレームでも可)。

> X <- matrix(rnorm(9,0,1),3,3)
> X
          [,1]       [,2]       [,3]
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[2,3]
[1] -0.3638491

複数のデータを同時に指定する場合は、”1:3″(1から3まで)や”c(1,3)”(1と3を指定)などを使います。また、空欄”[1,]”とすると、その空欄の部分(行もしくは列)のすべての要素が指定されます。

> X[c(1,2),1]
[1] 2.3587491 0.8162371

> X[1:3,1]
[1] 2.3587491 0.8162371 1.8183986

> X[1,c(1,2)]
[1]  2.3587491 -0.8524698

> X[1,]
[1]  2.3587491 -0.8524698  1.8697275

特定の行や列を取り除きたいときは、数字を負にします。

> X[,-2]
          [,1]       [,2]
[1,] 2.3587491  1.8697275
[2,] 0.8162371 -0.3638491
[3,] 1.8183986 -0.7109136

> X[,c(-2,-3)]
[1] 2.3587491 0.8162371 1.8183986

要素名による行列指定

要素名を指定することも可能です。

> colnames(X)<-c("a","b","c")
> X
             a          b          c
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[,"a"]
[1] 2.3587491 0.8162371 1.8183986

複数ある行のうち、特定の列名(行名)に一致する列番号を取り出したい場合は、colnames(rownames)を使います。複数の列名(行名)に一致する番号を知りたい場合は、”==”ではなく”%in%”を使います。

> colnames(X)
[1] "a" "b" "c"

> which(colnames(X)=="a")
[1] 1

> which(colnames(X) %in% c("a","b"))
[1] 1 2

> X[,which(colnames(X) %in% c("a","b"))]
             a          b
[1,] 2.3587491 -0.8524698
[2,] 0.8162371  0.4589377
[3,] 1.8183986  1.2378938

> X[,-which(colnames(X) %in% c("a","b"))]
[1]  1.8697275 -0.3638491 -0.7109136

複数ある行のうち、特定の列名(行名)に一致する列番号だけをのぞきたい場合は、”!=”も使えます。ただし、これは指定する行もしくは列が一つだけの場合に限られるので、「-which(colnames(X) %in% c(“a”,”b”))」などとするほうが汎用性はあるように思います。

function関数の使い方

function関数

Rで複数の関数を組み合わせて値を算出しなければならないとき、function関数で定義してしまうと楽です。しかし、いったい何をやってるのかわかりにくいところもあります。ここでは、function関数を使った関数の定義について簡単に紹介します。

変動係数

変動係数(CV)とは、標準偏差を平均で除したもので、変動性の比較などでよく使われる指標です。RのデフォルトでCVは定義されていないのですが、比較的よく使うものなので、例として挙げてみます。まず、function関数の仕組みですが、引数の指定と、その引数を使った計算式の指定、という形になっています。

CV <- function(x){ sd(x)/mean(x)}

#example
> x <- rnorm(100,10,10)
> CV(x)
[1] 1.042615

上記のfunction()のかっこ内に引数(今回はx)として何を使うのか指定し、{}内に引数を使った計算式を定義します。

ベクター中の0要素をカウントする関数

少し込み入ったものとして、ベクター中に0がいくつかあるのかカウントし、そして0が何番目の要素かを返す関数を定義してみます。

zerocount <- function(x){
 n_zero <- length(which(x==0))
 id_zero <- which(x==0)
 return(list(n_zero=n_zero, id_zero=id_zero))
}

#example
>x <- rpois(100, 2)
> zerocount(x)
$n_zero
[1] 17

$id_zero
 [1] 12 16 21 22 25 31 35 40 43 67 74 77 78 79 86 87 95
 

こちらでは、返す要素が二つあるので、返り値をリスト化して返すよう指定しています(4行目)。

応用

発展的には、乱数発生などのプロセスを組み合わせることで、簡単なシミュレーションモデルを作ることができます。そのほか、自作の作図関数なども作ることができます。興味のある人はぜひトライしてみてください。

最頻値の推定関数

最頻値

平均、中央値、分散など、確率変数を特徴づける指標はいくつもあるけれども、多くのものについてはRで標準実装されている。しかし、地味に知りたい「最頻値」なるものは実装されていない。

最頻値の推定の仕方

離散値の場合は以下のコードで簡単に推定できる。

x <- rpois(100,10)
names(which.max(table(x)))

しかし、サンプルサイズが比較的小さい場合や、確率変数が連続値の場合にはとてもいい方法とは思えない。そこで、table関数の代わりに、density関数をつかった関数を定義する。

mfv <- function(x,from=min(x,na.rm=T),to=max(x,na.rm = T)){
  if(is.numeric(x)==F)stop("use numeric")
  tmp <- density(x,from=from,to=to)
  Mode <- tmp$x[which.max(tmp$y)]
  return(Mode)
}

# example
> x <- rnorm(1000,10,20)
> mfv(x)
[1] 12.16744

density関数は、確率変数の確率密度(のようなもの;実際はカーネル)を変数の傾度にそって推定してくれるもの。上で定義したものは、このカーネルの値が最も高くなる確率変数xの値を返す、という関数になる。細かいけど、いちいち計算するのは面倒なので、関数化してしまうのが楽。。。ちなみに、fromとtoに値を指定すると、その範囲の中だけでカーネルを推定する。例えば、カウントのような正の値しかとらない変数の最頻値を知りたいときは、from=0とすることで、正の実数の範囲の中だけで最頻値を計算してくれる。また、確率など0-1しかとらない、とわかっているのであれば、from=0、to=1とする。

だが、そもそもdensity関数がなにをどうやって推定しているのかようわかっていないので、正確ではない可能性がある。確認次第追記します。

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

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

関数の実装には、以下のスクリプトを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を描きます(下図)。

#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)