演習の概要

注意事項

  • 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

データ分析の準備

パッケージのインストール

  • 今回の演習では、以下のパッケージを使います。
install.packages("ggplot2")
install.packages("betareg")
install.packages("Zelig")
install.packages("lme4")
install.packages("gamlss.dist")
install.packages("gamlss")
install.packages("gamlss.util")
install.packages("brms")
library(ggplot2)
library(betareg)
library(Zelig)
library(lme4)
library(gamlss.dist)
library(gamlss)
library(gamlss.util)
library(brms)

モデル推定

ベータ回帰モデル

事例①

  • betaregパッケージのGasolineYieldデータを用いた分析を行います。
  • データの概要はこちらを参照してください。
  • データを読み込みます。
data("GasolineYield", package = "betareg")
  • 可視化します。ここでは、非説明変数をyeild(蒸留分別後に原油からガソリンに転換された割合)、説明変数をtemp(ガソリンが帰化する温度)としています。
ggplot(GasolineYield) + geom_point(aes(x = temp, y = yield, col = 2)) + 
  guides(colour = "none")

  • betareg()関数を用いて、ベータ回帰モデルを推定します。説明変数には、batchtempを用います。
gs.beta <- betareg(yield ~ batch + temp, data = GasolineYield)
summary(gs.beta)

Call:
betareg(formula = yield ~ batch + temp, data = GasolineYield)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-2.8750 -0.8149  0.1601  0.8384  2.0483 

Coefficients (mean model with logit link):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -6.1595710  0.1823247 -33.784  < 2e-16 ***
batch1       1.7277289  0.1012294  17.067  < 2e-16 ***
batch2       1.3225969  0.1179020  11.218  < 2e-16 ***
batch3       1.5723099  0.1161045  13.542  < 2e-16 ***
batch4       1.0597141  0.1023598  10.353  < 2e-16 ***
batch5       1.1337518  0.1035232  10.952  < 2e-16 ***
batch6       1.0401618  0.1060365   9.809  < 2e-16 ***
batch7       0.5436922  0.1091275   4.982 6.29e-07 ***
batch8       0.4959007  0.1089257   4.553 5.30e-06 ***
batch9       0.3857930  0.1185933   3.253  0.00114 ** 
temp         0.0109669  0.0004126  26.577  < 2e-16 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)    440.3      110.0   4.002 6.29e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood:  84.8 on 12 Df
Pseudo R-squared: 0.9617
Number of iterations: 51 (BFGS) + 3 (Fisher scoring) 
  • 非説明変数の予測精度を、計算します。
gs.beta.prd <- 
  cbind(
    Y=GasolineYield$yield,
    X=GasolineYield$temp,
    response=predict(gs.beta, type = "response"),
    link=predict(gs.beta, type = "link"),
    precision=predict(gs.beta, type = "precision"),
    variance=predict(gs.beta, type = "variance"),
    predict(gs.beta, type = "quantile", at = c(0.25, 0.5, 0.75))
  )
colnames(gs.beta.prd)<-c("Y","X","response", "link", "precision", "variance",
                      "q_0.25", "q_0.5", "q_0.75")
gs.beta.prd <- as.data.frame(gs.beta.prd)
  • 予測精度を可視化する。
ggplot(gs.beta.prd, aes(x = Y, y = response, col = 2)) + geom_point() + 
  geom_line(aes(x=Y, y=response), col="darkblue", lwd=1) +
  geom_ribbon(aes(ymin = q_0.25, ymax = q_0.75), col="darkblue", alpha = .2, lwd=0) +
  guides(colour = "none")

ガンマ回帰モデル

  • 分析には、Zeligパッケージのcoalitionデータを用います。
  • 参考文献1, 参考文献2を参考にしています。
data(coalition, package = "Zelig")
summary(coalition)
    duration         ciep12           invest           fract           polar      
 Min.   : 0.50   Min.   :0.0000   Min.   :0.0000   Min.   :349.0   Min.   : 0.00  
 1st Qu.: 6.00   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:677.0   1st Qu.: 3.00  
 Median :14.00   Median :1.0000   Median :0.0000   Median :719.0   Median :14.50  
 Mean   :18.44   Mean   :0.8631   Mean   :0.4522   Mean   :718.8   Mean   :15.29  
 3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:788.0   3rd Qu.:25.00  
 Max.   :59.00   Max.   :1.0000   Max.   :1.0000   Max.   :868.0   Max.   :43.00  
     numst2           crisis      
 Min.   :0.0000   Min.   :  0.00  
 1st Qu.:0.0000   1st Qu.:  0.00  
 Median :1.0000   Median : 10.00  
 Mean   :0.6306   Mean   : 22.38  
 3rd Qu.:1.0000   3rd Qu.: 29.75  
 Max.   :1.0000   Max.   :274.00  
  • 変数densityの分布を可視化します。
ggplot(coalition, aes(x=duration)) + 
  geom_density(color="darkblue", fill="lightblue") +
  geom_vline(aes(xintercept=mean(duration)),
               color="blue", linetype="dashed")

  • 変数fractdurationとの関係を可視化します。
ggplot(coalition) + geom_point(aes(x = fract, y = duration, col = 2)) + 
  guides(colour = "none")