学校データ(クラスサイズとSES)のデータの整形

データを読み込む

library(readxl)
hensei2012 <- read_excel("191023_Hensei2012.xlsx", sheet = 1)
setwd("/Users/koyo/Dropbox/000078_CSKAKEN/01_CSNC")
hensei2013 <- read_excel("sho_csnc.xlsx", sheet = 1)
setwd("/Users/koyo/Dropbox/000078_CSKAKEN/191023_CS_SES")

データを整形する

colnames(hensei2012) <- c("shu", "city", "schl", "schlname", "n.g1", 
                          "v01", "v02", "nc.g1", "csmean.g1")
hensei2012 <- subset(hensei2012, v01 != "対象外") 
hensei2012 <- subset(hensei2012, v01 != "複式")
hensei2012 <- hensei2012[c("schlname", "nc.g1", "csmean.g1")]
hensei2012$nc.g1 <- as.numeric(hensei2012$nc.g1)
hensei2012$csmean.g1 <- as.numeric(hensei2012$csmean.g1)

hensei2013 <- hensei2013[,-c(9:10)]

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
csnc.6y <- full_join(hensei2013, hensei2012, by = "schlname")
csnc.6y <- csnc.6y[,c(1:8, 19:20, 9:18)]

SESデータを結合する

setwd("/Users/koyo/Dropbox/191004_SES_Trial/000_Gakku_SES")
ses <- read.csv("Gakku_SES.csv", fileEncoding = "Shift_JIS", stringsAsFactors=FALSE)
setwd("/Users/koyo/Dropbox/000078_CSKAKEN/191023_CS_SES")
ses <- ses[c(5:13)]
colnames(ses) <- c("schlname", "setai.n", "nenshu.m", "nenshu.sd", 
                   "sotsu.n", "tan.n", "dai.n", "tan.p", "dai.p")

csnc.ses <- full_join(csnc.6y, ses, by = "schlname")
write.csv(csnc.ses, "csnc.ses.csv", fileEncoding = "Shift_JIS")

学力データの整形

NRTのデータを読み込む

setwd("/Users/koyo/Dropbox/000078_CSKAKEN/04_NRT/SY2018")
nrt.all <- read.csv("nrt.csv", fileEncoding = "Shift_JIS")
setwd("/Users/koyo/Dropbox/000078_CSKAKEN/191023_CS_SES")

国語だけを取り出す

nrt.koku <- nrt.all %>%
            dplyr::select(c("renban", "sho.sid",
                            "koku.ss.1213", #1-2年生追加
                            "koku.ss.1314", 
                            "koku.ss.1415", 
                            "koku.ss.1516",  
                            "koku.ss.1617",
                            "koku.ss.1718" #6-7年生追加
                            ))

colnames(nrt.koku) <- c("renban", "sid",
                        "g12.koku",
                        "g23.koku",
                        "g34.koku",
                        "g45.koku",
                        "g56.koku",
                        "g67.koku"
                        )

nrt.koku <- na.omit(nrt.koku)

nrt.koku <- nrt.koku %>%
  mutate(id = row_number())

社会だけを取り出す

nrt.shak <- nrt.all %>%
            dplyr::select(c("renban", "sho.sid",
                            "shak.ss.1415", 
                            "shak.ss.1516",  
                            "shak.ss.1617",
                            "shak.ss.1718" #6-7年生追加
                            ))

colnames(nrt.shak) <- c("renban", "sid",
                        "g34.shak",
                        "g45.shak",
                        "g56.shak",
                        "g67.shak"
                        )

nrt.shak <- na.omit(nrt.shak)

nrt.shak <- nrt.shak %>%
  mutate(id = row_number())

NRTデータのある学校の番号を取り出す

nrt.koku.schl <- data.frame(unique(nrt.koku$sid)); colnames(nrt.koku.schl) <- c("sid.new")
nrt.shak.schl <- data.frame(unique(nrt.shak$sid)); colnames(nrt.shak.schl) <- c("sid.new")

分析対象校のクラスサイズなどの表を作る

統廃合校を除外するなどの前処理

csnc.ses.taisho <- dplyr::filter(csnc.ses, taisho.g1 == 1 &
                       togo == 0 & nonrt == 0 & fuku == 0)
csnc.ses.taisho <- na.omit(csnc.ses.taisho)

国語,社会の教科別の対象校に限定した分析対象校のクラスサイズなどの表

library(dplyr)
# 国語対象校
koku.csnc.ses_ <- inner_join(nrt.koku.schl, csnc.ses.taisho, by = "sid.new")

## 必要な列だけを取り出して列名を付け替える
koku.csnc.ses <- koku.csnc.ses_[,c(1, 9:28)]
colnames(koku.csnc.ses) <- c("sid", "nc.g1", "cs.g1", "nc.g2", "cs.g2", 
                             "nc.g3", "cs.g3", "nc.g4", "cs.g4", "nc.g5", 
                             "cs.g5", "nc.g6", "cs.g6", 
                             "setai.n", "nenshu.m", "nenshu.sd", "sotsu.n", 
                             "tan.n", "dai.n", "tan.p", "dai.p")


# 社会対象校
shak.csnc.ses_ <- inner_join(nrt.shak.schl, csnc.ses.taisho, by = "sid.new")
## 必要な列だけを取り出して列名を付け替える
shak.csnc.ses <- shak.csnc.ses_[,c(1, 9:28)]
colnames(shak.csnc.ses)  <- c("sid", "nc.g1", "cs.g1", "nc.g2", "cs.g2", 
                             "nc.g3", "cs.g3", "nc.g4", "cs.g4", "nc.g5", 
                             "cs.g5", "nc.g6", "cs.g6", 
                             "setai.n", "nenshu.m", "nenshu.sd", "sotsu.n", 
                             "tan.n", "dai.n", "tan.p", "dai.p")

2019/10/28 修正

1-6年の学年学級数変動のない学校に限定

# 学級数変動なし
# 国語
koku.csnc.ses$nc.d12 <- koku.csnc.ses$nc.g2 - koku.csnc.ses$nc.g1 
koku.csnc.ses$nc.d23 <- koku.csnc.ses$nc.g3 - koku.csnc.ses$nc.g2 
koku.csnc.ses$nc.d34 <- koku.csnc.ses$nc.g4 - koku.csnc.ses$nc.g3 
koku.csnc.ses$nc.d45 <- koku.csnc.ses$nc.g5 - koku.csnc.ses$nc.g4 
koku.csnc.ses$nc.d56 <- koku.csnc.ses$nc.g6 - koku.csnc.ses$nc.g5 

koku.csnc.ses <- subset(koku.csnc.ses, abs(koku.csnc.ses$nc.d12) == 0 &
                                       abs(koku.csnc.ses$nc.d23) == 0 &
                                       abs(koku.csnc.ses$nc.d34) == 0 &
                                       abs(koku.csnc.ses$nc.d45) == 0 &
                                       abs(koku.csnc.ses$nc.d56) == 0)

# 社会
shak.csnc.ses$nc.d34 <- shak.csnc.ses$nc.g4 - shak.csnc.ses$nc.g3 
shak.csnc.ses$nc.d45 <- shak.csnc.ses$nc.g5 - shak.csnc.ses$nc.g4 
shak.csnc.ses$nc.d56 <- shak.csnc.ses$nc.g6 - shak.csnc.ses$nc.g5 

shak.csnc.ses <- subset(shak.csnc.ses, abs(shak.csnc.ses$nc.d34) == 0 &
                                       abs(shak.csnc.ses$nc.d45) == 0 &
                                       abs(shak.csnc.ses$nc.d56) == 0)

各校の1-6年のクラスサイズの平均

koku.csnc.ses <- dplyr::mutate(koku.csnc.ses, 
                                   csm = (cs.g1 + cs.g2 + cs.g3 + cs.g4 + cs.g5 + cs.g6)/6)
shak.csnc.ses <- dplyr::mutate(shak.csnc.ses, 
                                   csm = (cs.g3 + cs.g4 + cs.g5 + cs.g6)/4)

クラスサイズなどの中心化

# 国語
### 中心化
koku.csm <- mean(koku.csnc.ses$csm)

koku.csnc.ses$cs.c <- koku.csnc.ses$csm -  mean(koku.csm)

## 世帯年収を中心化する
koku.nenshu.mean <- mean(koku.csnc.ses$nenshu.m)
koku.csnc.ses$nenshu.c <- koku.csnc.ses$nenshu.m - koku.nenshu.mean

# 短大高専卒業率を中心化する
koku.tan.p.mean <- mean(koku.csnc.ses$tan.p)
koku.csnc.ses$tan.c <- koku.csnc.ses$tan.p - koku.tan.p.mean

# 大学大学院卒業率を中心化する
koku.dai.p.mean <- mean(koku.csnc.ses$dai.p)
koku.csnc.ses$dai.c <- koku.csnc.ses$dai.p - koku.dai.p.mean

# 社会
## クラスサイズを中心化する
shak.csm <- mean(shak.csnc.ses$csm)

shak.csnc.ses$cs.c <- shak.csnc.ses$csm -  mean(shak.csm)

## 世帯年収を中心化する
shak.nenshu.mean <- mean(shak.csnc.ses$nenshu.m)
shak.csnc.ses$nenshu.c <- shak.csnc.ses$nenshu.m - shak.nenshu.mean

# 短大高専卒業率を中心化する
shak.tan.p.mean <- mean(shak.csnc.ses$tan.p)
shak.csnc.ses$tan.c <- shak.csnc.ses$tan.p - shak.tan.p.mean

# 大学大学院卒業率を中心化する
shak.dai.p.mean <- mean(shak.csnc.ses$dai.p)
shak.csnc.ses$dai.c <- shak.csnc.ses$dai.p - shak.dai.p.mean

縦断データの作成 2019/10/28 修正

国語

# 2019/10/28 修正
# 1-2年生
koku.nrt.12_ <- nrt.koku[c("id", "sid", "g12.koku")]
koku.nrt.12_$year <- 0
koku.csnc.1 <- koku.csnc.ses[c("sid",  "nc.g1", "cs.g1", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.csnc.1$cs.d01 <- 0 # ここは1年生時用の特殊処理
koku.nrt.12 <- dplyr::inner_join(koku.csnc.1, koku.nrt.12_, by = "sid")
koku.nrt.12 <- koku.nrt.12[c("id", "sid", "nc.g1" , "cs.g1", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g12.koku", "year")]
colnames(koku.nrt.12) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

#2-3年生
koku.nrt.23_ <- nrt.koku[c("id", "sid", "g23.koku")]
koku.nrt.23_$year <- 1
koku.csnc.2 <- koku.csnc.ses[c("sid", "nc.g2", "cs.g2", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.nrt.23 <- dplyr::inner_join(koku.csnc.2, koku.nrt.23_, by = "sid")
koku.nrt.23 <- koku.nrt.23[c("id", "sid", "nc.g2" , "cs.g2", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g23.koku", "year")]
colnames(koku.nrt.23) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

#3-4年生
koku.nrt.34_ <- nrt.koku[c("id", "sid", "g34.koku")]
koku.nrt.34_$year <- 2
koku.csnc.3 <- koku.csnc.ses[c("sid", "nc.g3", "cs.g3", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.nrt.34 <- dplyr::inner_join(koku.csnc.3, koku.nrt.34_, by = "sid")
koku.nrt.34 <- koku.nrt.34[c("id", "sid", "nc.g3", "cs.g3", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g34.koku", "year")]
colnames(koku.nrt.34) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

#4-5年生
koku.nrt.45_ <- nrt.koku[c("id", "sid", "g45.koku")]
koku.nrt.45_$year <- 3
koku.csnc.4 <- koku.csnc.ses[c("sid", "nc.g4", "cs.g4", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.nrt.45 <- dplyr::inner_join(koku.csnc.4, koku.nrt.45_, by = "sid")
koku.nrt.45 <- koku.nrt.45[c("id", "sid", "nc.g4", "cs.g4", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g45.koku", "year")]
colnames(koku.nrt.45) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

#5-6年生
koku.nrt.56_ <- nrt.koku[c("id", "sid", "g56.koku")]
koku.nrt.56_$year <- 4
koku.csnc.5 <- koku.csnc.ses[c("sid",  "nc.g5", "cs.g5", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.nrt.56 <- dplyr::inner_join(koku.csnc.5, koku.nrt.56_, by = "sid")
koku.nrt.56 <- koku.nrt.56[c("id", "sid", "nc.g5" , "cs.g5", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g56.koku", "year")]
colnames(koku.nrt.56) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

#6-7年生
koku.nrt.67_ <- nrt.koku[c("id", "sid", "g67.koku")]
koku.nrt.67_$year <- 5
koku.csnc.6 <- koku.csnc.ses[c("sid",  "nc.g6", "cs.g6", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
koku.nrt.67 <- dplyr::inner_join(koku.csnc.6, koku.nrt.67_, by = "sid")
koku.nrt.67 <- koku.nrt.67[c("id", "sid", "nc.g6", "cs.g6", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g67.koku", "year")]
colnames(koku.nrt.67) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "koku", "year")

# うぃーんがっしゃん
koku.nrt.lt <- dplyr::bind_rows(koku.nrt.12, koku.nrt.23, koku.nrt.34, koku.nrt.45, koku.nrt.56, koku.nrt.67)

社会

# 2019/10/28 修正
#3-4年生
shak.nrt.34_ <- nrt.shak[c("id", "sid", "g34.shak")]
shak.nrt.34_$year <- 0
shak.csnc.3 <- shak.csnc.ses[c("sid",  "nc.g3", "cs.g3", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
shak.csnc.3$cs.d23 <- 0 # ここは3年生時用の特殊処理
shak.nrt.34 <- dplyr::inner_join(shak.csnc.3, shak.nrt.34_, by = "sid")
shak.nrt.34 <- shak.nrt.34[c( "id", "sid", "nc.g3", "cs.g3", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g34.shak", "year")]
colnames(shak.nrt.34) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "shak", "year")

#4-5年生
shak.nrt.45_ <- nrt.shak[c("id", "sid", "g45.shak")]
shak.nrt.45_$year <- 1
shak.csnc.4 <- shak.csnc.ses[c("sid", "nc.g4", "cs.g4", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
shak.nrt.45 <- dplyr::inner_join(shak.csnc.4, shak.nrt.45_, by = "sid")
shak.nrt.45 <- shak.nrt.45[c("id", "sid", "nc.g4" , "cs.g4", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g45.shak", "year")]
colnames(shak.nrt.45) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "shak", "year")

#5-6年生
shak.nrt.56_ <- nrt.shak[c("id", "sid", "g56.shak")]
shak.nrt.56_$year <- 2
shak.csnc.5 <- shak.csnc.ses[c("sid", "nc.g5", "cs.g5", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
shak.nrt.56 <- dplyr::inner_join(shak.csnc.5, shak.nrt.56_, by = "sid")
shak.nrt.56 <- shak.nrt.56[c("id", "sid", "nc.g5", "cs.g5", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g56.shak", "year")]
colnames(shak.nrt.56) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "shak", "year")

#6-7年生
shak.nrt.67_ <- nrt.shak[c("id", "sid", "g67.shak")]
shak.nrt.67_$year <- 3
shak.csnc.6 <- shak.csnc.ses[c("sid",  "nc.g6", "cs.g6", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c")]
shak.nrt.67 <- dplyr::inner_join(shak.csnc.6, shak.nrt.67_, by = "sid")
shak.nrt.67 <- shak.nrt.67[c("id", "sid", "nc.g6", "cs.g6", "csm", "cs.c", "nenshu.c", "tan.c", "dai.c", "g67.shak", "year")]
colnames(shak.nrt.67) <- c("id", "sid", "nc", "cs", "cs.m", "cs.c", "nenshu.c", "tan.c", "dai.c", "shak", "year")

# うぃーんがっしゃん
shak.nrt.lt <- dplyr::bind_rows(shak.nrt.34, shak.nrt.45, shak.nrt.56, shak.nrt.67)

分析する

library(brms)
## Loading required package: Rcpp
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Loading 'brms' package (version 2.9.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

国語

# クラスサイズだけ投入
koku.res1 <- brm(koku ~ year + 
                           cs.c + 
                           year:cs.c + 
                           (1 + year|id) + (1 + cs.c|sid), 
                         data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234,
               )
## Compiling the C++ model
## Start sampling
## 
## SAMPLING FOR MODEL 'fec4628b2169c0d0c95a80cd2e9c35b0' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.007761 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 77.61 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2074.24 seconds (Warm-up)
## Chain 1:                6657.08 seconds (Sampling)
## Chain 1:                8731.32 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'fec4628b2169c0d0c95a80cd2e9c35b0' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.003232 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 32.32 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2055.95 seconds (Warm-up)
## Chain 2:                6489.97 seconds (Sampling)
## Chain 2:                8545.92 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'fec4628b2169c0d0c95a80cd2e9c35b0' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.003079 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 30.79 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 5210.32 seconds (Warm-up)
## Chain 3:                1697.16 seconds (Sampling)
## Chain 3:                6907.48 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'fec4628b2169c0d0c95a80cd2e9c35b0' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.002891 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 28.91 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1817.41 seconds (Warm-up)
## Chain 4:                61.7885 seconds (Sampling)
## Chain 4:                1879.2 seconds (Total)
## Chain 4:
## Warning: There were 2751 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4758 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.6, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
# モデル1 + 初期値にだけ大卒率が影響
koku.res2<- brm(koku ~ year + 
                         cs.c + 
                         dai.c + 
                         year:cs.c + 
                         (1 + year|id) + (1 + cs.c * dai.c|sid), 
                         data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234
                         )
## Compiling the C++ model
## Start sampling
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.008297 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 82.97 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 3230.83 seconds (Warm-up)
## Chain 1:                2452.03 seconds (Sampling)
## Chain 1:                5682.86 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.004487 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 44.87 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3072.26 seconds (Warm-up)
## Chain 2:                2441.02 seconds (Sampling)
## Chain 2:                5513.27 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.004274 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 42.74 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 4159.54 seconds (Warm-up)
## Chain 3:                2453.63 seconds (Sampling)
## Chain 3:                6613.17 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.004156 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 41.56 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2935.41 seconds (Warm-up)
## Chain 4:                2432.18 seconds (Sampling)
## Chain 4:                5367.59 seconds (Total)
## Chain 4:
## Warning: There were 62 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
# モデル1 + 初期値と推移に大卒率が影響
koku.res3 <- brm(koku ~ year + 
                         cs.c + 
                         dai.c + 
                         year:cs.c + 
                         year:dai.c + 
                         (1 + year|id) + (1 + cs.c * dai.c|sid), 
                         data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234
                         )
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.00856 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 85.6 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 3356.17 seconds (Warm-up)
## Chain 1:                2493.17 seconds (Sampling)
## Chain 1:                5849.34 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.004232 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 42.32 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3400.14 seconds (Warm-up)
## Chain 2:                2509.23 seconds (Sampling)
## Chain 2:                5909.36 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.004244 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 42.44 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 3223.92 seconds (Warm-up)
## Chain 3:                9868 seconds (Sampling)
## Chain 3:                13091.9 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.004323 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 43.23 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 3198.24 seconds (Warm-up)
## Chain 4:                2474.9 seconds (Sampling)
## Chain 4:                5673.14 seconds (Total)
## Chain 4:
## Warning: There were 214 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 2341 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
# 初期値に大卒率とクラスサイズの交互作用が影響
koku.res4 <- brm(koku ~ year + 
                         cs.c + 
                         dai.c + 
                         dai.c:cs.c + 
                         year:cs.c + 
                         year:dai.c + 
                         (1 + year|id) + (1 + cs.c * dai.c|sid), 
                        data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234
                         )
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.005086 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 50.86 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 3649.34 seconds (Warm-up)
## Chain 1:                9673.2 seconds (Sampling)
## Chain 1:                13322.5 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.004401 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 44.01 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3247.52 seconds (Warm-up)
## Chain 2:                9526.89 seconds (Sampling)
## Chain 2:                12774.4 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.004162 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 41.62 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 3403.45 seconds (Warm-up)
## Chain 3:                2556.46 seconds (Sampling)
## Chain 3:                5959.91 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.004416 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 44.16 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 3550.77 seconds (Warm-up)
## Chain 4:                2529.38 seconds (Sampling)
## Chain 4:                6080.16 seconds (Total)
## Chain 4:
## Warning: There were 752 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4266 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
# 推移に大卒率とクラスサイズの交互作用が影響
koku.res5 <- brm(koku ~ year + 
                         cs.c +
                         dai.c + 
                         year:cs.c + 
                         year:dai.c + 
                         year:cs.c:dai.c + 
                         (1 + year|id) + (1 + cs.c * dai.c|sid), 
                         data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234
                         )
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.004631 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 46.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 3424.84 seconds (Warm-up)
## Chain 1:                9548.24 seconds (Sampling)
## Chain 1:                12973.1 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.004151 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 41.51 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3473.97 seconds (Warm-up)
## Chain 2:                2552.15 seconds (Sampling)
## Chain 2:                6026.12 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.004293 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 42.93 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 3392.71 seconds (Warm-up)
## Chain 3:                5050.24 seconds (Sampling)
## Chain 3:                8442.95 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.004326 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 43.26 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 3386.78 seconds (Warm-up)
## Chain 4:                2507.78 seconds (Sampling)
## Chain 4:                5894.57 seconds (Total)
## Chain 4:
## Warning: There were 415 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 2241 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
# 初期値と推移に大卒率とクラスサイズの交互作用が影響
koku.res6 <- brm(koku ~ year + 
                         dai.c + 
                         cs.c +
                         dai.c:cs.c + 
                         year:cs.c + 
                         year:dai.c + 
                         year:cs.c:dai.c + 
                         (1 + year|id) + (1 + cs.c * dai.c|sid), 
                         data = koku.nrt.lt,
               prior   = NULL, 
               chains  = 4,
               iter    = 5000,
               warmup  = 2500,
               see     = 1234
                         )
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.006233 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 62.33 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 4183.83 seconds (Warm-up)
## Chain 1:                2616.81 seconds (Sampling)
## Chain 1:                6800.65 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.004334 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 43.34 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3725.76 seconds (Warm-up)
## Chain 2:                2606.24 seconds (Sampling)
## Chain 2:                6332 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.004272 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 42.72 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 3355.07 seconds (Warm-up)
## Chain 3:                2602.38 seconds (Sampling)
## Chain 3:                5957.44 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'e95c3961c8fbd2991e14129d68c8ae77' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.004613 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 46.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 3579.41 seconds (Warm-up)
## Chain 4:                5007.04 seconds (Sampling)
## Chain 4:                8586.44 seconds (Total)
## Chain 4:
## Warning: There were 187 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
print(koku.res1, digits = 3)
## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results! 
## We recommend running more iterations and/or setting stronger priors.
## Warning: There were 2751 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + year:cs.c + (1 + year | id) + (1 + cs.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample   Rhat
## sd(Intercept)          6.379     1.711    2.822    7.553          2  8.778
## sd(year)               1.148     0.509    0.808    2.036          2 24.094
## cor(Intercept,year)   -0.087     0.118   -0.318    0.034          2  4.327
## 
## ~sid (Number of levels: 102) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          1.856     0.430    1.153    2.548          3 1.932
## sd(cs.c)               0.072     0.055    0.005    0.206         15 1.099
## cor(Intercept,cs.c)    0.014     0.521   -0.738    0.931          5 1.327
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept   53.717     0.246   53.343   54.241          7 1.178
## year        -0.212     0.027   -0.253   -0.157          8 1.158
## cs.c        -0.088     0.038   -0.151   -0.005         11 1.124
## year:cs.c   -0.007     0.005   -0.017    0.001          6 1.244
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample   Rhat
## sigma    4.695     0.777    4.200    6.076          2 33.454
## 
## 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).
print(koku.res2, digits = 3)
## Warning: There were 62 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          7.349     0.107    7.140    7.562       4351 1.001
## sd(year)               0.855     0.026    0.804    0.905       3359 1.001
## cor(Intercept,year)   -0.022     0.030   -0.080    0.039       7128 1.000
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                1.673     0.291    1.081    2.223       1422
## sd(cs.c)                     0.090     0.065    0.004    0.239       1212
## sd(dai.c)                    7.521     6.323    0.321   23.931       4685
## sd(cs.c:dai.c)               2.511     1.870    0.106    7.021       2032
## cor(Intercept,cs.c)          0.083     0.379   -0.672    0.789       4851
## cor(Intercept,dai.c)        -0.039     0.434   -0.807    0.780      11347
## cor(cs.c,dai.c)             -0.046     0.446   -0.825    0.793      10731
## cor(Intercept,cs.c:dai.c)   -0.228     0.417   -0.862    0.679       5112
## cor(cs.c,cs.c:dai.c)         0.001     0.438   -0.795    0.792       7239
## cor(dai.c,cs.c:dai.c)       -0.013     0.447   -0.809    0.803       7737
##                            Rhat
## sd(Intercept)             1.002
## sd(cs.c)                  1.001
## sd(dai.c)                 1.001
## sd(cs.c:dai.c)            1.002
## cor(Intercept,cs.c)       1.000
## cor(Intercept,dai.c)      1.000
## cor(cs.c,dai.c)           1.001
## cor(Intercept,cs.c:dai.c) 1.000
## cor(cs.c,cs.c:dai.c)      1.000
## cor(dai.c,cs.c:dai.c)     1.000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept   53.787     0.246   53.304   54.276       5962 1.000
## year        -0.205     0.026   -0.258   -0.153      14065 1.000
## cs.c        -0.078     0.045   -0.167    0.010       6158 1.003
## dai.c       -1.697     9.377  -20.071   16.642       5634 1.001
## year:cs.c   -0.009     0.004   -0.017    0.000      14379 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    4.247     0.025    4.197    4.297       5727 1.000
## 
## 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).
print(koku.res3, digits = 3)
## Warning: There were 214 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + year:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          7.347     0.107    7.139    7.560       4147 1.000
## sd(year)               0.851     0.026    0.800    0.902       3360 1.000
## cor(Intercept,year)   -0.019     0.031   -0.079    0.042       7375 1.000
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                1.676     0.284    1.106    2.227       1922
## sd(cs.c)                     0.089     0.063    0.004    0.233       1666
## sd(dai.c)                    7.636     6.456    0.300   24.084       5085
## sd(cs.c:dai.c)               2.473     1.840    0.103    6.779       2708
## cor(Intercept,cs.c)          0.096     0.382   -0.685    0.797       4579
## cor(Intercept,dai.c)        -0.041     0.436   -0.815    0.782      11481
## cor(cs.c,dai.c)             -0.043     0.447   -0.825    0.796      10899
## cor(Intercept,cs.c:dai.c)   -0.224     0.416   -0.865    0.680       6754
## cor(cs.c,cs.c:dai.c)        -0.013     0.447   -0.826    0.808       8273
## cor(dai.c,cs.c:dai.c)       -0.006     0.448   -0.819    0.804       8620
##                            Rhat
## sd(Intercept)             1.002
## sd(cs.c)                  1.001
## sd(dai.c)                 1.000
## sd(cs.c:dai.c)            1.000
## cor(Intercept,cs.c)       1.001
## cor(Intercept,dai.c)      1.000
## cor(cs.c,dai.c)           1.000
## cor(Intercept,cs.c:dai.c) 1.000
## cor(cs.c,cs.c:dai.c)      1.001
## cor(dai.c,cs.c:dai.c)     1.000
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept    53.817     0.244   53.328   54.289       8044 1.000
## year         -0.224     0.027   -0.277   -0.170      12300 1.000
## cs.c         -0.069     0.046   -0.160    0.021       6970 1.000
## dai.c        -5.821     9.634  -24.617   12.760       7877 1.000
## year:cs.c    -0.015     0.005   -0.024   -0.006      13730 1.000
## year:dai.c    2.722     0.808    1.127    4.282      13347 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    4.247     0.025    4.198    4.297       5063 1.000
## 
## 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).
print(koku.res4, digits = 3)
## Warning: There were 752 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + dai.c:cs.c + year:cs.c + year:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          7.352     0.106    7.143    7.564       4292 1.000
## sd(year)               0.852     0.026    0.800    0.901       3033 1.001
## cor(Intercept,year)   -0.020     0.031   -0.080    0.041       7028 1.000
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                1.675     0.293    1.065    2.242       1628
## sd(cs.c)                     0.092     0.064    0.004    0.233       1503
## sd(dai.c)                    7.685     6.444    0.274   24.122       5374
## sd(cs.c:dai.c)               2.599     1.875    0.114    6.918       2382
## cor(Intercept,cs.c)          0.099     0.380   -0.656    0.797       4124
## cor(Intercept,dai.c)        -0.050     0.433   -0.817    0.777       9510
## cor(cs.c,dai.c)             -0.046     0.448   -0.830    0.798       9915
## cor(Intercept,cs.c:dai.c)   -0.224     0.415   -0.866    0.668       5305
## cor(cs.c,cs.c:dai.c)        -0.007     0.439   -0.803    0.804       7164
## cor(dai.c,cs.c:dai.c)       -0.012     0.448   -0.811    0.807       6691
##                            Rhat
## sd(Intercept)             1.002
## sd(cs.c)                  1.001
## sd(dai.c)                 1.000
## sd(cs.c:dai.c)            1.001
## cor(Intercept,cs.c)       1.000
## cor(Intercept,dai.c)      1.000
## cor(cs.c,dai.c)           1.000
## cor(Intercept,cs.c:dai.c) 1.000
## cor(cs.c,cs.c:dai.c)      1.000
## cor(dai.c,cs.c:dai.c)     1.000
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept    53.790     0.269   53.265   54.325       6839 1.001
## year         -0.222     0.027   -0.274   -0.170      10863 1.000
## cs.c         -0.067     0.046   -0.158    0.021       6379 1.000
## dai.c        -6.632    10.617  -28.023   13.824       6561 1.001
## cs.c:dai.c    0.387     1.709   -2.929    3.819       7047 1.000
## year:cs.c    -0.015     0.005   -0.024   -0.006      11354 1.000
## year:dai.c    2.717     0.818    1.120    4.323      10599 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    4.247     0.025    4.197    4.297       4899 1.001
## 
## 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).
print(koku.res5, digits = 3)
## Warning: There were 415 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + year:dai.c + year:cs.c:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          7.350     0.106    7.147    7.562       4110 1.001
## sd(year)               0.852     0.026    0.802    0.901       3455 1.000
## cor(Intercept,year)   -0.020     0.030   -0.078    0.040       7591 1.000
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                1.672     0.278    1.116    2.216       1879
## sd(cs.c)                     0.089     0.063    0.003    0.230       1544
## sd(dai.c)                    7.740     6.493    0.276   24.358       5567
## sd(cs.c:dai.c)               2.510     1.897    0.091    7.055       1658
## cor(Intercept,cs.c)          0.099     0.381   -0.662    0.810       4609
## cor(Intercept,dai.c)        -0.036     0.432   -0.812    0.778       8512
## cor(cs.c,dai.c)             -0.053     0.447   -0.839    0.797      10392
## cor(Intercept,cs.c:dai.c)   -0.223     0.414   -0.856    0.673       5893
## cor(cs.c,cs.c:dai.c)        -0.018     0.439   -0.814    0.796       7526
## cor(dai.c,cs.c:dai.c)       -0.012     0.452   -0.835    0.812       6714
##                            Rhat
## sd(Intercept)             1.003
## sd(cs.c)                  1.001
## sd(dai.c)                 1.001
## sd(cs.c:dai.c)            1.003
## cor(Intercept,cs.c)       1.001
## cor(Intercept,dai.c)      1.000
## cor(cs.c,dai.c)           1.000
## cor(Intercept,cs.c:dai.c) 1.001
## cor(cs.c,cs.c:dai.c)      1.000
## cor(dai.c,cs.c:dai.c)     1.000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept         53.814     0.247   53.332   54.297       6513 1.000
## year              -0.228     0.029   -0.283   -0.171       8745 1.000
## cs.c              -0.067     0.045   -0.158    0.021       7051 1.000
## dai.c             -5.687     9.593  -24.736   12.778       7624 1.000
## year:cs.c         -0.015     0.005   -0.024   -0.006      11718 1.000
## year:dai.c         2.301     1.051    0.217    4.346      10211 1.001
## year:cs.c:dai.c    0.100     0.161   -0.218    0.415      10034 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    4.246     0.026    4.196    4.297       4579 1.000
## 
## 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).
print(koku.res6, digits = 3)
## Warning: There were 187 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + dai.c + cs.c + dai.c:cs.c + year:cs.c + year:dai.c + year:cs.c:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)          7.351     0.106    7.148    7.559       4287 1.000
## sd(year)               0.852     0.026    0.801    0.902       3361 1.002
## cor(Intercept,year)   -0.020     0.030   -0.078    0.040       6580 1.000
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                1.681     0.285    1.111    2.238       1395
## sd(cs.c)                     0.087     0.062    0.003    0.227       1206
## sd(dai.c)                    7.651     6.316    0.302   23.567       4676
## sd(cs.c:dai.c)               2.642     1.941    0.107    7.131       1796
## cor(Intercept,cs.c)          0.098     0.382   -0.672    0.794       4363
## cor(Intercept,dai.c)        -0.049     0.430   -0.807    0.774       9417
## cor(cs.c,dai.c)             -0.049     0.443   -0.823    0.775      10805
## cor(Intercept,cs.c:dai.c)   -0.219     0.414   -0.856    0.682       5354
## cor(cs.c,cs.c:dai.c)        -0.009     0.437   -0.796    0.796       7119
## cor(dai.c,cs.c:dai.c)       -0.005     0.444   -0.798    0.803       6988
##                            Rhat
## sd(Intercept)             1.005
## sd(cs.c)                  1.003
## sd(dai.c)                 1.001
## sd(cs.c:dai.c)            1.002
## cor(Intercept,cs.c)       1.001
## cor(Intercept,dai.c)      1.000
## cor(cs.c,dai.c)           1.000
## cor(Intercept,cs.c:dai.c) 1.000
## cor(cs.c,cs.c:dai.c)      1.000
## cor(dai.c,cs.c:dai.c)     1.001
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept         53.803     0.270   53.278   54.337       6355 1.000
## year              -0.229     0.028   -0.283   -0.173      11823 1.000
## dai.c             -6.483    10.747  -27.710   14.981       5493 1.000
## cs.c              -0.065     0.046   -0.156    0.024       5410 1.000
## dai.c:cs.c         0.236     1.792   -3.300    3.720       4992 1.000
## year:cs.c         -0.015     0.005   -0.024   -0.005      12588 1.000
## year:dai.c         2.327     1.086    0.211    4.448      12636 1.000
## year:dai.c:cs.c    0.099     0.165   -0.221    0.424      11296 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    4.247     0.026    4.196    4.298       5082 1.001
## 
## 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).
print(koku.res1, digits = 2)
## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results! 
## We recommend running more iterations and/or setting stronger priors.
## Warning: There were 2751 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + year:cs.c + (1 + year | id) + (1 + cs.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)           6.38      1.71     2.82     7.55          2  8.78
## sd(year)                1.15      0.51     0.81     2.04          2 24.09
## cor(Intercept,year)    -0.09      0.12    -0.32     0.03          2  4.33
## 
## ~sid (Number of levels: 102) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           1.86      0.43     1.15     2.55          3 1.93
## sd(cs.c)                0.07      0.06     0.00     0.21         15 1.10
## cor(Intercept,cs.c)     0.01      0.52    -0.74     0.93          5 1.33
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept    53.72      0.25    53.34    54.24          7 1.18
## year         -0.21      0.03    -0.25    -0.16          8 1.16
## cs.c         -0.09      0.04    -0.15    -0.00         11 1.12
## year:cs.c    -0.01      0.00    -0.02     0.00          6 1.24
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma     4.69      0.78     4.20     6.08          2 33.45
## 
## 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).
print(koku.res2, digits = 2)
## Warning: There were 62 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           7.35      0.11     7.14     7.56       4351 1.00
## sd(year)                0.86      0.03     0.80     0.91       3359 1.00
## cor(Intercept,year)    -0.02      0.03    -0.08     0.04       7128 1.00
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                 1.67      0.29     1.08     2.22       1422
## sd(cs.c)                      0.09      0.06     0.00     0.24       1212
## sd(dai.c)                     7.52      6.32     0.32    23.93       4685
## sd(cs.c:dai.c)                2.51      1.87     0.11     7.02       2032
## cor(Intercept,cs.c)           0.08      0.38    -0.67     0.79       4851
## cor(Intercept,dai.c)         -0.04      0.43    -0.81     0.78      11347
## cor(cs.c,dai.c)              -0.05      0.45    -0.83     0.79      10731
## cor(Intercept,cs.c:dai.c)    -0.23      0.42    -0.86     0.68       5112
## cor(cs.c,cs.c:dai.c)          0.00      0.44    -0.80     0.79       7239
## cor(dai.c,cs.c:dai.c)        -0.01      0.45    -0.81     0.80       7737
##                           Rhat
## sd(Intercept)             1.00
## sd(cs.c)                  1.00
## sd(dai.c)                 1.00
## sd(cs.c:dai.c)            1.00
## cor(Intercept,cs.c)       1.00
## cor(Intercept,dai.c)      1.00
## cor(cs.c,dai.c)           1.00
## cor(Intercept,cs.c:dai.c) 1.00
## cor(cs.c,cs.c:dai.c)      1.00
## cor(dai.c,cs.c:dai.c)     1.00
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept    53.79      0.25    53.30    54.28       5962 1.00
## year         -0.21      0.03    -0.26    -0.15      14065 1.00
## cs.c         -0.08      0.05    -0.17     0.01       6158 1.00
## dai.c        -1.70      9.38   -20.07    16.64       5634 1.00
## year:cs.c    -0.01      0.00    -0.02     0.00      14379 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     4.25      0.03     4.20     4.30       5727 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).
print(koku.res3, digits = 2)
## Warning: There were 214 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + year:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           7.35      0.11     7.14     7.56       4147 1.00
## sd(year)                0.85      0.03     0.80     0.90       3360 1.00
## cor(Intercept,year)    -0.02      0.03    -0.08     0.04       7375 1.00
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                 1.68      0.28     1.11     2.23       1922
## sd(cs.c)                      0.09      0.06     0.00     0.23       1666
## sd(dai.c)                     7.64      6.46     0.30    24.08       5085
## sd(cs.c:dai.c)                2.47      1.84     0.10     6.78       2708
## cor(Intercept,cs.c)           0.10      0.38    -0.68     0.80       4579
## cor(Intercept,dai.c)         -0.04      0.44    -0.81     0.78      11481
## cor(cs.c,dai.c)              -0.04      0.45    -0.82     0.80      10899
## cor(Intercept,cs.c:dai.c)    -0.22      0.42    -0.86     0.68       6754
## cor(cs.c,cs.c:dai.c)         -0.01      0.45    -0.83     0.81       8273
## cor(dai.c,cs.c:dai.c)        -0.01      0.45    -0.82     0.80       8620
##                           Rhat
## sd(Intercept)             1.00
## sd(cs.c)                  1.00
## sd(dai.c)                 1.00
## sd(cs.c:dai.c)            1.00
## cor(Intercept,cs.c)       1.00
## cor(Intercept,dai.c)      1.00
## cor(cs.c,dai.c)           1.00
## cor(Intercept,cs.c:dai.c) 1.00
## cor(cs.c,cs.c:dai.c)      1.00
## cor(dai.c,cs.c:dai.c)     1.00
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept     53.82      0.24    53.33    54.29       8044 1.00
## year          -0.22      0.03    -0.28    -0.17      12300 1.00
## cs.c          -0.07      0.05    -0.16     0.02       6970 1.00
## dai.c         -5.82      9.63   -24.62    12.76       7877 1.00
## year:cs.c     -0.02      0.00    -0.02    -0.01      13730 1.00
## year:dai.c     2.72      0.81     1.13     4.28      13347 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     4.25      0.03     4.20     4.30       5063 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).
print(koku.res4, digits = 2)
## Warning: There were 752 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + dai.c:cs.c + year:cs.c + year:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           7.35      0.11     7.14     7.56       4292 1.00
## sd(year)                0.85      0.03     0.80     0.90       3033 1.00
## cor(Intercept,year)    -0.02      0.03    -0.08     0.04       7028 1.00
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                 1.68      0.29     1.06     2.24       1628
## sd(cs.c)                      0.09      0.06     0.00     0.23       1503
## sd(dai.c)                     7.69      6.44     0.27    24.12       5374
## sd(cs.c:dai.c)                2.60      1.87     0.11     6.92       2382
## cor(Intercept,cs.c)           0.10      0.38    -0.66     0.80       4124
## cor(Intercept,dai.c)         -0.05      0.43    -0.82     0.78       9510
## cor(cs.c,dai.c)              -0.05      0.45    -0.83     0.80       9915
## cor(Intercept,cs.c:dai.c)    -0.22      0.42    -0.87     0.67       5305
## cor(cs.c,cs.c:dai.c)         -0.01      0.44    -0.80     0.80       7164
## cor(dai.c,cs.c:dai.c)        -0.01      0.45    -0.81     0.81       6691
##                           Rhat
## sd(Intercept)             1.00
## sd(cs.c)                  1.00
## sd(dai.c)                 1.00
## sd(cs.c:dai.c)            1.00
## cor(Intercept,cs.c)       1.00
## cor(Intercept,dai.c)      1.00
## cor(cs.c,dai.c)           1.00
## cor(Intercept,cs.c:dai.c) 1.00
## cor(cs.c,cs.c:dai.c)      1.00
## cor(dai.c,cs.c:dai.c)     1.00
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept     53.79      0.27    53.26    54.33       6839 1.00
## year          -0.22      0.03    -0.27    -0.17      10863 1.00
## cs.c          -0.07      0.05    -0.16     0.02       6379 1.00
## dai.c         -6.63     10.62   -28.02    13.82       6561 1.00
## cs.c:dai.c     0.39      1.71    -2.93     3.82       7047 1.00
## year:cs.c     -0.02      0.00    -0.02    -0.01      11354 1.00
## year:dai.c     2.72      0.82     1.12     4.32      10599 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     4.25      0.03     4.20     4.30       4899 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).
print(koku.res5, digits = 2)
## Warning: There were 415 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + cs.c + dai.c + year:cs.c + year:dai.c + year:cs.c:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           7.35      0.11     7.15     7.56       4110 1.00
## sd(year)                0.85      0.03     0.80     0.90       3455 1.00
## cor(Intercept,year)    -0.02      0.03    -0.08     0.04       7591 1.00
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                 1.67      0.28     1.12     2.22       1879
## sd(cs.c)                      0.09      0.06     0.00     0.23       1544
## sd(dai.c)                     7.74      6.49     0.28    24.36       5567
## sd(cs.c:dai.c)                2.51      1.90     0.09     7.05       1658
## cor(Intercept,cs.c)           0.10      0.38    -0.66     0.81       4609
## cor(Intercept,dai.c)         -0.04      0.43    -0.81     0.78       8512
## cor(cs.c,dai.c)              -0.05      0.45    -0.84     0.80      10392
## cor(Intercept,cs.c:dai.c)    -0.22      0.41    -0.86     0.67       5893
## cor(cs.c,cs.c:dai.c)         -0.02      0.44    -0.81     0.80       7526
## cor(dai.c,cs.c:dai.c)        -0.01      0.45    -0.83     0.81       6714
##                           Rhat
## sd(Intercept)             1.00
## sd(cs.c)                  1.00
## sd(dai.c)                 1.00
## sd(cs.c:dai.c)            1.00
## cor(Intercept,cs.c)       1.00
## cor(Intercept,dai.c)      1.00
## cor(cs.c,dai.c)           1.00
## cor(Intercept,cs.c:dai.c) 1.00
## cor(cs.c,cs.c:dai.c)      1.00
## cor(dai.c,cs.c:dai.c)     1.00
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept          53.81      0.25    53.33    54.30       6513 1.00
## year               -0.23      0.03    -0.28    -0.17       8745 1.00
## cs.c               -0.07      0.05    -0.16     0.02       7051 1.00
## dai.c              -5.69      9.59   -24.74    12.78       7624 1.00
## year:cs.c          -0.02      0.00    -0.02    -0.01      11718 1.00
## year:dai.c          2.30      1.05     0.22     4.35      10211 1.00
## year:cs.c:dai.c     0.10      0.16    -0.22     0.41      10034 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     4.25      0.03     4.20     4.30       4579 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).
print(koku.res6, digits = 2)
## Warning: There were 187 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: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: koku ~ year + dai.c + cs.c + dai.c:cs.c + year:cs.c + year:dai.c + year:cs.c:dai.c + (1 + year | id) + (1 + cs.c * dai.c | sid) 
##    Data: koku.nrt.lt (Number of observations: 20718) 
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 3453) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           7.35      0.11     7.15     7.56       4287 1.00
## sd(year)                0.85      0.03     0.80     0.90       3361 1.00
## cor(Intercept,year)    -0.02      0.03    -0.08     0.04       6580 1.00
## 
## ~sid (Number of levels: 102) 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                 1.68      0.29     1.11     2.24       1395
## sd(cs.c)                      0.09      0.06     0.00     0.23       1206
## sd(dai.c)                     7.65      6.32     0.30    23.57       4676
## sd(cs.c:dai.c)                2.64      1.94     0.11     7.13       1796
## cor(Intercept,cs.c)           0.10      0.38    -0.67     0.79       4363
## cor(Intercept,dai.c)         -0.05      0.43    -0.81     0.77       9417
## cor(cs.c,dai.c)              -0.05      0.44    -0.82     0.78      10805
## cor(Intercept,cs.c:dai.c)    -0.22      0.41    -0.86     0.68       5354
## cor(cs.c,cs.c:dai.c)         -0.01      0.44    -0.80     0.80       7119
## cor(dai.c,cs.c:dai.c)        -0.00      0.44    -0.80     0.80       6988
##                           Rhat
## sd(Intercept)             1.00
## sd(cs.c)                  1.00
## sd(dai.c)                 1.00
## sd(cs.c:dai.c)            1.00
## cor(Intercept,cs.c)       1.00
## cor(Intercept,dai.c)      1.00
## cor(cs.c,dai.c)           1.00
## cor(Intercept,cs.c:dai.c) 1.00
## cor(cs.c,cs.c:dai.c)      1.00
## cor(dai.c,cs.c:dai.c)     1.00
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept          53.80      0.27    53.28    54.34       6355 1.00
## year               -0.23      0.03    -0.28    -0.17      11823 1.00
## dai.c              -6.48     10.75   -27.71    14.98       5493 1.00
## cs.c               -0.07      0.05    -0.16     0.02       5410 1.00
## dai.c:cs.c          0.24      1.79    -3.30     3.72       4992 1.00
## year:cs.c          -0.01      0.00    -0.02    -0.01      12588 1.00
## year:dai.c          2.33      1.09     0.21     4.45      12636 1.00
## year:dai.c:cs.c     0.10      0.16    -0.22     0.42      11296 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     4.25      0.03     4.20     4.30       5082 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).
brms::waic(koku.res1)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -64276.3 140.6
## p_waic      5393.3  81.1
## waic      128552.6 281.2
## Warning: 3116 (15.0%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::waic(koku.res2)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -61778.2 126.9
## p_waic      4093.6  45.9
## waic      123556.4 253.8
## Warning: 2594 (12.5%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::waic(koku.res3)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -61775.7 126.9
## p_waic      4085.7  45.8
## waic      123551.4 253.7
## Warning: 2591 (12.5%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::waic(koku.res4)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -61776.1 126.8
## p_waic      4088.7  45.8
## waic      123552.3 253.6
## Warning: 2599 (12.5%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::waic(koku.res5)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -61771.3 126.9
## p_waic      4085.9  45.8
## waic      123542.6 253.7
## Warning: 2579 (12.4%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::waic(koku.res6)
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -61776.0 126.9
## p_waic      4088.4  45.9
## waic      123552.1 253.8
## Warning: 2594 (12.5%) p_waic estimates greater than 0.4. We recommend
## trying loo instead.
brms::loo(koku.res1)
## Warning: Found 1861 observations with a pareto_k > 0.7 in model
## 'koku.res1'. With this many problematic observations, it may be more
## appropriate to use 'kfold' with argument 'K = 10' to perform 10-fold cross-
## validation rather than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -64244.9 130.1
## p_loo      5361.9  65.7
## looic    128489.8 260.2
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     15607 75.3%   0         
##  (0.5, 0.7]   (ok)        3250 15.7%   0         
##    (0.7, 1]   (bad)       1326  6.4%   0         
##    (1, Inf)   (very bad)   535  2.6%   0         
## See help('pareto-k-diagnostic') for details.
brms::loo(koku.res2)
## Warning: Found 124 observations with a pareto_k > 0.7 in model 'koku.res2'.
## With this many problematic observations, it may be more appropriate to use
## 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather
## than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -62000.8 129.1
## p_loo      4316.2  48.5
## looic    124001.5 258.3
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     19488 94.1%   534       
##  (0.5, 0.7]   (ok)        1106  5.3%   93        
##    (0.7, 1]   (bad)        120  0.6%   15        
##    (1, Inf)   (very bad)     4  0.0%   8         
## See help('pareto-k-diagnostic') for details.
brms::loo(koku.res3)
## Warning: Found 121 observations with a pareto_k > 0.7 in model 'koku.res3'.
## With this many problematic observations, it may be more appropriate to use
## 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather
## than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -61995.9 129.0
## p_loo      4305.9  48.3
## looic    123991.8 258.0
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     19521 94.2%   329       
##  (0.5, 0.7]   (ok)        1076  5.2%   116       
##    (0.7, 1]   (bad)        118  0.6%   19        
##    (1, Inf)   (very bad)     3  0.0%   9         
## See help('pareto-k-diagnostic') for details.
brms::loo(koku.res4)
## Warning: Found 129 observations with a pareto_k > 0.7 in model 'koku.res4'.
## With this many problematic observations, it may be more appropriate to use
## 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather
## than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -61998.1 129.0
## p_loo      4310.6  48.3
## looic    123996.1 258.0
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     19524 94.2%   284       
##  (0.5, 0.7]   (ok)        1065  5.1%   112       
##    (0.7, 1]   (bad)        126  0.6%   11        
##    (1, Inf)   (very bad)     3  0.0%   27        
## See help('pareto-k-diagnostic') for details.
brms::loo(koku.res5)
## Warning: Found 131 observations with a pareto_k > 0.7 in model 'koku.res5'.
## With this many problematic observations, it may be more appropriate to use
## 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather
## than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -61994.2 129.2
## p_loo      4308.8  48.5
## looic    123988.4 258.4
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     19551 94.4%   384       
##  (0.5, 0.7]   (ok)        1036  5.0%   79        
##    (0.7, 1]   (bad)        123  0.6%   11        
##    (1, Inf)   (very bad)     8  0.0%   5         
## See help('pareto-k-diagnostic') for details.
brms::loo(koku.res6)
## Warning: Found 139 observations with a pareto_k > 0.7 in model 'koku.res6'.
## With this many problematic observations, it may be more appropriate to use
## 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather
## than LOO.
## 
## Computed from 10000 by 20718 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -61999.5 129.2
## p_loo      4311.8  48.6
## looic    123999.0 258.5
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     19494 94.1%   408       
##  (0.5, 0.7]   (ok)        1085  5.2%   90        
##    (0.7, 1]   (bad)        134  0.6%   13        
##    (1, Inf)   (very bad)     5  0.0%   2         
## See help('pareto-k-diagnostic') for details.