Applied Bayesian Method
Expert information - prior information
From doing research and running the regression on the actual dataset, we have
- age : a unit increase in age, increases the insurance charges by 250 dollars.
- sex : If the client is a female, the charges will be 2000 AUD less than that of a male on the average. This is a very strong expert knowledge.
- bmi : There is no expert knowledge on the body mass index.
- children: Having one more child can cost a client 700 AUD more. This is a weak expert information.
- smoker: The clients who smoke can cost more than those who do not. This is a weak expert knowledge.
- region: People living in the southest have to pay more money on medical insurance than those who live in southwest and northwest. This is a strong knowledge.
The prediction
we use the model to predict with these three queries:
The query predictions
graphics.off()
rm(list=ls())
setwd("C:/Users/mlcl.local/Downloads")library(magrittr)
library(rjags)
library(coda)
library(runjags)
library(tufte)data <- read.csv("insurance_data.csv",header = TRUE )
head(data)## age sex bmi children smoker southwest southeast northeast northwest
## 1 19 1 27.900 0 1 1 0 0 0
## 2 18 0 33.770 1 0 0 1 0 0
## 3 28 0 33.000 3 0 0 1 0 0
## 4 33 0 22.705 0 0 0 0 0 1
## 5 32 0 28.880 0 0 0 0 0 1
## 6 31 1 25.740 0 0 0 1 0 0
## charges
## 1 16884.924
## 2 1725.552
## 3 4449.462
## 4 21984.471
## 5 3866.855
## 6 3756.622Charges distribution:
hist(data$charges, main="Histogram of Medical Charges",
xlab="Charges",
col="darkmagenta",
freq=FALSE, breaks = 20)
## Multiple predictors:
yName <- "charges"
xName <- c("age","sex","bmi","children","smoker","southwest","southeast",
"northeast")
fileNameRoot1 <- "charges_reg"
numSavedSteps <- 3000
thinSteps <- 4The model
Below is the JAGS model specification:
modelString = "
# Standardize the data:
data {
ym <- mean(y)
ysd <- sd(y)
for ( i in 1:Ntotal ) {
zy[i] <- ( y[i] - ym ) / ysd
}
for ( j in 1:Nx ) {
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for ( i in 1:Ntotal ) {
zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
}
}
# Specify the priors for original beta parameters
# Prior locations to reflect the expert information
mu0 <- ym # Set to overall mean a priori based on the interpretation of constant term in regression
mu[1] <- 250 # age
mu[2] <- -2000 # sex=female
mu[3] <- 350 # bmi
mu[4] <- 700 # number of children
mu[5] <- 20000 # smoker=1
mu[6] <- -1000 # southwest
mu[7] <- 2000 # southeast
mu[8] <- 200 # northeast
# Prior variances to reflect the expert information
Var0 <- (ysd^2)*10
Var[1] <- (xsd[1]^2)*10 # age
Var[2] <- (xsd[2]^2)*10 # sex=female
Var[3] <- (xsd[3]^2)*10 # bmi
Var[4] <- (xsd[4]^2)*10 # number of children
Var[5] <- (xsd[5]^2)*10 # smoker=1
Var[6] <- (xsd[6]^2)*10 # southwest
Var[7] <- (xsd[7]^2)*10 # southeast
Var[8] <- (xsd[8]^2)*10 # northeast
# Compute corresponding prior means and variances for the standardised parameters
muZ[1:Nx] <- mu[1:Nx] * xsd[1:Nx] / ysd
muZ0 <- (mu0 + sum( mu[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd - ym) / ysd
# Compute corresponding prior variances and variances for the standardised parameters
VarZ[1:Nx] <- Var[1:Nx] * ( xsd[1:Nx]/ ysd )^2
VarZ0 <- Var0 / (ysd^2)
}
# Specify the model for standardized data:
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( sum(zbeta[1:Nx]*xm[1:Nx]/xsd[1:Nx] - ym/ysd), 1/VarZ0 )
for ( j in 1:Nx ) {
zbeta[j] ~ dnorm( muZ[j] , 1/VarZ[j] )
}
zsigma ~ dunif( 1.0E-5 , 1.0E+1 )
nu ~ dexp(1/30.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
# Compute predictions at every step of the MCMC
pred1 <- beta0 + beta[1] * xPred1[1] + beta[2] * xPred1[2] + beta[3] * xPred1[3] +
beta[4] * xPred1[4] + beta[5] * xPred1[5] + beta[6]*xPred1[6] +
beta[7]*xPred1[7] + beta[8]*xPred1[8]
pred2 <- beta0 + beta[1] * xPred2[1] + beta[2] * xPred2[2] + beta[3] * xPred2[3] +
beta[4] * xPred2[4] + beta[5] * xPred2[5] + beta[6]*xPred2[6] +
beta[7]*xPred2[7] + beta[8]*xPred2[8]
pred3 <- beta0 + beta[1] * xPred3[1] + beta[2] * xPred3[2] + beta[3] * xPred3[3] +
beta[4] * xPred3[4] + beta[5] * xPred3[5] + beta[6]*xPred3[6] +
beta[7]*xPred3[7] + beta[8]*xPred3[8]
}
"Chain Generation
names <- c("beta[1]","beta[2]","beta[3]","beta[4]","beta[5]","beta[6]",
"beta[7]","beta[8]")
xpred1 <- c(18,1,35,0,1,0,0,1)
names(xpred1) <- names
xpred2 <- c(40,0,25,2,1,1,0,0)
names(xpred2) <- names
xpred3 <- c(30,1,20,1,0,0,1,0)
names(xpred3) <- namessource("C:/Users/mlcl.local/Downloads/fpfunctions - Copy.R")##
## *********************************************************************
## Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
## A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
## *********************************************************************source("C:/Users/mlcl.local/Downloads/DBDA2E-utilities.R")##
## *********************************************************************
## Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
## A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
## *********************************************************************mcmccoda1 <- genMCMC(data=data, xName=xName, yName=yName, numSavedSteps=numSavedSteps,
thinSteps=thinSteps, saveName=fileNameRoot1, xPred1 = xpred1,
xPred2=xpred2, xPred3=xpred3)##
## CORRELATION MATRIX OF PREDICTORS:
## age sex bmi children smoker southwest southeast
## age 1.000 0.021 0.109 0.042 -0.025 0.010 -0.012
## sex 0.021 1.000 -0.046 -0.017 -0.076 0.004 -0.017
## bmi 0.109 -0.046 1.000 0.013 0.004 -0.006 0.270
## children 0.042 -0.017 0.013 1.000 0.008 0.022 -0.023
## smoker -0.025 -0.076 0.004 0.008 1.000 -0.037 0.068
## southwest 0.010 0.004 -0.006 0.022 -0.037 1.000 -0.346
## southeast -0.012 -0.017 0.270 -0.023 0.068 -0.346 1.000
## northeast 0.002 0.002 -0.138 -0.023 0.003 -0.320 -0.346
## northeast
## age 0.002
## sex 0.002
## bmi -0.138
## children -0.023
## smoker 0.003
## southwest -0.320
## southeast -0.346
## northeast 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 4.3.0 on Sat Oct 06 13:26:02 2018
## 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 information:
## Observed stochastic nodes: 1338
## Unobserved stochastic nodes: 11
## Total graph size: 29581
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 1000
## -------------------------------------------------| 1000
## ************************************************** 100%
## . . . . . . . . . . . Updating 4000
## -------------------------------------------------| 4000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Finished running the simulation
