In this chapter, we are concerned with situations in which the value to be predicted is on a metric scale, and there are several predictors, each of which is also on a metric scale.
y ∼ normal(μ, σ ) and μ = β0+β1x1+β2x2.
The data block of the JAGS code then standardizes the predictors by looping through the columns of the x matrix, as follows:
data {
ym <- mean(y)
ysd <- sd(y)
for ( i in 1:Ntotal ) { # Ntotal is the number of data rows
zy[i] <- ( y[i] - ym ) / ysd
}
for ( j in 1:Nx ) {# Nx is the number of x predictors
xm[j] <- mean(x[,j]) # x is a matrix, each column a predictor
xsd[j] <- sd(x[,j])
for ( i in 1:Ntotal ) {
zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
}
}
}
The standardized parameters are then transformed to the original scale by generalizing Equation 17.2 (p. 485) to multiple predictors:
The JAGS model specification looks like this:
model {
for ( i in 1:Ntotal ) {
zy[i] ̃ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigmaˆ2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 )
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/2ˆ2 )
}
zsigma ̃ dunif( 1.0E-5 , 1.0E+1 )
nu <- nuMinusOne+1
nuMinusOne ̃ dexp(1/29.0)
# Transform to original scale:
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
}
# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업
cur_dir = getwd()
setwd(sprintf("%s/%s", cur_dir, 'data'))
# Example for Jags-Ymet-XmetMulti-Mrobust.R
#-------------------------------------------------------------------------------
# Optional generic preliminaries:
#graphics.off() # This closes all of R's graphics windows.
#rm(list=ls()) # Careful! This clears all of R's memory!
#.............................................................................
# # Two predictors:
myData = read.csv( file="Guber1999data.csv" )
myData
State | Spend | StuTeaRat | Salary | PrcntTake | SATV | SATM | SATT | |
---|---|---|---|---|---|---|---|---|
1 | Alabama | 4.405 | 17.2 | 31.144 | 8 | 491 | 538 | 1029 |
2 | Alaska | 8.963 | 17.6 | 47.951 | 47 | 445 | 489 | 934 |
3 | Arizona | 4.778 | 19.3 | 32.175 | 27 | 448 | 496 | 944 |
4 | Arkansas | 4.459 | 17.1 | 28.934 | 6 | 482 | 523 | 1005 |
5 | California | 4.992 | 24 | 41.078 | 45 | 417 | 485 | 902 |
6 | Colorado | 5.443 | 18.4 | 34.571 | 29 | 462 | 518 | 980 |
7 | Connecticut | 8.817 | 14.4 | 50.045 | 81 | 431 | 477 | 908 |
8 | Delaware | 7.03 | 16.6 | 39.076 | 68 | 429 | 468 | 897 |
9 | Florida | 5.718 | 19.1 | 32.588 | 48 | 420 | 469 | 889 |
10 | Georgia | 5.193 | 16.3 | 32.291 | 65 | 406 | 448 | 854 |
11 | Hawaii | 6.078 | 17.9 | 38.518 | 57 | 407 | 482 | 889 |
12 | Idaho | 4.21 | 19.1 | 29.783 | 15 | 468 | 511 | 979 |
13 | Illinois | 6.136 | 17.3 | 39.431 | 13 | 488 | 560 | 1048 |
14 | Indiana | 5.826 | 17.5 | 36.785 | 58 | 415 | 467 | 882 |
15 | Iowa | 5.483 | 15.8 | 31.511 | 5 | 516 | 583 | 1099 |
16 | Kansas | 5.817 | 15.1 | 34.652 | 9 | 503 | 557 | 1060 |
17 | Kentucky | 5.217 | 17 | 32.257 | 11 | 477 | 522 | 999 |
18 | Louisiana | 4.761 | 16.8 | 26.461 | 9 | 486 | 535 | 1021 |
19 | Maine | 6.428 | 13.8 | 31.972 | 68 | 427 | 469 | 896 |
20 | Maryland | 7.245 | 17 | 40.661 | 64 | 430 | 479 | 909 |
21 | Massachusetts | 7.287 | 14.8 | 40.795 | 80 | 430 | 477 | 907 |
22 | Michigan | 6.994 | 20.1 | 41.895 | 11 | 484 | 549 | 1033 |
23 | Minnesota | 6 | 17.5 | 35.948 | 9 | 506 | 579 | 1085 |
24 | Mississippi | 4.08 | 17.5 | 26.818 | 4 | 496 | 540 | 1036 |
25 | Missouri | 5.383 | 15.5 | 31.189 | 9 | 495 | 550 | 1045 |
26 | Montana | 5.692 | 16.3 | 28.785 | 21 | 473 | 536 | 1009 |
27 | Nebraska | 5.935 | 14.5 | 30.922 | 9 | 494 | 556 | 1050 |
28 | Nevada | 5.16 | 18.7 | 34.836 | 30 | 434 | 483 | 917 |
29 | New Hampshire | 5.859 | 15.6 | 34.72 | 70 | 444 | 491 | 935 |
30 | New Jersey | 9.774 | 13.8 | 46.087 | 70 | 420 | 478 | 898 |
31 | New Mexico | 4.586 | 17.2 | 28.493 | 11 | 485 | 530 | 1015 |
32 | New York | 9.623 | 15.2 | 47.612 | 74 | 419 | 473 | 892 |
33 | North Carolina | 5.077 | 16.2 | 30.793 | 60 | 411 | 454 | 865 |
34 | North Dakota | 4.775 | 15.3 | 26.327 | 5 | 515 | 592 | 1107 |
35 | Ohio | 6.162 | 16.6 | 36.802 | 23 | 460 | 515 | 975 |
36 | Oklahoma | 4.845 | 15.5 | 28.172 | 9 | 491 | 536 | 1027 |
37 | Oregon | 6.436 | 19.9 | 38.555 | 51 | 448 | 499 | 947 |
38 | Pennsylvania | 7.109 | 17.1 | 44.51 | 70 | 419 | 461 | 880 |
39 | Rhode Island | 7.469 | 14.7 | 40.729 | 70 | 425 | 463 | 888 |
40 | South Carolina | 4.797 | 16.4 | 30.279 | 58 | 401 | 443 | 844 |
41 | South Dakota | 4.775 | 14.4 | 25.994 | 5 | 505 | 563 | 1068 |
42 | Tennessee | 4.388 | 18.6 | 32.477 | 12 | 497 | 543 | 1040 |
43 | Texas | 5.222 | 15.7 | 31.223 | 47 | 419 | 474 | 893 |
44 | Utah | 3.656 | 24.3 | 29.082 | 4 | 513 | 563 | 1076 |
45 | Vermont | 6.75 | 13.8 | 35.406 | 68 | 429 | 472 | 901 |
46 | Virginia | 5.327 | 14.6 | 33.987 | 65 | 428 | 468 | 896 |
47 | Washington | 5.906 | 20.2 | 36.151 | 48 | 443 | 494 | 937 |
48 | West Virginia | 6.107 | 14.8 | 31.944 | 17 | 448 | 484 | 932 |
49 | Wisconsin | 6.93 | 15.9 | 37.746 | 9 | 501 | 572 | 1073 |
50 | Wyoming | 6.16 | 14.9 | 31.285 | 10 | 476 | 525 | 1001 |
summary(myData)
State Spend StuTeaRat Salary Alabama : 1 Min. :3.656 Min. :13.80 Min. :25.99 Alaska : 1 1st Qu.:4.882 1st Qu.:15.22 1st Qu.:30.98 Arizona : 1 Median :5.768 Median :16.60 Median :33.29 Arkansas : 1 Mean :5.905 Mean :16.86 Mean :34.83 California: 1 3rd Qu.:6.434 3rd Qu.:17.57 3rd Qu.:38.55 Colorado : 1 Max. :9.774 Max. :24.30 Max. :50.05 (Other) :44 PrcntTake SATV SATM SATT Min. : 4.00 Min. :401.0 Min. :443.0 Min. : 844.0 1st Qu.: 9.00 1st Qu.:427.2 1st Qu.:474.8 1st Qu.: 897.2 Median :28.00 Median :448.0 Median :497.5 Median : 945.5 Mean :35.24 Mean :457.1 Mean :508.8 Mean : 965.9 3rd Qu.:63.00 3rd Qu.:490.2 3rd Qu.:539.5 3rd Qu.:1032.0 Max. :81.00 Max. :516.0 Max. :592.0 Max. :1107.0
yName = "SATT" ; xName = c("Spend","PrcntTake")
fileNameRoot = "Guber1999data-Jags-"
numSavedSteps=15000 ; thinSteps=5
graphFileType = "png"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Jags-Ymet-XmetMulti-Mrobust.R")
********************************************************************* Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier. *********************************************************************
Loading required package: coda
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
#startTime = proc.time()
mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
saveName=fileNameRoot )
CORRELATION MATRIX OF PREDICTORS: Spend PrcntTake Spend 1.000 0.593 PrcntTake 0.593 1.000 Calling 3 simulations using the parallel method... Following the progress of chain 1 (the program will wait for all chains to finish before continuing): Welcome to JAGS 3.4.0 on Wed Sep 16 20:42:43 2015 JAGS is free software and comes with ABSOLUTELY NO WARRANTY Loading module: basemod: ok Loading module: bugs: ok . . Reading data file data.txt . Compiling data graph Resolving undeclared variables Allocating nodes Initializing Reading data back into data table Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 535 . Reading parameter file inits1.txt . Initializing model . Adapting 500 -------------------------------------------------| 500 ++++++++++++++++++++++++++++++++++++++++++++++++++ 100% Adaptation successful . Updating 1000 -------------------------------------------------| 1000 ************************************************** 100% . . . . . . . . Updating 25000 -------------------------------------------------| 25000 ************************************************** 100% . . . . Updating 0 . Deleting model . All chains have finished Simulation complete. Reading coda files... Coda files loaded successfully Finished running the simulation
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in parameterNames ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
# 창으로 뜨는 것들을 없앤다.
graphics.off()
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda ,
saveName=fileNameRoot )
show(summaryInfo)
Mean Median Mode ESS HDImass HDIlow CHAIN 2.000000000 2.000000000 1.997413e+00 1.5 0.95 1.0000000 beta0 991.492221267 991.336500000 9.899081e+02 13808.0 0.95 948.3840000 beta[1] 12.832703606 12.858900000 1.273112e+01 13153.9 0.95 4.5361400 beta[2] -2.878617251 -2.879090000 -2.881414e+00 13227.0 0.95 -3.3193300 sigma 31.530577473 31.329000000 3.106441e+01 13692.3 0.95 23.9116000 zbeta0 -0.001200131 -0.001465985 -9.973416e-05 15000.0 0.95 -0.1229360 zbeta[1] 0.233739136 0.234216000 2.318889e-01 13153.9 0.95 0.0826227 zbeta[2] -1.029646875 -1.029815000 -1.030645e+00 13227.0 0.95 -1.1872800 zsigma 0.421415965 0.418722000 4.151860e-01 13692.3 0.95 0.3195860 nu 33.397469008 24.824800000 1.103249e+01 10801.1 0.95 1.8336700 log10(nu) 1.375887230 1.394885757 1.485261e+00 12329.1 0.95 0.6664712 HDIhigh CompVal PcntGtCompVal ROPElow ROPEhigh PcntLtROPE CHAIN 3.000000 NA NA NA NA NA beta0 1035.490000 NA NA NA NA NA beta[1] 21.497200 NA NA NA NA NA beta[2] -2.442310 NA NA NA NA NA sigma 39.230700 NA NA NA NA NA zbeta0 0.121780 NA NA NA NA NA zbeta[1] 0.391557 NA NA NA NA NA zbeta[2] -0.873587 NA NA NA NA NA zsigma 0.524331 NA NA NA NA NA nu 90.746200 NA NA NA NA NA log10(nu) 2.072673 NA NA NA NA NA PcntInROPE PcntGtROPE CHAIN NA NA beta0 NA NA beta[1] NA NA beta[2] NA NA sigma NA NA zbeta0 NA NA zbeta[1] NA NA zbeta[2] NA NA zsigma NA NA nu NA NA log10(nu) NA NA
# Display posterior information:
plotMCMC( mcmcCoda , data=myData , xName=xName , yName=yName ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
# 창으로 뜨는 것들을 없앤다.
graphics.off()
can compute a surrogate. At each step in the MCMC chain, a credible value of R^2 is computed as
where
is the standardized regression coefficient for the jth predictor at that step in the MCMC chain,
and
is the correlation of the predicted values,
with the jth predictor values,
we conceptualize the interaction term of Equation 18.2 as an additional additive predictor, like this:
# Read in data:
myData = read.csv( file="Guber1999data.csv" )
# Append the new interaction variable:
myData = cbind( myData , SpendXPrcnt = myData[,"Spend"] * myData[,"PrcntTake"] )
# Specify names of data columns to use in the analysis:
yName = "SATT" ; xName = c("Spend","PrcntTake","SpendXPrcnt")
#### 수정된 버전 Jags-Ymet-XmetMulti-Mrobust-Example.R
# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업
cur_dir = getwd()
setwd(sprintf("%s/%s", cur_dir, 'data'))
#.............................................................................
# # Two predictors:
myData = read.csv( file="Guber1999data.csv" )
# Append the new interaction variable:
myData = cbind( myData , SpendXPrcnt = myData[,"Spend"] * myData[,"PrcntTake"] )
# Specify names of data columns to use in the analysis:
yName = "SATT" ; xName = c("Spend","PrcntTake","SpendXPrcnt")
myData
State | Spend | StuTeaRat | Salary | PrcntTake | SATV | SATM | SATT | SpendXPrcnt | |
---|---|---|---|---|---|---|---|---|---|
1 | Alabama | 4.405 | 17.2 | 31.144 | 8 | 491 | 538 | 1029 | 35.24 |
2 | Alaska | 8.963 | 17.6 | 47.951 | 47 | 445 | 489 | 934 | 421.261 |
3 | Arizona | 4.778 | 19.3 | 32.175 | 27 | 448 | 496 | 944 | 129.006 |
4 | Arkansas | 4.459 | 17.1 | 28.934 | 6 | 482 | 523 | 1005 | 26.754 |
5 | California | 4.992 | 24 | 41.078 | 45 | 417 | 485 | 902 | 224.64 |
6 | Colorado | 5.443 | 18.4 | 34.571 | 29 | 462 | 518 | 980 | 157.847 |
7 | Connecticut | 8.817 | 14.4 | 50.045 | 81 | 431 | 477 | 908 | 714.177 |
8 | Delaware | 7.03 | 16.6 | 39.076 | 68 | 429 | 468 | 897 | 478.04 |
9 | Florida | 5.718 | 19.1 | 32.588 | 48 | 420 | 469 | 889 | 274.464 |
10 | Georgia | 5.193 | 16.3 | 32.291 | 65 | 406 | 448 | 854 | 337.545 |
11 | Hawaii | 6.078 | 17.9 | 38.518 | 57 | 407 | 482 | 889 | 346.446 |
12 | Idaho | 4.21 | 19.1 | 29.783 | 15 | 468 | 511 | 979 | 63.15 |
13 | Illinois | 6.136 | 17.3 | 39.431 | 13 | 488 | 560 | 1048 | 79.768 |
14 | Indiana | 5.826 | 17.5 | 36.785 | 58 | 415 | 467 | 882 | 337.908 |
15 | Iowa | 5.483 | 15.8 | 31.511 | 5 | 516 | 583 | 1099 | 27.415 |
16 | Kansas | 5.817 | 15.1 | 34.652 | 9 | 503 | 557 | 1060 | 52.353 |
17 | Kentucky | 5.217 | 17 | 32.257 | 11 | 477 | 522 | 999 | 57.387 |
18 | Louisiana | 4.761 | 16.8 | 26.461 | 9 | 486 | 535 | 1021 | 42.849 |
19 | Maine | 6.428 | 13.8 | 31.972 | 68 | 427 | 469 | 896 | 437.104 |
20 | Maryland | 7.245 | 17 | 40.661 | 64 | 430 | 479 | 909 | 463.68 |
21 | Massachusetts | 7.287 | 14.8 | 40.795 | 80 | 430 | 477 | 907 | 582.96 |
22 | Michigan | 6.994 | 20.1 | 41.895 | 11 | 484 | 549 | 1033 | 76.934 |
23 | Minnesota | 6 | 17.5 | 35.948 | 9 | 506 | 579 | 1085 | 54 |
24 | Mississippi | 4.08 | 17.5 | 26.818 | 4 | 496 | 540 | 1036 | 16.32 |
25 | Missouri | 5.383 | 15.5 | 31.189 | 9 | 495 | 550 | 1045 | 48.447 |
26 | Montana | 5.692 | 16.3 | 28.785 | 21 | 473 | 536 | 1009 | 119.532 |
27 | Nebraska | 5.935 | 14.5 | 30.922 | 9 | 494 | 556 | 1050 | 53.415 |
28 | Nevada | 5.16 | 18.7 | 34.836 | 30 | 434 | 483 | 917 | 154.8 |
29 | New Hampshire | 5.859 | 15.6 | 34.72 | 70 | 444 | 491 | 935 | 410.13 |
30 | New Jersey | 9.774 | 13.8 | 46.087 | 70 | 420 | 478 | 898 | 684.18 |
31 | New Mexico | 4.586 | 17.2 | 28.493 | 11 | 485 | 530 | 1015 | 50.446 |
32 | New York | 9.623 | 15.2 | 47.612 | 74 | 419 | 473 | 892 | 712.102 |
33 | North Carolina | 5.077 | 16.2 | 30.793 | 60 | 411 | 454 | 865 | 304.62 |
34 | North Dakota | 4.775 | 15.3 | 26.327 | 5 | 515 | 592 | 1107 | 23.875 |
35 | Ohio | 6.162 | 16.6 | 36.802 | 23 | 460 | 515 | 975 | 141.726 |
36 | Oklahoma | 4.845 | 15.5 | 28.172 | 9 | 491 | 536 | 1027 | 43.605 |
37 | Oregon | 6.436 | 19.9 | 38.555 | 51 | 448 | 499 | 947 | 328.236 |
38 | Pennsylvania | 7.109 | 17.1 | 44.51 | 70 | 419 | 461 | 880 | 497.63 |
39 | Rhode Island | 7.469 | 14.7 | 40.729 | 70 | 425 | 463 | 888 | 522.83 |
40 | South Carolina | 4.797 | 16.4 | 30.279 | 58 | 401 | 443 | 844 | 278.226 |
41 | South Dakota | 4.775 | 14.4 | 25.994 | 5 | 505 | 563 | 1068 | 23.875 |
42 | Tennessee | 4.388 | 18.6 | 32.477 | 12 | 497 | 543 | 1040 | 52.656 |
43 | Texas | 5.222 | 15.7 | 31.223 | 47 | 419 | 474 | 893 | 245.434 |
44 | Utah | 3.656 | 24.3 | 29.082 | 4 | 513 | 563 | 1076 | 14.624 |
45 | Vermont | 6.75 | 13.8 | 35.406 | 68 | 429 | 472 | 901 | 459 |
46 | Virginia | 5.327 | 14.6 | 33.987 | 65 | 428 | 468 | 896 | 346.255 |
47 | Washington | 5.906 | 20.2 | 36.151 | 48 | 443 | 494 | 937 | 283.488 |
48 | West Virginia | 6.107 | 14.8 | 31.944 | 17 | 448 | 484 | 932 | 103.819 |
49 | Wisconsin | 6.93 | 15.9 | 37.746 | 9 | 501 | 572 | 1073 | 62.37 |
50 | Wyoming | 6.16 | 14.9 | 31.285 | 10 | 476 | 525 | 1001 | 61.6 |
fileNameRoot = "Guber1999data-interaction-Jags-"
numSavedSteps=15000 ; thinSteps=5
graphFileType = "png"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Jags-Ymet-XmetMulti-Mrobust.R")
********************************************************************* Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier. *********************************************************************
Loading required package: coda
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
#startTime = proc.time()
mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
saveName=fileNameRoot )
CORRELATION MATRIX OF PREDICTORS: Spend PrcntTake SpendXPrcnt Spend 1.000 0.593 0.775 PrcntTake 0.593 1.000 0.951 SpendXPrcnt 0.775 0.951 1.000 Calling 3 simulations using the parallel method... Following the progress of chain 1 (the program will wait for all chains to finish before continuing): Welcome to JAGS 3.4.0 on Wed Oct 7 02:03:51 2015 JAGS is free software and comes with ABSOLUTELY NO WARRANTY Loading module: basemod: ok Loading module: bugs: ok . . Reading data file data.txt . Compiling data graph Resolving undeclared variables Allocating nodes Initializing Reading data back into data table Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 638 . Reading parameter file inits1.txt . Initializing model . Adapting 500 -------------------------------------------------| 500 ++++++++++++++++++++++++++++++++++++++++++++++++++ 100% Adaptation successful . Updating 1000 -------------------------------------------------| 1000 ************************************************** 100% . . . . . . . . Updating 25000 -------------------------------------------------| 25000 ************************************************** 100% . . . . Updating 0 . Deleting model All chains have finished Simulation complete. Reading coda files... Coda files loaded successfully Finished running the simulation
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in parameterNames ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
# 창으로 뜨는 것들을 없앤다.
graphics.off()
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda ,
saveName=fileNameRoot )
show(summaryInfo)
Mean Median Mode ESS HDImass HDIlow CHAIN 2.000000e+00 2.00000000 1.997413e+00 1.5 0.95 1.0000000 beta0 1.046328e+03 1046.55500000 1.045477e+03 1192.8 0.95 964.1770000 beta[1] 2.706973e+00 2.64031500 2.081181e+00 1215.5 0.95 -13.2015000 beta[2] -4.064498e+00 -4.06793500 -4.089410e+00 898.6 0.95 -5.6597000 beta[3] 2.037744e-01 0.20437850 2.117125e-01 813.3 0.95 -0.0561615 sigma 3.104436e+01 30.81855000 3.040303e+01 13947.4 0.95 23.9590000 zbeta0 -1.564266e-03 -0.00124368 3.764121e-03 15000.0 0.95 -0.1236310 zbeta[1] 4.930570e-02 0.04809150 3.790802e-02 1215.5 0.95 -0.2404560 zbeta[2] -1.453822e+00 -1.45505000 -1.462733e+00 898.6 0.95 -2.0244100 zbeta[3] 5.615321e-01 0.56319650 5.834098e-01 813.3 0.95 -0.1547620 zsigma 4.149175e-01 0.41189950 4.063463e-01 13947.4 0.95 0.3202200 nu 3.490538e+01 26.42730000 1.158489e+01 10637.7 0.95 2.1464200 log10(nu) 1.400149e+00 1.42205278 1.457488e+00 12121.0 0.95 0.7047432 HDIhigh CompVal PcntGtCompVal ROPElow ROPEhigh PcntLtROPE CHAIN 3.000000 NA NA NA NA NA beta0 1129.180000 NA NA NA NA NA beta[1] 17.740100 NA NA NA NA NA beta[2] -2.461580 NA NA NA NA NA beta[3] 0.468960 NA NA NA NA NA sigma 38.820800 NA NA NA NA NA zbeta0 0.120886 NA NA NA NA NA zbeta[1] 0.323124 NA NA NA NA NA zbeta[2] -0.880479 NA NA NA NA NA zbeta[3] 1.292290 NA NA NA NA NA zsigma 0.518852 NA NA NA NA NA nu 94.077300 NA NA NA NA NA log10(nu) 2.096834 NA NA NA NA NA PcntInROPE PcntGtROPE CHAIN NA NA beta0 NA NA beta[1] NA NA beta[2] NA NA beta[3] NA NA sigma NA NA zbeta0 NA NA zbeta[1] NA NA zbeta[2] NA NA zbeta[3] NA NA zsigma NA NA nu NA NA log10(nu) NA NA
# Display posterior information:
plotMCMC( mcmcCoda , data=myData , xName=xName , yName=yName ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
# 창으로 뜨는 것들을 없앤다.
graphics.off()
When the analysis is run, the first thing it does is display the correlations of the predictors:
To properly understand the credible slopes on the two predictors, we must consider the credible slopes on each predictor as a function of the value of the other predictor.
In summary,
With so many candidate predictors of noisy data, there may be some regression coefficients that are spuriously estimated to be non-zero.
coefficients.
To illustrate these ideas with a concrete example, consider again the SAT data of Figure 18.3, but now supplemented with 12 more randomly generated predictors.
vague normal priors.
The resulting posterior distribution is shown in Figure 18.11.
the same extent as Spend fell above zero. - This apparent relation of xRand10 with SAT score is spurious, an accident of random sampling.
Now we repeat the analysis using the hierarchical model of Figure 18.10, with νβ fixed at one (i.e., heavy tails), and with σβ given a gamma prior that has mode 1.0 and standard deviation 1.0 (i.e., broad for the standardized data).
The resulting posterior distribution is shown in Figure 18.12.
regression coefficients are near zero. - Indeed, the estimate of σβ (not displayed) has its posterior mode around 0.05, even though its prior mode was at 1.0. - The shrinkage also causes the estimate of the regression coefficient on Spend to shift toward zero. - but it has also suppressed what might be a real but small regression coefficient on Spend.
. You can imagine that the posterior distribution zero and long concave tails, like this:
Deciding which predictors to include is often called variable selection.
have relations small enough to be counted as insignificant.
This section introduces some basic ideas and methods of Bayesian variable selection, using some simple illustrative examples.
reference for the latest and greatest methods.
from an independent Bernoulli prior, such as
model {
for ( i in 1:Ntotal ) {
zy[i] ̃ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigmaˆ2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 )
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/2ˆ2 )
}
zsigma ̃ dunif( 1.0E-5 , 1.0E+1 )
nu <- nuMinusOne+1
nuMinusOne ̃ dexp(1/29.0)
# Transform to original scale:
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
}
# Nx : number of predictors.
model {
for ( i in 1:Ntotal ) {
zy[i] ̃ dt( zbeta0 + sum( delta[1:Nx] * zbeta[1:Nx] * zx[i,1:Nx] ) ,
1/zsigmaˆ2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 )
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/2ˆ2 )
delta[j] ̃ dbern( 0.5 )
}
zsigma ̃ dunif( 1.0E-5 , 1.0E+1 )
nu <- nuMinusOne+1
nuMinusOne ̃ dexp(1/29.0)
# Transform to original scale:
beta[1:Nx] <- ( delta[1:Nx] * zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( delta[1:Nx] * zbeta[1:Nx] * xm[1:Nx]
/ xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
}
As a first example of applying the variable-selection method, recall the SAT data of Figure 18.3 and the posterior distribution shown in Figure 18.5.
With two predictors, there are four possible models involving different subsets of predictors.
Figure 18.13 shows the results.
We will now see that the degree of vagueness of the prior on the regression coefficient can have an enormous influence on the inclusion probability, even though the degree of vagueness has little influence on the estimate of the regression coefficient itself.
Recall that the prior in the model was specified as a generic broad distribution on the standardized regression coefficients, like this:
model {
...
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 ) # SD=2
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/2ˆ2 ) # SD=2
delta[j] ̃ dbern( 0.5 )
}
...
}
We now re-run the analysis using different arbitrary degrees of vagueness on the priors for the standardized regression coefficients.
We will illustrate with SD=1, like this:
model {
...
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/1ˆ2 ) # SD=1
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/1ˆ2 ) # SD=1
delta[j] ̃ dbern( 0.5 )
}
...
}
and with SD=10, like this:
model {
...
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/10ˆ2 ) # SD=10
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/10ˆ2 ) # SD=10
delta[j] ̃ dbern( 0.5 )
}
...
}
Figure 18.14 shows the results.
values may be minimally affected.
If you have strong previous research that can inform the prior, then it should be used. But if previous knowledge is weak, then the uncertainty should be expressed in the prior.
include a higher-level distribution to express your prior uncertainty.
braces.
In the present application, we have uncertainty, and therefore, we place a broad prior on σβ .
The code below shows a few different options in commented lines.
always have some arbitrary bounds. - Therefore, the next option, - not commented out of the code below, - is a gamma distribution that has a mode at 1.0 - but is very broad with a standard deviation of 10.0
model {
...
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 )
for ( j in 1:Nx ) {
zbeta[j] ̃ dt( 0 , 1/sigmaBetaˆ2 , 1 ) # Notice sigmaBeta
delta[j] ̃ dbern( 0.5 )
}
zsigma ̃ dunif( 1.0E-5 , 1.0E+1 )
## Uncomment one of the following specifications for sigmaBeta:
# sigmaBeta <- 2.0
# sigmaBeta ̃ dunif( 1.0E-5 , 1.0E+2 )
sigmaBeta ̃ dgamma(1.1051,0.1051) # mode 1.0, sd 10.0
# sigmaBeta <- 1/sqrt(tauBeta) ; tauBeta ̃ dgamma(0.001,0.001)
...
}
When it is run, it produces results very similar to those in Figure 18.13
# 18.13 관련 코드
# Nx : number of predictors.
model {
for ( i in 1:Ntotal ) {
zy[i] ̃ dt( zbeta0 + sum( delta[1:Nx] * zbeta[1:Nx] * zx[i,1:Nx] ) ,
1/zsigmaˆ2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ̃ dnorm( 0 , 1/2ˆ2 )
for ( j in 1:Nx ) {
zbeta[j] ̃ dnorm( 0 , 1/2ˆ2 )
delta[j] ̃ dbern( 0.5 )
}
zsigma ̃ dunif( 1.0E-5 , 1.0E+1 )
nu <- nuMinusOne+1
nuMinusOne ̃ dexp(1/29.0)
# Transform to original scale:
beta[1:Nx] <- ( delta[1:Nx] * zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( delta[1:Nx] * zbeta[1:Nx] * xm[1:Nx]
/ xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
}
# Example for Jags-Ymet-XmetMulti-MrobustVarSelect.R
# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업
cur_dir = getwd()
setwd(sprintf("%s/%s", cur_dir, 'data'))
#-------------------------------------------------------------------------------
# Load data file and specity column names of x (predictors) and y (predicted):
myData = read.csv( file="Guber1999data.csv" )
# UNCOMMENT ONE OF THE FOLLOWING SECTIONS (In RStudio, select and ctrl-shift-C)
#.............................................................................
# Two predictors:
yName = "SATT" ; xName = c("Spend","PrcntTake")
fileNameRoot = "Guber1999data-Jags-VarSelect-"
numSavedSteps=15000 ; thinSteps=25
#.............................................................................
# # Two predictors with redundant predictor:
# PropNotTake = (100-myData[,"PrcntTake"])/100
# myData = cbind( myData , PropNotTake )
# yName = "SATT" ; xName = c("Spend","PrcntTake","PropNotTake")
# fileNameRoot = "Guber1999data-Jags-Redund-VarSelect-"
# numSavedSteps=15000 ; thinSteps=30
#.............................................................................
# # Two predictors with two redundant predictors:
# PropNotTake = (100-myData[,"PrcntTake"])/100
# Partic = myData[,"PrcntTake"]/10
# myData = cbind( myData , PropNotTake , Partic )
# yName = "SATT" ; xName = c("Spend","PrcntTake","PropNotTake","Partic")
# fileNameRoot = "Guber1999data-Jags-Redund2-VarSelect-"
# numSavedSteps=15000 ; thinSteps=15
#.............................................................................
# # Four predictors:
# yName = "SATT" ; xName = c("Spend","PrcntTake","StuTeaRat","Salary")
# fileNameRoot = "Guber1999data-Jags-4X-VarSelect-"
# numSavedSteps=15000 ; thinSteps=20
#.............................................................................
# # Append columns of random predictors:
# standardize = function(v){(v-mean(v))/sd(v)}
# Ny=nrow(myData)
# NxRand = 12
# set.seed(47405)
# for ( xIdx in 1:NxRand ) {
# xRand = standardize(rnorm(Ny))
# myData = cbind( myData , xRand )
# colnames(myData)[ncol(myData)] = paste0("xRand",xIdx)
# }
# yName = "SATT" ; xName = c("Spend","PrcntTake", paste0("xRand",1:NxRand) )
# fileNameRoot = "Guber1999data-Jags-RandX-VarSelect-"
# numSavedSteps=15000 ; thinSteps=5
#.............................................................................
# # Two predictors with interaction.
# # Append named column with interaction product:
# myData = cbind( myData , SpendXPrcnt=myData[,"Spend"]*myData[,"PrcntTake"] )
# yName = "SATT" ; xName = c("Spend","PrcntTake","SpendXPrcnt")
# fileNameRoot = "Guber1999data-Jags-Inter-VarSelect-"
# numSavedSteps=15000 ; thinSteps=25
#.............................................................................
graphFileType = "png"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Jags-Ymet-XmetMulti-MrobustVarSelect.R")
********************************************************************* Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier. *********************************************************************
Installing packages into ‘/home/moosung/R/x86_64-pc-linux-gnu-library/3.2’ (as ‘lib’ is unspecified)
Error in contrib.url(repos, type): trying to use CRAN without setting a mirror
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
#startTime = proc.time()
mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
saveName=fileNameRoot )
#stopTime = proc.time()
#duration = stopTime - startTime
#show(duration)
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in parameterNames ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
# 창으로 뜨는 것들을 없앤다.
graphics.off()
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda ,
saveName=fileNameRoot )
show(summaryInfo)
# Display posterior information:
plotMCMC( mcmcCoda , data=myData , xName=xName , yName=yName ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
# 창으로 뜨는 것들을 없앤다.
graphics.off()
First we consider the correlations of the candidate predictors:
one.
Should only one or the other be included, or both, or neither?
Figure 18.15 shows the results.
be reported and what should be concluded?
Unfortunately, there is no singular “correct” answer.
책의 설명을 더 읽어보자.
MrobustVarSelect.R. The code is meant primarily for pedagogical purposes, and it does not scale well for larger applications, for reasons explained presently.
To conclude this section regarding variable selection, - it is appropriate to recapitulate the considerations at the beginning of the section. - Variable selection is a reasonable approach only if it is genuinely plausible and meaningful that candidate predictors have zero relation to the predicted variable. - The results can be surprisingly sensitive to the seemingly innocuous choice of prior for the regression coefficients, and, of course, the prior for the inclusion probability. - Because of these limitations, hierarchical shrinkage priors may be a more meaningful approach.