理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

Rstudioの小ネタ (パッケージやファイルの保存方法)

Rstudioを閉じても、PC再起動してもファイルの読み込みやインストールしたパッケージは残せます (ただしパッケージは休んでいますので、起動するときにはlibraryで起こしましょう)

データ処理する前に必ず行う作業は以下の通りです
まずRstudioを起動させた状態で・・・
FileメニュからNewProject → NewProject → NewProjectと進んでいきます
Projectの名前を適当に書いてください(例:test001)、次に場所をしてしてください(as you like).

以下のように一つのフォルダの中に「---.Rhistory」と「---.Rproj」が入ってたら OK
データファイルも同じフォルダに収納します
f:id:yoshida931:20180604094238p:plain

次回から「---.Rproj」をダブルクリックしてRstudioを起動させてください
前回までに取り込んだデータは右上に表示され、そのまま使用可能です.
また右下に同じフォルダにあるデータファイルが表示されます.ここから直接データを取り込むこともできます.
またインストールしたパッケージも記憶されていますので、library( )で起動させるだけで使用できます.
f:id:yoshida931:20180604094641p:plain

Rで簡単 多重ロジスティック回帰分析

サンプルはRのmtcarsを使用します.車に詳しくないので回帰分析に相応しいサンプルか分かりませんが、学習のため数値のみ使用させていただきます.

サンプルから直接取り込み

dat <- mtcars[,c(8,2,5,6)]    #ここがポイント8列目を応答変数として、2列目・5列目・6列目を説明変数として設定
ans <- glm(dat$vs~.,data=dat,family=binomial) 
ans <- glm(dat$vs~dat$cyl+dat$drat+dat$wt,family=binomial)  

これで終了です.

summary(ans)

Call:
glm(formula = dat$vs ~ dat$cyl + dat$drat + dat$wt, family = binomial)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.43462  -0.15612  -0.06808   0.22595   1.35668  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   22.302     11.577   1.926   0.0541 .
dat$cyl       -3.375      1.544  -2.186   0.0288 *
dat$drat      -2.189      1.874  -1.168   0.2428  
dat$wt         1.777      1.660   1.070   0.2845  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.86  on 31  degrees of freedom
Residual deviance: 13.92  on 28  degrees of freedom
AIC: 21.92

Number of Fisher Scoring iterations: 7

オッズ比

exp(ans$coefficients)

 (Intercept)      dat$cyl     dat$drat       dat$wt 
4.849125e+09 3.423400e-02 1.120633e-01 5.911782e+00 

ロジスティック回帰分析(説明変数が単一かつ連続の場合)

投稿日:2018.2.13 最終更新日:2018.5.17

ロジット関数とロジスティック関数

準備として関数の特徴を押さえておきます

ロジット関数 \ \  y=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

(標準)ロジスティック関数 \ \ ロジット関数の逆関数=P = \frac{1}{1+exp(-y)}

mathwords.net

サンプルirisより

品種"virginica=1"、"別の品種=0"という2値を目的変数、”Sepal.Length(がく片の長さ)(cm)”を説明変数とします. ロジスティック回帰で求めることは、”Sepal.Length(がく片の長さ)(cm)”から予測される、品種が"virginica"になる確率です.

class(iris$Species)  #属性を確認
is <- as.numeric(iris$Species)   #  iris$Speciesが要因になっているので数字に変換(as.integerでも同じ)
y <- gsub(1, 0, is)  # is の1を0に変換
y <- gsub(2, 0, y)  # is の2を0に変換
y <- gsub(3, 1, y)  # is の3を1に変換

x <- iris$Sepal.Length

#実際のデータは以下のようになっています
> y
  [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [24] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [47] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [70] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [93] "0" "0" "0" "0" "0" "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[116] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[139] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
> x
  [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4 5.1 5.7 5.1 5.4 5.1 4.6
 [24] 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8
 [47] 5.1 4.6 5.3 5.0 7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8 6.2
 [70] 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4 6.0 6.7 6.3 5.6 5.5 5.5 6.1
 [93] 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8
[116] 6.4 6.5 7.7 7.7 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7 6.3 6.4
[139] 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9

ちなみに、そのまま単回帰分析を行ってみます.

summary(lm(formula=y~x) )
plot(x,y)
abline(lm(y~x)) #回帰直線
f:id:yoshida931:20180208184829p:plain:w300

目的変数は0、1なのですが、上図ではかなりはみ出してしまいます.0より少なくても、1より大きくても、良い予測とは言えません.したがって、この回帰直線は適当なモデルとは言えません.そこで指数を用いたロジスティックモデルを使用します.

ロジスティック変換 logistic transformation

目的変数が2値変数(品種=virginica=1、品種=別の品種=0)の回帰分析をおこなうためには、2値変数に与えられる上限と下限を取り除けば解決します.ロジスティック変換は、0~1の区間制限を取り除くことで、統計的推測に柔軟性を持たせた方法といえます.

品種が"virginica(y=1)"になる確率をPとした場合のロジスティック(ロジット)変換 を考えてみます.

オッズ=\frac{品種=virginica=1 である確率}{品種=他の品種=0 である確率}=\frac{p}{1-p}

ロジット(logit)とは、0から1の値をとる確率p に対し、そのオッズ\frac{p}{1-p}の対数から計算されます.

Pのロジスティック変換\ \  logit(p)=log\frac{P}{1-P}

f(P)=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

グラフで確認してみます

pl <- function(p){ return(log(p/(1-p))) }
plot(pl,xlab="P",ylab = "f(P)")
abline(0,0)
ロジットf(p)がどのような値になっても、0<p< 1となっています.
f:id:yoshida931:20180213133241p:plain:w300

\ \  f(P) = b0+b1*x

ロジスティック回帰モデル

2値の変数(virginica=1、別の品種=0)が1になる確率のロジットを目的変数にしたものが、ロジスティック回帰分析になります。

ロジスティック回帰モデル

\ \  f(P) = b0+b1*x

P = \frac{1}{1+exp(-(b0+b1*x))}

パラメータの推定(最尤推定法) P(yi=1|b0,b1)

尤度関数
L(θ)= \prod_{i=0}^{150} (\frac{1}{1+exp(-(b0+b1x))})^{yi})*(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}

このままでは大変なので対数をとります
対数尤度関数 L(θ)= \sum_{i=0}^{150} log(\frac{1}{1+exp(-(b0+b1x))})^{yi}log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi} 

yiは0か1です.例えば例題のy1と y150を見てみます.それぞれy1=1, y150=0となっています.
y1のときは
f\theta (y_1)=log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}
y150のときは
f\theta (y_{150})=log(\frac{1}{1+exp(-(b0+b1x))})^{yi}
となります.これを、1~150全て掛け合わせることになります.

ここからNewton法により最尤法を解いていくことになります.
私には解説が困難ので以下のサイトをご覧ください

https://qiita.com/NaoyaOura/items/6ad5142b0306476d9293qiita.com

Rを使ったロジスティック回帰分析

ここではRを使って楽させていただきます.一般化線形モデル(glm:Generalized Linear Models)は最尤法によりパラメータを推定する関数です! Rの関数はglmそのままです.

# glm関数で、応答変数の分布は二項分布なので引数をbinomialにすればOKのはず・・・
(ans <- glm(y~x,family=binomial))

Error in weights * y : non-numeric argument to binary operator
というエラーが出たので、例題で使用したxyの属性を調べます

typeof(x);typeof(y)
[1] "double"
[1] "character"

xは数で、yが文字になっていました.ということでyを数に置きなおして再トライ

y <- as.numeric(y)
(ans <- glm(y~x,family=binomial))

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

Coefficients:
(Intercept)            x  
    -16.320        2.592  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3   AIC: 121.3

#やっぱりRは便利です

b0=-16.32, b1=2.596という解になりました

#対数尤度を求めてみます
x1 <- x[101:150]
x2 <- x[1:100]
e1 <- exp(-(-16.3198+2.592*x1))
e2 <- exp(-(-16.3198+2.592*x2))
(ly <- sum(log(1/(1+e1)))+sum(log(1-1/(1+e2))))

[1] -58.67273

視覚的にイメージしておきます.

par(mfrow = c(1,2))
plot(x,y,yaxt="n",
     ylab = "目的変数(0,1)",
     xlab = "Sepal.Length(がく片の長さ)(cm)")     #まずは散布図
name<-c("0","1")
axis(side=2,at=c(0,1),labels=name) 
abline(lm(y~x))    #回帰直線
p <- function(x){ return(1/(1+exp(-(-16.320+2.592*x)))) } #pの関数を作ります
plot(p,-16.320+2.592*x,ylim=c(0,1),xlim = c(3.3,9.3),
     ylab = "virginicaである確率",
     xlab="b0+b1*x") 
f:id:yoshida931:20180517150938p:plain:w700

Rを使用して最尤法で推定したパラメータの有意性も確認することが可能です. 関数summaryを使用します.

summary(ans)

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9870  -0.5520  -0.2614   0.5832   2.7001  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -16.3198     2.6581  -6.140 8.27e-10 ***
x             2.5921     0.4316   6.006 1.90e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 190.95  on 149  degrees of freedom
Residual deviance: 117.35  on 148  degrees of freedom
AIC: 121.35     # -2*対数尤度+2*パラメータ数

Number of Fisher Scoring iterations: 5

標準誤差の算出については、またいつか・・・

オッズについて

ロジスティック変換より

がく片の長さxcmのオッズ=\frac{P}{1-P} = exp(-16.3198+2.5921x)

がく片の長さxcmから1cm増えた場合のオッズ =\frac{P}{1-P} = exp(-16.3198+2.5921(x+1))

オッズ比= \frac{がく片の長さがx+1cmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}{ がく片の長さがxcmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}= exp(-16.3198+2.5921(x+1))-exp(-16.3198+2.5921x)=exp(2.5921)=13.35646

xが1増えた場合の対数オッズ比は2.59、オッズ比は13.36になります.

 summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900 

xには大きな変動がないので、xが0.5増えた場合のオッズ比を求めてみます

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}=exp(2.5921×0.5)=3.654832

つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます.オッズ比にはxが含まれていないので、どの幅(xcm)のに対しても同じことが言えます.


参考・引用

回帰分析入門 (Rで学ぶ最新データ解析)

回帰分析入門 (Rで学ぶ最新データ解析)

  • 作者: 豊田秀樹,尾崎幸謙,川端一光,岩間徳兼,鈴川由美,久保沙織,池原一哉,阿部昌利,大橋洸太郎,秋山 隆,大久保治信
  • 出版社/メーカー: 東京図書
  • 発売日: 2012/01/13
  • メディア: 単行本(ソフトカバー)
  • 購入: 2人 クリック: 36回
  • この商品を含むブログ (6件) を見る

相関係数のイメージ

パッケージmvtnormを使用して相関係数0.0, 0.2, 0.5, 0.7, 0.8, 0.9のグラフを作成してみます

install.packages("mvtnorm")
library(mvtnorm)

共分散行列.分散を全て1に設定しているので共分散=相関係数となります.

sigma00 <- matrix(c(1,0,0,1), ncol=2)  
sigma02 <- matrix(c(1,0.2,0.2,1), ncol=2)
sigma05 <- matrix(c(1,0.5,0.5,1), ncol=2)
sigma07 <- matrix(c(1,0.7,0.7,1), ncol=2)
sigma08 <- matrix(c(1,0.8,0.8,1), ncol=2)
sigma09 <- matrix(c(1,0.9,0.9,1), ncol=2)

乱数から各相関図に300個プロットします

n00 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma00)  
n02 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma02)  
n05 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma05)  
n07 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma07)  
n08 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma08)  
n09 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma09)  

図を作成

par(mfrow = c(2,3))
plot(n00,xlab = "ρ=0.0",ylab="")  
plot(n02,xlab = "ρ=0.2",ylab="") 
plot(n05,xlab = "ρ=0.5",ylab="")    
plot(n07,xlab = "ρ=0.7",ylab="")
plot(n08,xlab = "ρ=0.8",ylab="")
plot(n09,xlab = "ρ=0.9",ylab="")

f:id:yoshida931:20180514142213p:plain:w800

変数の呼称について(目的変数と説明変数)

それぞれの研究界のご意見はあると思うのですが・・・
ややこしや

目的変数 Yは以下のように呼ばれています

目的変数 objective variable
応答変数 response variable 反応変数 reaction variable(response variable ) 結果変数 outcome variable
従属変数 dependent variable
基準変数 criterion variable
外的基準 external criterion
被説明変数 explained variable

説明変数 Xは以下のように呼ばれています

説明変数 explanatory variable
予測変数 predictor variable
独立変数 independent variable

パターンとしては

「目的変数&説明変数」
「従属変数&独立変数」
「被説明変数&説明変数」

個人的には…
応答変数&説明変数が理解しやすいかな?

2変量の正規分布をグラフでイメージ(persp)

また、ここで勉強させていただきました.
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html
忘れないように要点のみ転記させていただます.まさに備忘録.

今回はRの関数perspを使用して、密度関数の数式から3Dのグラフを描いてみます
確率変数x1、x2が正規分布に従い、無相関であることを仮定して進めていきます.

x1 <- seq(-3, 3, length=50)      # -3~3の範囲を50分割

head(x1)  #先頭部分だけ確認してみます
[1] -3.000000 -2.877551 -2.755102 -2.632653 -2.510204 -2.387755

x2 <- x1  

これで変数の設定は完了です.次に、分散1共分散0のマトリックスを作成します.

sigma.zero <- matrix(c(1,0,0,1), ncol=2)

     [,1] [,2]
[1,]    1    0
[2,]    0    1

Rの関数functionにx1、x2の同時確率密度を定義します. 上記sigma.zeroを分散共分散行列とする密度関数です  

f <- function(x1,x2) { dmvnorm(matrix(c(x1,x2), ncol=2), mean=c(0,0), sigma=sigma.zero) }

#6ポイントのみ確認してみます
f(-3,-3);f(0,0);f(3,3)
[1] 1.964128e-05
[1] 0.1591549
[1] 1.964128e-05

f(-1,-1);f(0,0);f(1,1)
[1] 0.05854983
[1] 0.1591549
[1] 0.05854983

関数outerを使用して、(x2,x2)であらわされる座標に対しFUN(f <- function(x1,x2))を適用します.
つまり、(x2,x2)に該当する値(f)が定まることになり3次元の描画を可能にします.

z <- outer(x1, x2, f)  
#50×50=2500個のデータが生成されました(length(z)で確認)
#10個だけ見てみましょう
z[1:10]
 [1] 1.964128e-05 2.814820e-05 3.973927e-05 5.526846e-05
 [5] 7.572219e-05 1.022015e-04 1.358875e-04 1.779878e-04
 [9] 2.296621e-04 2.919285e-04

http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html に書かれあります,
z の欠損値の置換は省略しております.
色塗りはRのヘルプ persp {graphics} をそのまま使用しました.

nrz <- nrow(z)
ncz <- ncol(z)
# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(x1, x2, z, col = color[facetcol], phi = 30, theta = -30)
f:id:yoshida931:20180512100244p:plain:w500

2変量の正規分布をグラフでイメージ(scatterplot3d)

ここで勉強させていただきました.
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html
忘れないように要点のみ転記させていただます.

必要なパッケージをインストールします

install.packages("mvtnorm")
library(mvtnorm)
install.packages("scatterplot3d")
library(scatterplot3d)

パッケージmvtnormは以下を参照
yoshida931.hatenablog.com

次に散布図を描きます

#共分散が0なので2変数には相関がありません
sigma.zero <- matrix(c(1,0,0,1), ncol=2)

     [,1] [,2]
[1,]    1    0
[2,]    0    1

#ともに平均=0,共分散=0となる2変数(乱数)を3組生成(n=100, n=1000, n=10000)
x100 <- rmvnorm(n=100, mean=c(0,0), sigma=sigma.zero)  
x1000 <- rmvnorm(n=1000, mean=c(0,0), sigma=sigma.zero)  
x10000 <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero)  

#3組並べて散布図
par(mfrow = c(1,3))
plot(x100)  
plot(x1000) 
plot(x10000)  

f:id:yoshida931:20180511164429p:plain:w500
確かに共分散はゼロです.

3Dで2変量の正規分布の同時確率密度関数を描いてみましょう

par(mfrow = c(1,3))
scatterplot3d(x100[,1], x100[,2], dmvnorm(x100, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)
scatterplot3d(x1000[,1], x1000[,2], dmvnorm(x1000, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)
scatterplot3d(x10000[,1], x10000[,2], dmvnorm(x10000, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)

f:id:yoshida931:20180511165033p:plain:w600