統計解析R演習8

注意事項

• 慶應義塾大学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()関数を用いて、ベータ回帰モデルを推定します。説明変数には、batchとtempを用います。

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

• 変数fractとdurationとの関係を可視化します。

ggplot(coalition) + geom_point(aes(x = fract, y = duration, col = 2)) + 
  guides(colour = "none")

• ガンマ回帰モデルを推定します。

coa.gam <- glm(duration ~ fract + numst2, data=coalition,
                family=Gamma)
summary(coa.gam)

Call:
glm(formula = duration ~ fract + numst2, family = Gamma, data = coalition)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2510  -0.9112  -0.2278   0.4132   1.5360  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.296e-02  1.329e-02  -0.975  0.33016    
fract        1.149e-04  1.723e-05   6.668 1.19e-10 ***
numst2      -1.739e-02  5.881e-03  -2.957  0.00335 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.6291004)

    Null deviance: 300.74  on 313  degrees of freedom
Residual deviance: 272.19  on 311  degrees of freedom
AIC: 2428.1

Number of Fisher Scoring iterations: 6

• 回帰係数の指数をとった値を示します。

exp(coef(coa.gam))

(Intercept)       fract      numst2 
  0.9871239   1.0001149   0.9827628 

• zellig()関数を用いてモデル推定を行うこともできます。

z.out <- zelig(duration ~ fract + numst2, model = "gamma", data = coalition)

How to cite this model in Zelig:
  R Core Team. 2007.
  gamma: Gamma Regression for Continuous, Positive Dependent Variables
  in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
  "Zelig: Everyone's Statistical Software," http://zeligproject.org/

x.low <- setx(z.out, numst2 = 0)
x.high <- setx(z.out, numst2 = 1)
s.out <- sim(z.out, x = x.low, x1 = x.high)
summary(s.out)

 sim x :
 -----
ev
         mean       sd      50%     2.5%    97.5%
[1,] 14.42859 1.059792 14.38851 12.61026 16.58989
pv
         mean       sd      50%      2.5%    97.5%
[1,] 14.64059 12.75808 10.76264 0.7475632 46.15597

 sim x1 :
 -----
ev
         mean       sd     50%     2.5%    97.5%
[1,] 19.25073 1.109373 19.1745 17.29701 21.62251
pv
         mean       sd      50%     2.5%    97.5%
[1,] 19.24491 16.58495 15.00363 1.028835 59.09073
fd
         mean       sd      50%     2.5%    97.5%
[1,] 4.822135 1.518249 4.851093 1.766683 7.658619

• 推定結果を可視化します。

plot(s.out)

指数-ガウス回帰モデル

• まず、指数-ガウス分布データを作成してみましょう。MrUnadonさんのコードを使わせていただいています。
• rnorm()関数で正規分布（ガウス分布）の乱数を、rexp()関数で指数分布の乱数を生成します。

df_norm<-data.frame(norm = rnorm(1000, mean = 10, sd =1), exp = rexp(1000,rate = 0.5))

• 次に、gamlss.distパッケージのrexGAUS()関数を使って、指数-ガウス分布の乱数を生成します。

df_exGauss<-data.frame(exGauss = rexGAUS(2000,mu=10,sigma=1,nu=1/0.5)) #nu=1/rate

• ２つのデータフレームを結合します。
• パイプ%>%を有効化するには、needs::prioritize(magrittr)と入力してください。

df_bind<-cbind(df_norm,df_exGauss)%>%
  dplyr::mutate(norm_exp = norm + exp)%>%
  tidyr::gather()

• 作成した乱数データを可視化します。

df_bind2<-transform(df_bind,key = factor(key,levels = c("norm","exp","norm_exp","exGauss")))
ggplot(df_bind2,aes(x=value,colour=key,fill=key))+#ggplot2で描画
  geom_histogram(alpha=0.5)+
  geom_vline(xintercept = 0,linetype=2)+
  theme_bw()+facet_wrap(~key)

• 指数-ガウス回帰モデルの推定例として、lme4パッケージのsleepstudyデータを用います。データの概要はここを参考にしてください。

data(sleepstudy)
summary(sleepstudy)

    Reaction          Days        Subject   
 Min.   :194.3   Min.   :0.0   308    : 10  
 1st Qu.:255.4   1st Qu.:2.0   309    : 10  
 Median :288.7   Median :4.5   310    : 10  
 Mean   :298.5   Mean   :4.5   330    : 10  
 3rd Qu.:336.8   3rd Qu.:7.0   331    : 10  
 Max.   :466.4   Max.   :9.0   332    : 10  
                               (Other):120  

• 変数DaysとReactionとの関係を可視化します。

base.plot <- ggplot(sleepstudy) + geom_point(aes(x = Days, y = Reaction, col = 2)) + 
  guides(colour = "none")
base.plot

• 変数Reactionのヒストグラムを作成します。

hist.plot <- ggplot(sleepstudy, aes(x = Reaction, fill=2))+
  geom_histogram(alpha=0.5)+
  geom_vline(xintercept = 0,linetype=2, col="darkblue")+ 
  guides(colour = "none")
hist.plot

• ここでは、①glmssパッケージのglmss()関数を用いる方法と、②brmsパッケージを用いてベイズ推定する方法の２つを紹介します。
• ①glmssパッケージのglmss()関数を用いる方法

sleep.model<-gamlss(Reaction~Days,
            data=sleepstudy,
            family=exGAUS, 
            mu.fix=FALSE,
            sigma.fix = FALSE,
            nu.fix = FALSE)

GAMLSS-RS iteration 1: Global Deviance = 1900.86 
GAMLSS-RS iteration 2: Global Deviance = 1900.163 
GAMLSS-RS iteration 3: Global Deviance = 1900.13 
GAMLSS-RS iteration 4: Global Deviance = 1900.105 
GAMLSS-RS iteration 5: Global Deviance = 1900.086 
GAMLSS-RS iteration 6: Global Deviance = 1900.075 
GAMLSS-RS iteration 7: Global Deviance = 1900.067 
GAMLSS-RS iteration 8: Global Deviance = 1900.061 
GAMLSS-RS iteration 9: Global Deviance = 1900.055 
GAMLSS-RS iteration 10: Global Deviance = 1900.051 
GAMLSS-RS iteration 11: Global Deviance = 1900.046 
GAMLSS-RS iteration 12: Global Deviance = 1900.043 
GAMLSS-RS iteration 13: Global Deviance = 1900.04 
GAMLSS-RS iteration 14: Global Deviance = 1900.037 
GAMLSS-RS iteration 15: Global Deviance = 1900.034 
GAMLSS-RS iteration 16: Global Deviance = 1900.032 
GAMLSS-RS iteration 17: Global Deviance = 1900.029 
GAMLSS-RS iteration 18: Global Deviance = 1900.026 
GAMLSS-RS iteration 19: Global Deviance = 1900.024 
GAMLSS-RS iteration 20: Global Deviance = 1900.023 

Algorithm RS has not yet converged

summary(sleep.model)

******************************************************************
Family:  c("exGAUS", "ex-Gaussian") 

Call:  gamlss(formula = Reaction ~ Days, family = exGAUS,  
    data = sleepstudy, mu.fix = FALSE, sigma.fix = FALSE,      nu.fix = FALSE) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  identity
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  233.888     10.175  22.987  < 2e-16 ***
Days          10.013      1.365   7.333 7.87e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.7680     0.1136   33.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Nu link function:  log 
Nu Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    2.970      0.548    5.42 1.94e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
No. of observations in the fit:  180 
Degrees of Freedom for the fit:  4
      Residual Deg. of Freedom:  176 
                      at cycle:  20 
 
Global Deviance:     1900.022 
            AIC:     1908.022 
            SBC:     1920.794 
******************************************************************

• 結果を可視化します。

plot(sleep.model)

******************************************************************
          Summary of the Quantile Residuals
                           mean   =  0.001632282 
                       variance   =  1.006378 
               coef. of skewness  =  -0.09752817 
               coef. of kurtosis  =  3.193978 
Filliben correlation coefficient  =  0.9929095 
******************************************************************

• ②brmsパッケージを用いてベイズ推定する方法

set.seed(1234)
sleep.exGauss.brms <- brm(Reaction ~ Days, data = sleepstudy,
                         family = exgaussian(),
                         prior = c(prior('normal(1, 1)', class = 'Intercept'),
                                   prior('normal(1, 1)', class = 'sigma'),
                                   prior('normal(1, 1)', class = 'beta')),
                         seed = 1, chains = 1)

Compiling the C++ model
Start sampling

SAMPLING FOR MODEL '64e67dadfe14ae649762403a4b205f37' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 5.7e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.656754 seconds (Warm-up)
Chain 1:                0.177539 seconds (Sampling)
Chain 1:                0.834293 seconds (Total)
Chain 1: 

There were 525 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems

summary(sleep.exGauss.brms)

There were 525 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: exgaussian 
  Links: mu = identity; sigma = identity; beta = identity 
Formula: Reaction ~ Days 
   Data: sleepstudy (Number of observations: 180) 
Samples: 1 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 1000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   -18.65    137.48  -246.51   214.72        168 1.00
Days          5.47     30.56   -46.55    56.24        168 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.49      0.81     0.24     3.33        231 1.00
beta     36.26      0.57    35.06    37.41        380 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).