지난 번의 문제를 다시 한 번 보자.

다음의 R 코드는 데이터 dat를 생성한다. 이때 설명변수는 x1에서 x10은 상관관계 0.1인 표준정규분포를 따른다. 실제 데이터생성모형은 다음과 같다.

\[y = 3 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_{10} x_{10} + e,\ \ e \sim \mathcal{N}(0,1)\]

정확한 모형을 알고 있다면 자료를 모형에 적합하여 모수( \(\beta_1, \beta_2, \cdots\) )을 추정하면 된다. 하지만 문제가 있다. 모수의 갯수가 늘어날 수록 모수를 정확하게 추정하기 위해 필요한 표본의 크기도 증가한다. 아래에서 표본의 크기(N)는 40이다.

set.seed(102)
library(data.table)
library(dplyr)
library(ggplot2)
N <- 40
N2 <- 2*N; 
ncolMatX <- 11

sigma.cor = 0.1

library(mvtnorm)
sigma <- matrix(
  sigma.cor,
  ncolMatX, ncolMatX)
diag(sigma) = 1

matX <- rmvnorm(N2, mean=rep(0,ncolMatX), sigma=sigma)
colnames(matX) <- c(paste0("x", 1:10), "e")
X <- data.table(matX)

paramTrue <- c(3, rnorm(5,1,0.2),
               rnorm(5,-1.4,0.2))
names(paramTrue) <- paste0("beta", 0:10)

lv <- X[, .(f1 = paramTrue['beta1']*x1 + 
              paramTrue['beta2']*x2 + 
              paramTrue['beta3']*x3 + 
              paramTrue['beta4']*x4 + 
              paramTrue['beta5']*x5,
            f2 = paramTrue['beta6']*x6 + 
              paramTrue['beta7']*x7 + 
              paramTrue['beta8']*x8 + 
              paramTrue['beta9']*x9 + 
              paramTrue['beta10']*x10)]

y <- 3 + lv$f1 + lv$f2 + X$e

dfParamTrue <- data.frame(key=factor(names(paramTrue), 
                                     levels= paste0("beta", 0:10)),
                          value=paramTrue)

datAll <- data.table(X,y)
dat <- datAll[1:N]
newdat <- datAll[(N+1):N2]

앞의 글에서 우리는 \(\beta_1, \beta_2, \beta_3, \beta_4, \beta_5\) 가 모두 비슷하고, \(\beta_6, \cdots, \beta_{10}\) 이 모두 비슷하다는 사전 지식을 활용하여 새로운 변수 \(f_1 = \beta_1 + \cdots + \beta_5\) 와 \(f_2=\beta_6 + \cdots +\beta_{10}\) 을 만들었다.

하지만 이 방법은 \(\beta_1 = \beta_2 = \cdots = \beta_5\) 라는 가정과 \(\beta_6 = \cdots = \beta_{10}\) 이라는 가정이 모두 성립할 때 가장 좋은 성능을 보일 것이다. 그리고 그 가정에서 어긋나는 정도에 따라 \(f_1\) , \(f_2\) 만의 선형모형은 성능이 저하될 것이다.

만약 \(\beta_1 = \beta_2 = \cdots = \beta_5\) 과 \(\beta_6 = \cdots = \beta_{10}\) 의 강한 가정을 조금 누그러트린다면 어떨까?

다음의 블로그는 The Best Of Both Worlds를 얘기하고 있다. 다시 말해 모든 베타를 자유롭게 두는 방법과 한 그룹에 속하는 모든 베타를 동일하게 설정하는 두 가지 방법의 장점을 모두 취하는 방법이 존재한다!

베이지안 방법

다음의 베이지안 모형은 \(\beta_1, \beta_2, \beta_3, \beta_4, \beta_5\) 가 모두 서로 비슷하고, \(\beta_6, \cdots, \beta_{10}\) 이 모두 서로 비슷하다는 사전 지식을 활용하는 모형이다.

\[y = 3 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_{10} x_{10} + e,\ \ e \sim \mathcal{N}(0,1)\]

• 사전 분포
  • \(\beta_1, \beta_2, \cdots, \beta_5 \sim \mathcal{N}(\mu_1, \sigma^2_1)\)
  • \(\beta_6, \beta_7, \cdots, \beta_{10} \sim \mathcal{N}(\mu_2, \sigma^2_2)\)
  • \(\mu_1, \mu_2 \sim \textrm{Unif}(-100,100)\)
  • \(\sigma^2_1, \sigma^2_2 \sim \textrm{Unif}(0, 100)\)

다음의 코드는 이를 구현한다. (JAGS로 MCMC 샘플링을 하기 전에 효율적인 샘플링을 위해 표준화(scale)를 실시했다.)

jags.model <- "


model {
  for (i in 1:N) {
    ymean[i] <- beta0 + beta1*x1[i] + beta2*x2[i] + beta3*x3[i] + beta4*x4[i] + beta5*x5[i] + 
            beta6*x6[i] + beta7*x7[i] + beta8*x8[i] + beta9*x9[i] + beta10*x10[i]
    y[i] ~ dnorm(ymean[i], tau)
  }

  # priors
  tau <- pow(sigma2, -2)
  sigma2 ~ dunif(0,100)
  beta0 ~ dnorm(0, 0.01)
  beta1 ~ dnorm(betaAmean, betaAtau)
  beta2 ~ dnorm(betaAmean, betaAtau)
  beta3 ~ dnorm(betaAmean, betaAtau)
  beta4 ~ dnorm(betaAmean, betaAtau)
  beta5 ~ dnorm(betaAmean, betaAtau)
  beta6 ~ dnorm(betaBmean, betaBtau)
  beta7 ~ dnorm(betaBmean, betaBtau)
  beta8 ~ dnorm(betaBmean, betaBtau)
  beta9 ~ dnorm(betaBmean, betaBtau)
  beta10 ~ dnorm(betaBmean, betaBtau)
  betaAmean ~ dunif(-100,100)
  betaAtau <- pow(betaAsigma2, -2)
  betaBmean ~ dunif(-100,100)
  betaBtau <- pow(betaBsigma2, -2)
  betaAsigma2 ~ dunif(0,100)
  betaBsigma2 ~ dunif(0,100)
}
"

cat(jags.model, file='currentModel.j')
#cat(jags.model2, file='currentModel.j')

require(rjags)

dat_ <- (dat_mat <- 
  dat %>% select(-e) %>% scale) %>% data.frame

N = nrow(dat_)


myj <- jags.model('currentModel.j',
                  data = c(as.list(dat_), 'N'=N),
                  n.chains = 3,
                  quiet = TRUE)

update(myj, n.iter=2000, progress.bar = 'none')

n.thin=5


mcmcsamp <- 
  coda.samples(myj, 
               c('beta0', 'beta1', 'beta2', 'beta3', 'beta4', 'beta5',
                 'beta6', 'beta7', 'beta8', 'beta9', 'beta10',
                 'betaAmean', 'betaBmean', 'betaAsigma2', 'betaBsigma2', 'sigma2'),
               1000*n.thin, thin=n.thin,
               progress.bar = 'none')  

우선 MCMC 샘플의 수렴성을 확인해야 한다. 보통 Rhat이 1.2보다 작으면 수렴한 것으로 본다.

#summary(mcmcsamp)
#plot(mcmcsamp)

#mcmcr::rhat(mcmcsamp)
#mcmcr::rhat(mcmcsamp, by='parameter', as_df=TRUE)

print(MCMCvis::MCMCsummary(mcmcsamp, n.eff=TRUE)) # install.packages('MCMCvis')

##                      mean         sd         2.5%           50%
## beta0       -0.0008153451 0.04730933 -0.093582641  0.0000556605
## beta1        0.3002368789 0.03549022  0.237489178  0.2981886062
## beta10      -0.3662336392 0.04543191 -0.455571464 -0.3664033830
## beta2        0.3097945925 0.03735937  0.243492538  0.3062498109
## beta3        0.3044016981 0.03615626  0.235866801  0.3021763979
## beta4        0.2887275983 0.03600711  0.214269698  0.2897465993
## beta5        0.2616867426 0.04282079  0.163398803  0.2679280557
## beta6       -0.2482596931 0.05724510 -0.359905930 -0.2478469606
## beta7       -0.3653164949 0.04829268 -0.459222861 -0.3656849740
## beta8       -0.4928264247 0.05339144 -0.599683460 -0.4940052212
## beta9       -0.4181135511 0.04756135 -0.515210669 -0.4167642616
## betaAmean    0.2925825249 0.03680493  0.223940666  0.2923181691
## betaAsigma2  0.0522688007 0.05154199  0.001763328  0.0384909401
## betaBmean   -0.3777158395 0.08670358 -0.540695126 -0.3791240688
## betaBsigma2  0.1553488341 0.11108145  0.040603461  0.1263819110
## sigma2       0.2924111077 0.03858603  0.227463442  0.2883110324
##                   97.5% Rhat n.eff
## beta0        0.09180819 1.00  3195
## beta1        0.37752391 1.00  1932
## beta10      -0.27775564 1.00  2894
## beta2        0.39141868 1.00  1870
## beta3        0.38345474 1.00  1857
## beta4        0.36016644 1.00  2468
## beta5        0.33151414 1.00  1547
## beta6       -0.13748495 1.00  2113
## beta7       -0.27138717 1.00  2753
## beta8       -0.38564016 1.00  2669
## beta9       -0.32821133 1.00  3512
## betaAmean    0.36672810 1.01  2273
## betaAsigma2  0.18400535 1.01   941
## betaBmean   -0.20465755 1.00  3077
## betaBsigma2  0.46000681 1.00  1265
## sigma2       0.37912164 1.00  2892

print(gelman.diag(mcmcsamp)) # coda::gelman.diag

## Potential scale reduction factors:
## 
##             Point est. Upper C.I.
## beta0            1.000       1.00
## beta1            1.001       1.00
## beta10           0.999       1.00
## beta2            1.000       1.00
## beta3            1.001       1.00
## beta4            1.000       1.00
## beta5            1.000       1.00
## beta6            1.000       1.00
## beta7            1.003       1.01
## beta8            1.002       1.01
## beta9            1.000       1.00
## betaAmean        1.005       1.01
## betaAsigma2      1.007       1.01
## betaBmean        1.003       1.00
## betaBsigma2      1.001       1.00
## sigma2           1.001       1.00
## 
## Multivariate psrf
## 
## 1.01

#https://www.reddit.com/r/rstats/comments/884rlc/how_to_calculate_gelmanrubin_convergence_from/
#https://avehtari.github.io/rhat_ess/rhat_ess.html

이제 새로운 데이터에 적용하여 예측 성능을 비교해보자.

newdat_ <- sweep(newdat %>% select(-y), 2, attr(dat_mat, "scaled:center"), "-")
newdat_ <- sweep(newdat_ , 2, attr(dat_mat, "scaled:scale"), "/")

#mcmcsamp
#str(mcmcsamp)
bparm <- rbind(mcmcsamp[[1]], mcmcsamp[[2]], mcmcsamp[[3]])[,paste0("beta", 0:10)]
Abparm <- array(bparm, c(nrow(bparm), ncol(bparm), nrow(newdat_))) # mcmcsamp, parameter, datasamp

newdat_[,"x0"] = 1
newdat_ <- as.matrix(newdat_)
newdat_ <- newdat_[,paste0("x",0:10)]

Anewdat_ <- array(newdat_, c(nrow(newdat_), ncol(newdat_), nrow(bparm))) # datasamp, parameter, mcmcsamp
Anewdat_ <- aperm(Anewdat_, c(3,2,1)) # mcmcsamp, parameter, datasamp

newpred_ <- apply((Anewdat_*Abparm), c(1,3), sum) # mcmcsamp, datasamp

newpred <- attr(dat_mat, "scaled:scale")['y'] * newpred_ + attr(dat_mat, "scaled:center")['y'] 

newpredMean <- apply(newpred, 2, mean)

dat2All <- X[, .(f1= x1+x2+x3+x4+x5, 
              f2= x6+x7+x8+x9+x10,
              y = y)]
dat2 <- dat2All[1:N]
newdat2 <- dat2All[(N+1):N2]
fit <- lm(y ~ . -e , dat)
fit2 <- lm(y ~ . , dat2)
print(mean((predict(fit, newdat)- newdat$y)^2))

## [1] 1.868142

print(mean((predict(fit2, newdat2)-newdat2$y)^2))

## [1] 1.347253

print(mean((newpredMean-newdat2$y)^2))

## [1] 1.59107

단 한번의 분석결과로 확신할 수 없지만 베이지안 분석은 두 극단(모두 계수를 자유롭게 두는 모형과 그룹 내의 모든 계수가 같다고 두는 모형)의 중간쯤 된다. 아니, 중간 쯤 되는 듯 보인다. (하지만 위의 코드를 랜덤 시드(set.seed)를 바꿔가며 돌려보면 할 때마다 매우 다른 결과가 나옴을 알 수 있다. 평균적 결과는 아래 부록을 확인하자.)

베이지안 추가 분석

계수의 사후 분포

점선은 1과 -1.4. 동그란 점은 실제 모수. 바이올린 플롯은 계수의 사후 분포(posterior distribution)

library(tidyr)
sdy = attr(dat_mat, "scaled:scale")['y']
meany = attr(dat_mat, "scaled:center")['y']

beta0hat <- meany + sdy*bparm[, 'beta0'] -
  apply(sdy*
  sweep(sweep(bparm[, paste0("beta", 1:10)],2,attr(dat_mat, "scaled:center")[paste0("x", 1:10)],"*"),
        2,
        attr(dat_mat, "scaled:scale")[paste0("x", 1:10)],"/"), 1, sum)

dfBparmOrigScale <- sweep(bparm[, paste0("beta", 1:10)]*sdy, 2, 
    attr(dat_mat, "scaled:scale")[paste0("x", 1:10)], "/") %>% as.data.frame 

dfBparmOrigScale$beta0 <- beta0hat

dfBparmOrigScale %>%
    gather(key='key', value='value', beta1:beta0) %>% 
#dfBparmOrigScale %>% #mutate(key=factor(key, levels=paste0("beta", 0:10))) %>%
  ggplot(aes(x=key, y=value)) + geom_violin(fill='gray80') + 
  geom_point(data=dfParamTrue, aes(x=key, y=value), size=2) +
  scale_x_discrete(limits=paste0("beta", 0:10)) + 
  geom_hline(yintercept=1.1, linetype='dotted') + 
  geom_hline(yintercept=-1.4, linetype='dotted') + 
  theme_light() 

계수의 사후 결합 분포

dfBparam <- bparm[, paste0("beta", 1:10)] %>% as.data.frame 

g2 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta2), alpha=0.01) + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g3 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta3), alpha=0.01) + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g4 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta4), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g5 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta5), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g6 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta6), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g7 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta7), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g8 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta8), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g9 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta9), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()
g10 <- ggplot(dfBparam) + geom_point(aes(x=beta1, y=beta10), alpha=0.01)  + coord_cartesian(xlim=range(bparm[,'beta1'])) + coord_equal() + theme_light()

library(patchwork)
(g2 + g3 + g4) /
  (g5 + g6 + g7) /
  (g8 + g9 + g10)

마무리

만약 계수를 모두 풀어지는 모형과 그룹내의 계수를 모두 같게 놓는 두 모형만이 존재한다면, 그 둘 사이에 무엇을 선택할지에 대해 고민할 수 밖에 없다. 이때 가장 좋은 방법은 각 집단 내의 모수가 평균과 차이가 나는 정도, 그리고 표본크기를 고려하여 가장 좋은 모형을 선택하는 것이다.

이에 대한 적절한 정보가 없다면 베이지안으로 분석하는 것이 속편할 수 있다.

부록

부록 A.

선형 모형을 위해 주어진 변수를 모두 표준화하고, 예측값에서 다시 원래 값을 얻는 과정은 다음과 같다.

\(y = \beta_0+ \beta_1 x_1 + \beta_2 x_2 + \cdots\) 에서 \(x_1, x_2, \cdots, x_{10}\) 을 모두 표준화를 하면 다음과 같다. 이때 계수가 \(\beta_1\) 에서 $\beta_1'$로 변화되었음을 유의한다.

\[\frac{y – \bar{y}}{\textrm{sd}(y)} = \beta_0′ + \beta_1’\frac{x_1-\bar{x_1}}{\textrm{sd}(x_1)} + \beta_2’\frac{x_2-\bar{x_2}}{\textrm{sd}(x_2)} + \cdots\]

이를 아래와 같이 잘 정리하면 \(\beta_1’\) 과 \(\beta_1\) 의 관계를 유추해낼 수 있다.

\[ y = \bar{y} + \textrm{sd}(y)\beta’_0 + \textrm{sd}(y)\beta_1’\frac{x_1-\bar{x_1}}{\textrm{sd}(x_1)} + \textrm{sd}(y)\beta_2’\frac{x_2-\bar{x_2}}{\textrm{sd}(x_2)} \]

\(\beta_0’\) 에서 \(\beta_0\) 을 계산해내는 식도 다음과 같이 유도할 수 있다.

\[ y =\left( \bar{y} + \textrm{sd}(y)\beta’_0 – \frac{\textrm{sd}(y)\beta’_1\bar{x_1}}{\textrm{sd}(x_1)} – \frac{\textrm{sd}(y)\beta’_2\bar{x_2}}{\textrm{sd}(x_2)} + \cdots \right) + \frac{\textrm{sd}(y)\beta_1′}{\textrm{sd}(x_1)}x_1 + \frac{\textrm{sd}(y)\beta_2′}{\textrm{sd}(x_2)}x_2 + \cdots \]

부록 B. 시뮬레이션 연구

func = function(f1,f2, e) {
  return(2 + f1 - 1.4*f2 + e)
  #return(1 + 2*sin(f1/2) + 3 * cos(f2/2))
  #return(1 + sin(f1/2) + cos(f2/2))
  #return(1 + 2 * sin(f1/2) * cos(f2/3))
  #return(1 + sin(f1) * cos(f2))
}

jags.model <- "


model {
  for (i in 1:N) {
    ymean[i] <- beta0 + beta1*x1[i] + beta2*x2[i] + beta3*x3[i] + beta4*x4[i] + beta5*x5[i] + 
            beta6*x6[i] + beta7*x7[i] + beta8*x8[i] + beta9*x9[i] + beta10*x10[i]
    y[i] ~ dnorm(ymean[i], tau)
  }

  # priors
  tau <- pow(sigma2, -2)
  sigma2 ~ dunif(0,100)
  beta0 ~ dnorm(0, 0.01)
  beta1 ~ dnorm(betaAmean, betaAtau)
  beta2 ~ dnorm(betaAmean, betaAtau)
  beta3 ~ dnorm(betaAmean, betaAtau)
  beta4 ~ dnorm(betaAmean, betaAtau)
  beta5 ~ dnorm(betaAmean, betaAtau)
  beta6 ~ dnorm(betaBmean, betaBtau)
  beta7 ~ dnorm(betaBmean, betaBtau)
  beta8 ~ dnorm(betaBmean, betaBtau)
  beta9 ~ dnorm(betaBmean, betaBtau)
  beta10 ~ dnorm(betaBmean, betaBtau)
  betaAmean ~ dunif(-100,100)
  betaAtau <- pow(betaAsigma2, -2)
  betaBmean ~ dunif(-100,100)
  betaBtau <- pow(betaBsigma2, -2)
  betaAsigma2 ~ dunif(0,100)
  betaBsigma2 ~ dunif(0,100)
}
"

cat(jags.model, file='currentModel.j')

samplesize = c(20,40,60,80,100,200,300)
nrep = 100
resultSim <- expand.grid(irep=1:nrep, samplesize=samplesize)

resultSim <- resultSim %>% mutate(mod_full=NA, mod_factor=NA, mod_bayes=NA,
                                  t_mod_full = NA, t_mod_factor = NA, t_mod_bayes = NA,
                                  m_mod_full = NA, m_mod_factor = NA, m_mod_bayes = NA)
# 1. MSE
# 2. time
# 3. memory

for (isim in 1:nrow(resultSim)) {
  N <- resultSim[isim, "samplesize"]
  rm(dat,newdat, dat2, newdat2)

  # Training set
  ncolMatX <- 11
  sigma.cor = 0.1
  sd.paramTrue = 0.2

  sigma <- matrix(
    sigma.cor,
    ncolMatX, ncolMatX)
  diag(sigma) = 1

  matX <- rmvnorm(N, mean=rep(0,ncolMatX), sigma=sigma)
  colnames(matX) <- c(paste0("x", 1:10), "e")
  X <- data.table(matX)

  paramTrue <- c(3, rnorm(5,1,sd.paramTrue),
                 rnorm(5,-1.4,sd.paramTrue))
  names(paramTrue) <- paste0("beta", 0:10)

  lv <- X[, .(f1 = paramTrue['beta1']*x1 + 
                paramTrue['beta2']*x2 + 
                paramTrue['beta3']*x3 + 
                paramTrue['beta4']*x4 + 
                paramTrue['beta5']*x5,
              f2 = paramTrue['beta6']*x6 + 
                paramTrue['beta7']*x7 + 
                paramTrue['beta8']*x8 + 
                paramTrue['beta9']*x9 + 
                paramTrue['beta10']*x10)]

  y <- func(lv$f1, lv$f2, X$e)

  dfParamTrue <- data.frame(key=factor(names(paramTrue), 
                                       levels= paste0("beta", 0:10)),
                            value=paramTrue)

  dat <- data.table(X,y)
  dat2 <- X[, .(f1= x1+x2+x3+x4+x5, 
                f2= x6+x7+x8+x9+x10,
                y = y)]
  ###

  ### Training
  t1 <- Sys.time()
  # full linear model
  fit <- lm(y ~ . - e, dat)
  t2 <- Sys.time()
  # factor linear model
  fit2 <- lm(y ~ ., dat2)
  t3 <- Sys.time()
  # bayesian linear model
  dat_ <- (dat_mat <- 
             dat %>% select(-e) %>% scale) %>% data.frame

  N = nrow(dat_)

  myj <- jags.model('currentModel.j',
                    data = c(as.list(dat_), 'N'=N),
                    n.chains = 3)

  update(myj, n.iter=2000)

  n.thin=5

  mcmcsamp <- 
    coda.samples(myj, 
                 c('beta0', 'beta1', 'beta2', 'beta3', 'beta4', 'beta5',
                   'beta6', 'beta7', 'beta8', 'beta9', 'beta10',
                   'betaAmean', 'betaBmean', 'betaAsigma2', 'betaBsigma2', 'sigma2'),
                 1000*n.thin, thin=n.thin)  # if there no betaAmean, warning will be issue

  t4 <- Sys.time()

  rm(dat, dat2)

  resultSim[isim, "t_mod_full"] = t2 - t1
  resultSim[isim, "t_mod_factor"] = t3 - t2
  resultSim[isim, "t_mod_bayes"] = t4 - t3

  resultSim[isim, "m_mod_full"] = as.numeric(object.size(fit))/1024/1024 # MB
  resultSim[isim, "m_mod_factor"] = as.numeric(object.size(fit2))/1024/1024 # MB
  resultSim[isim, "m_mod_bayes"] = as.numeric(object.size(mcmcsamp))/1024/1024 # MB

  #########################
  ## NEW DATA

  N <- 500

  ncolMatX <- 11; sigma.cor = 0.1

  sigma <- matrix(
    sigma.cor,
    ncolMatX, ncolMatX)
  diag(sigma) = 1

  matX <- rmvnorm(N, mean=rep(0,ncolMatX), sigma=sigma)
  colnames(matX) <- c(paste0("x", 1:10), "e")
  X <- data.table(matX)

  lv <- X[, .(f1 = paramTrue['beta1']*x1 + 
                paramTrue['beta2']*x2 + 
                paramTrue['beta3']*x3 + 
                paramTrue['beta4']*x4 + 
                paramTrue['beta5']*x5,
              f2 = paramTrue['beta6']*x6 + 
                paramTrue['beta7']*x7 + 
                paramTrue['beta8']*x8 + 
                paramTrue['beta9']*x9 + 
                paramTrue['beta10']*x10)]

  y <- func(lv$f1, lv$f2, X$e)

  newdat <- data.table(X,y)
  newdat2 <- X[, .(f1= x1+x2+x3+x4+x5, 
                   f2= x6+x7+x8+x9+x10,
                   y = y)]

  ####

  ## Bayesian model prediction

  newdat_ <- sweep(newdat %>% select(-y), 2, attr(dat_mat, "scaled:center"), "-")
  newdat_ <- sweep(newdat_ , 2, attr(dat_mat, "scaled:scale"), "/")

  #mcmcsamp
  #str(mcmcsamp)
  bparm <- rbind(mcmcsamp[[1]], mcmcsamp[[2]], mcmcsamp[[3]])[,paste0("beta", 0:10)]
  Abparm <- array(bparm, c(nrow(bparm), ncol(bparm), nrow(newdat_))) # mcmcsamp, parameter, datasamp

  newdat_[,"x0"] = 1
  newdat_ <- as.matrix(newdat_)
  newdat_ <- newdat_[,paste0("x",0:10)]

  Anewdat_ <- array(newdat_, c(nrow(newdat_), ncol(newdat_), nrow(bparm))) # datasamp, parameter, mcmcsamp
  Anewdat_ <- aperm(Anewdat_, c(3,2,1)) # mcmcsamp, parameter, datasamp

  newpred_ <- apply((Anewdat_*Abparm), c(1,3), sum) # mcmcsamp, datasamp

  newpred <- attr(dat_mat, "scaled:scale")['y'] * newpred_ + attr(dat_mat, "scaled:center")['y'] 

  newpredMean <- apply(newpred, 2, mean)
  ####

  resultSim[isim, "mod_full"] <- mean((predict(fit, newdat)- newdat$y)^2)
  resultSim[isim, "mod_factor"] <- mean((predict(fit2, newdat2)-newdat2$y)^2)
  resultSim[isim, "mod_bayes"] <- mean((newpredMean-newdat2$y)^2)

  if (isim %% 10 == 0) print(isim)
}

결과는 다음과 같다. (두 번째 그래프는 베이지안 모형과 다른 모형과의 비교를 위한 그래프이고, 세 번째 그래프는 각 모형의 결과를 따로 그린 그림. 결과적은 셋은 모두 같은 결과를 나타낸다.)

• mod_full: 전체 선형
• mod_factor: 요인 선형
• mod_bayes: 베이지언 선형

언제나 베이지안이 최고는 아니지만, 평균적으로 베이지안이 최고다!

물론 선형 모형이라는 한계가 있지만, 트리모형에도 BART(Bayesian Additive Regression Trees)라는 베이지안 모형이 존재한다! (표본 크기가 작거나, 컴퓨터에 남아 도는 core가 많을 때 한 번 시도를…)

• BART: Bayesian Additive Regression Trees, https://arxiv.org/abs/0806.3286
• https://blog.zakjost.com/post/bart/