Multilevel: Anorexia Data
Carolyn J. Anderson
10.14.2022
- Setup
- Examining the data
- Maximum Likelihood Estimation
- Bayesian estimation
Below is code for a number of different models. The results of some but not all models are run below. You can do these on your own.
There is code for obtaining estimate of the posterior distribution of the icc (intra-class correlation), as well as the posterior probabilities for specific hypotheses. I also put in a better way to extract posterior values of parameter from samples (compared to what I did first time I did this).
Setup
The packages that we’ll use are
library(lme4)
library(lmerTest)
library(rjags)
library(runjags)
library(coda)
library(texreg)Set working directory to where the data live and set up data for initial analyses:
setwd("D:/Dropbox/edps 590BAY/Lectures/8 Multilevel models")
ano.wide <- read.table("anorexia_wide.txt",header=TRUE)
ano.wide$change <- ano.wide$weight2 - ano.wide$weight1
(n <- nrow(ano.wide))## [1] 72head(ano.wide)## Rx person weight1 weight2 change
## 1 1 1 80.5 82.2 1.7
## 2 1 2 84.9 85.6 0.7
## 3 1 3 81.5 81.4 -0.1
## 4 1 4 82.6 81.9 -0.7
## 5 1 5 79.9 76.4 -3.5
## 6 1 6 88.7 103.6 14.9Examining the data
Let’s do some graphics to look at the data. First we’ll just look at histograms of before treatment, after treatment, and change. We’ll add a 2nd graph of change but looking at is a bit differently. For the 2nd graph of change we need to re-format the data, and create treatment dummy codes
ano <- read.table("anorexia_long.txt",header=TRUE)
ano <- ano[order(ano$girl),]
ano$Rx1 <- ifelse(ano$Rx==1,1,0)
ano$Rx2 <- ifelse(ano$Rx==2,1,0)
ano$Rx3 <- ifelse(ano$Rx==3,1,0)and now for the graphs
par(mfrow=c(2,2))
hist(ano.wide$weight1, main='Weight Before Treatment',col="darkseagreen1")
hist(ano.wide$weight2, main='Weight After Treatment',col="darkseagreen1")
hist(ano.wide$change, main='Change in Weight',col="darkseagreen1")
plot(ano$time,ano$weight,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="Change for Each Girl"
)
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time,g$weight,col=i,lwd=2)
}
This 4th plot show different intercepts and different slopes. To better see what’s going on, we’ll compute means for each treatment before and after treatment
# means for each treatment
byRx <- aggregate(ano$weight,list(ano$Rx),FUN="mean")
# mean for treatment x time
rx_by_time <- aggregate(ano$weight,list(ano$Rx,ano$time),FUN="mean")
names(rx_by_time) <-c("rx","time","weight")
rx1 <- subset(rx_by_time,rx==1) # Cognitive
rx2 <- subset(rx_by_time,rx==2) # Control
rx3 <- subset(rx_by_time,rx==3) # FamilyNow to look at these…what do you expect based on this?
par(mfrow=c(1,1))
plot(rx2$time,rx2$weight,
type='b',
col="cyan",
main="Weight Change by Treatment",
ylim=c(80,95),
xlim=c(1,2),
xlab="When Weight Measured",
ylab="Weight",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
lines(rx1$time,rx1$weight,type="b",col="blue")
lines(rx3$time,rx3$weight,type="b",col="red")
legend(1.02,95,title="Treatment",legend=c("Control", "Cognitive", "Family"), col=c("cyan","blue","red"),lty=c(1,1,1),cex=0.8)
Maximum Likelihood Estimation
If we are taking a Bayesian approach, we would not look at MLE; however, for pedagogical purposes, we do it here. Before using Bayesian estimation, we’ll see what restricted maximum likelihood (REML) yields. I choose REML since these are unbiased; however, if we’re only looking at lmer models from package lme4, I would normally choose MLE. Using the lme4 package.
########################################################
# Simple Model using lme4 --- just a random intercept #
########################################################
# Remember: treatment = 1 Cognitive
# = 2 Control
# = 3 Family
summary(model1reml <- lmer(weight ~ time + (1 | girl), data=ano, REML=TRUE))Here are the fitted values from model 0
girl <- data.frame(seq(1:n))
par(mfrow=c(1,1))
b1 <- data.frame(coef(model1reml)[[1]])
tmp1 <- cbind(b1,girl)
names(tmp1) <- c("b0","b1","girl")
ano <- merge(ano,tmp1,by="girl")
ano$y.1 <- ano$b0 + ano$b1*ano$time
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.1,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="weight ~ time + (1 | girl)",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time, g$y.1, col=i)
}
Now make the model a bit more complex by adding in treatment
########################################################
#
# Model 2: weight ~ time + treatment + (1 | girl) #
########################################################
ano$treatment <- as.factor(ano$Rx)
model2reml <- lmer(weight ~ time + treatment + (1 | girl), data=ano, REML=TRUE)
summary(model2reml)Now graphing conditional values
b2 <- data.frame(coef(model2reml)[[1]])
tmp2 <- cbind(b2,girl)
names(tmp2) <- c("b0.2","b1.2","b2.2","b3.2","girl")
ano <- merge(ano,tmp2,by="girl")
ano$y.2 <- ano$b0.2 + ano$b1.2*ano$time + ano$b2.2*ano$Rx2 + ano$b3.2*ano$Rx3
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.2,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="weight ~ time + treatment + (1 | girl)",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time, g$y.2, col=i)
}
Doesn’t look much different, so let’s add an interaction between time and treatment. We see this is in the data.
#######################################################
# Model 3: weight ~ time + treatment + time*treatment + ( 1 | girl)
#
########################################################
model3reml <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl),data=ano,REML=TRUE)
summary(model3reml)The fitted values are
b3 <- data.frame(coef(model3reml)[[1]])
tmp3 <- cbind(b3,girl)
names(tmp3) <- c("b0.3","b1.3","b2.3","b3.3","b4.3","b5.3","girl")
ano <- merge(ano,tmp3,by="girl")
ano$y.3 <- ano$b0.3 + ano$b1.3*ano$time + ano$b2.3*ano$Rx2 + ano$b3.3*ano$Rx3 + ano$b4.3*ano$time*ano$Rx2 + ano$b5.3*ano$time*ano$Rx3
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.3,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="~time+treatment+time*treatment+( 1 | girl)",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time, g$y.3, col=i)
}
Based on our initial look at the data, we might want a random slope but NO random intercept.
########################################################
# Model 4: weight ~ time + treatment + time*treatment + ( 0 +time|girl)
########################################################
model4reml <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl) + ( 0 + time |girl),data=ano,REML=TRUE)## boundary (singular) fit: see ?isSingularmodel4reml <- lmer(weight ~ time + treatment + time*treatment + ( 0 + time |girl),data=ano,REML=TRUE)
summary(model4reml)What are fitted values now? These “look” like the data.
b4 <- data.frame(coef(model4reml)[[1]])
tmp4 <- cbind(b4,girl)
names(tmp4) <- c("b0.4","b1.4","b2.4","b3.4","b4.4","b5.4","girl")
ano <- merge(ano,tmp4,by="girl")
ano$y.4 <- ano$b0.4 + ano$b1.4*ano$time + ano$b2.4*ano$Rx2 + ano$b3.4*ano$Rx3 + ano$b4.4*ano$time*ano$Rx2 + ano$b5.4*ano$time*ano$Rx3
plot(ano$time,ano$y.4,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="~ ... +time*treatment+(0+time|girl)",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time, g$y.4, col=i)
}
## Summary of model fit via MLE
screenreg(c(model1reml, model2reml, model3reml, model4reml))##
## =========================================================================
## Model 1 Model 2 Model 3 Model 4
## -------------------------------------------------------------------------
## (Intercept) 79.64 *** 80.05 *** 79.68 *** 79.68 ***
## (1.55) (1.70) (2.30) (1.99)
## time 2.76 ** 2.76 ** 3.01 * 3.01 *
## (0.94) (0.94) (1.40) (1.35)
## treatment2 -2.86 * 2.32 2.32
## (1.38) (3.35) (2.89)
## treatment3 2.67 -3.72 -3.72
## (1.56) (3.78) (3.27)
## time:treatment2 -3.46 -3.46
## (2.03) (1.96)
## time:treatment3 4.26 4.26
## (2.30) (2.22)
## -------------------------------------------------------------------------
## AIC 955.48 942.76 929.86 922.05
## BIC 967.35 960.58 953.61 945.81
## Log Likelihood -473.74 -465.38 -456.93 -453.02
## Num. obs. 144 144 144 144
## Num. groups: girl 72 72 72 72
## Var: girl (Intercept) 13.84 10.04 11.80
## Var: Residual 31.87 31.87 28.34 22.91
## Var: girl time 6.89
## =========================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Bayesian estimation
To illustrate how to fit models, we start simple and work to more complex.
Model 0: The random intercept model
This model is also called the Null HLM and is the same as random effects ANOVA.
########################################################
# Bayesian Model 0 : random intercept no predictors (null hlm)
########################################################
dataList <- list(
y = ano$weight, # sd of weight (response)
sdY = sd(ano$weight), # sd of response
n = length(ano$weight), # number obs
ng= length(ano$weight)/2, # number of girls
girl= ano$girl
)The model part
model0 <- "model {
for (i in 1:n) { # likelihood
y[i] ~ dnorm(mu[i],precision)
mu[i] <- U0j[girl[i]] # note index
}
for (j in 1:ng) { # random intercepts
U0j[j] ~ dnorm(g0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2)) # other priors, mean of random
tau ~ dunif(0.0001,200) # sd of random
ptau <- 1/tau^2 # precision random
sigma ~ dunif(0.0001,2000) # sd of response (weight)
precision <- 1/sigma^2 # precision of weight/response
total.var <- sigma^2 + tau^2
icc <- tau^2/total.var
}"
writeLines(model0, con="model0.txt")I used different random number generators to show you what is available. This shouldn’t make a difference in results.
start1 = list("g0"=mean(ano$weight), "sigma"=sd(ano$weight), "tau"=.5, .RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=rnorm(1,0,3), "sigma"=5, "tau"=1, .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4), "sigma"=10, "tau"=5, .RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10), "sigma"=50, "tau"=20, .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)And run jags
model0.runjags <- run.jags(model=model0,
method="parallel",
monitor=c("g0", "sigma", "tau","total.var", "icc"),
data=dataList,
n.chains=4,
sample=50000,
burnin=5000,
inits=start,
thin=15) ## Calling 4 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 Fri Oct 14 10:55:34 2022
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 144
## Unobserved stochastic nodes: 75
## Total graph size: 382
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 5000
## -------------------------------------------------| 5000
## ************************************************** 100%
## . . . . . . Updating 750000
## -------------------------------------------------| 750000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 5 variables....
## Finished running the simulationplot(model0.runjags)## Generating plots...
print(model0.runjags)##
## JAGS model summary statistics from 200000 samples (thin = 15; chains = 4; adapt+burnin = 6000):
##
## Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD
## g0 82.529 83.787 85.081 83.785 0.64938 -- 0.0033175 0.5
## sigma 5.117 6.0904 7.1781 6.1253 0.53426 -- 0.0032203 0.6
## tau 0.74594 3.3152 5.0795 3.193 1.0639 -- 0.0087614 0.8
## total.var 37.551 48.616 61.401 49.132 6.1899 -- 0.031094 0.5
## icc 3.786e-07 0.2302 0.41532 0.22729 0.1157 -- 0.00082215 0.7
##
## SSeff AC.150 psrf
## g0 38316 0.0067444 1.0001
## sigma 27524 0.053373 0.99999
## tau 14745 0.16135 1.0001
## total.var 39628 0.0030098 0.99999
## icc 19806 0.095698 1.0001
##
## Total time taken: 1.7 minutessummary(model1reml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ time + (1 | girl)
## Data: ano
##
## REML criterion at convergence: 947.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92295 -0.46924 -0.03478 0.52092 2.24652
##
## Random effects:
## Groups Name Variance Std.Dev.
## girl (Intercept) 13.84 3.720
## Residual 31.87 5.645
## Number of obs: 144, groups: girl, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.6444 1.5509 99.2652 51.353 < 2e-16 ***
## time 2.7639 0.9409 71.0000 2.938 0.00446 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## time -0.910Compute ICC (intra-class correlation) and credible intervals
There are multiple ways that we can do this: (1) within jags (see above), (2) monte carlo after obtaining posteriors, (3) multiple posteriors which is illustrated below.
samps <-combine.mcmc(mcmc.objects = model0.runjags$mcmc)
# I need to know which columns I need
colnames(samps)## [1] "g0" "sigma" "tau" "total.var" "icc"sigma2 <- samps[,2]**2
tau2 <- samps[,3]**2
icc <- tau2/(sigma2+tau2)
summary(icc)##
## Iterations = 6001:3005986
## Thinning interval = 15
## Number of chains = 1
## Sample size per chain = 2e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.2264698 0.1154773 0.0002582 0.0006845
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.008047 0.143948 0.228794 0.308976 0.446422Model 1: random intercept
Now for our next model where we add in a predictor, “time”:
####################################################################
# Bayesian Model 1 : random intercept and fixed time #
####################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model1 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- b0j[girl[i]] + g1*time[i]
}
for (j in 1:ng) {
b0j[j] ~ dnorm(g0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
tau ~ dunif(0.0001,200)
ptau <- 1/tau^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model1, con="model1.txt")
start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3), "sigma"=sd(ano$weight), "tau"=.5,
.RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,0,3), "sigma"=5, "tau"=1,
.RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,3,4), "sigma"=50, "tau"=3,
.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,-3,10), "sigma"=10, "tau"=10,
.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
start <- list(start1,start2,start3,start4)
model1.runjags <- run.jags(model=model1,
method="parallel",
monitor=c("g0", "g1", "sigma", "tau"),
data=dataList,
sample=20000,
n.chains=4,
thin=10,
inits=start)## Calling 4 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 Fri Oct 14 10:57:25 2022
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 144
## Unobserved stochastic nodes: 76
## Total graph size: 671
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 4000
## -------------------------------------------------| 4000
## ************************************************** 100%
## . . . . . Updating 200000
## -------------------------------------------------| 200000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 4 variables....
## Finished running the simulation# And checking a bit more whether model1.runjags converged
plot(model1.runjags)## Generating plots...
print(model1.runjags)##
## JAGS model summary statistics from 80000 samples (thin = 10; chains = 4; adapt+burnin = 5000):
##
## Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD SSeff
## g0 76.53 79.613 82.764 79.606 1.589 -- 0.012911 0.8 15147
## g1 0.85913 2.7778 4.6803 2.7836 0.97037 -- 0.007939 0.8 14940
## sigma 4.8576 5.7862 6.8531 5.8217 0.51778 -- 0.0038384 0.7 18197
## tau 1.5737 3.5984 5.381 3.5214 0.94995 -- 0.010339 1.1 8441
##
## AC.100 psrf
## g0 0.027198 1.0001
## g1 0.029574 1.0002
## sigma 0.050738 1.0004
## tau 0.14458 1.0012
##
## Total time taken: 38.5 seconds# For comparison
summary(model1reml) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ time + (1 | girl)
## Data: ano
##
## REML criterion at convergence: 947.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92295 -0.46924 -0.03478 0.52092 2.24652
##
## Random effects:
## Groups Name Variance Std.Dev.
## girl (Intercept) 13.84 3.720
## Residual 31.87 5.645
## Number of obs: 144, groups: girl, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.6444 1.5509 99.2652 51.353 < 2e-16 ***
## time 2.7639 0.9409 71.0000 2.938 0.00446 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## time -0.910After fitting model, I decided I want an ICC, so…
# Compute residual ICC (intra-class correlation) and credible interval
samps <-combine.mcmc(mcmc.objects = model1.runjags$mcmc)
colnames(samps)## [1] "g0" "g1" "sigma" "tau"sigma2 <- samps[,3]**2
tau2 <- samps[,4]**2
icc1 <- tau2/(sigma2+tau2)
summary(icc1)##
## Iterations = 5001:804991
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 80000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.2764488 0.1151774 0.0004072 0.0011031
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.03616 0.19857 0.28118 0.35806 0.49009This model is too simple for the data, so we’ll go on to our next model, model2 which adds treatment as a fixed effect.
Model 2: Random intercept + fixed time and Rx
##################################################################
# Bayesian Model 2 : random intercept + fixed time + Rx #
##################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
rx1 = ano$Rx1,
rx3 = ano$Rx3,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model2 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- b0j[girl[i]] + g1*time[i] + g2*rx1[i] + g3*rx3[i]
}
for (j in 1:ng) {
b0j[j] ~ dnorm(g0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
g2 ~ dnorm(0,1/(100*sdY^2))
g3 ~ dnorm(0,1/(100*sdY^2))
tau ~ dunif(0.0001,200)
ptau <- 1/tau^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model2, con="model2.txt")
start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3), "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),
"sigma"=sd(ano$weight), "tau"=.5, .RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),
"sigma"=5, "tau"=1, .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,3,4), "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),
"sigma"=50, "tau"=3, .RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,-3,10), "g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
"sigma"=10, "tau"=10, .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
model2.runjags <- run.jags(model=model2,
method="parallel",
monitor=c("g0", "g1", "g2","g3","sigma", "tau"),
data=dataList,
sample=20000,
n.chains=4,
thin=10,
inits=start)## Calling 4 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 Fri Oct 14 10:58:10 2022
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 144
## Unobserved stochastic nodes: 78
## Total graph size: 965
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 4000
## -------------------------------------------------| 4000
## ************************************************** 100%
## . . . . . . . Updating 200000
## -------------------------------------------------| 200000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 6 variables....
## Finished running the simulation# And we should do some more model checking...
plot(model2.runjags)## Generating plots...
print(model2.runjags)##
## JAGS model summary statistics from 80000 samples (thin = 10; chains = 4; adapt+burnin = 5000):
##
## Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD SSeff
## g0 73.755 77.122 80.688 77.117 1.7592 -- 0.018509 1.1 9033
## g1 0.89853 2.7946 4.7061 2.794 0.97031 -- 0.0095332 1 10360
## g2 0.17535 2.8963 5.6152 2.8894 1.3849 -- 0.008085 0.6 29339
## g3 2.4247 5.5586 8.7142 5.56 1.6013 -- 0.0092449 0.6 30000
## sigma 4.913 5.8095 6.8547 5.8355 0.50442 -- 0.0044172 0.9 13041
## tau 0.53586 2.9675 4.6853 2.8494 1.0304 -- 0.014391 1.4 5127
##
## AC.100 psrf
## g0 0.11538 1.0001
## g1 0.093218 1.0001
## g2 0.017826 1
## g3 0.017716 1
## sigma 0.075106 1.0001
## tau 0.24454 1.0006
##
## Total time taken: 34.4 seconds# For comparison
summary(model2reml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ time + treatment + (1 | girl)
## Data: ano
##
## REML criterion at convergence: 930.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0131 -0.5380 -0.0562 0.6090 2.3744
##
## Random effects:
## Groups Name Variance Std.Dev.
## girl (Intercept) 10.04 3.168
## Residual 31.87 5.645
## Number of obs: 144, groups: girl, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 80.0473 1.6992 123.5077 47.108 < 2e-16 ***
## time 2.7639 0.9409 71.0000 2.938 0.00446 **
## treatment2 -2.8604 1.3764 69.0000 -2.078 0.04141 *
## treatment3 2.6687 1.5567 69.0000 1.714 0.09096 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time trtmn2
## time -0.831
## treatment2 -0.383 0.000
## treatment3 -0.339 0.000 0.418Compute residual ICC (intra-class correlation) and credible interval
samps <-combine.mcmc(mcmc.objects = model2.runjags$mcmc)
colnames(samps) # to see which columns I need## [1] "g0" "g1" "g2" "g3" "sigma" "tau"sigma2 <- samps[,5]**2
tau2 <- samps[,6]**2
icc2 <- tau2/(sigma2+tau2)
summary(icc2)##
## Iterations = 5001:804991
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 80000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.2080618 0.1153097 0.0004077 0.0013918
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.004971 0.121053 0.208733 0.290632 0.432762Is there a difference between Rx1 and Rx3 (cognitive and family)?
g2 <- samps[,3]
g3 <- samps[,4]
diff <- g3-g2
summary(diff)##
## Iterations = 5001:804991
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 80000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 2.671969 1.566567 0.005539 0.005528
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## -0.3917 1.6203 2.6740 3.7198 5.7564hist(diff, main="Posterior of Differences between family & cognitive", col="cyan")
### What is the probability that family is better than cognitive?
prob <- mean(diff>0)These results should not be treated as final, because we have more models to fit. However, this does illustrate how to be get postiors samples, compute variaous posteriors, get probabilities for hypothses, etc.
Model 3: Random intercept + fixed time, Rx & time*Rx
Next model (this will run fine but it is not run here)
##################################################################################
# Bayesian Model 3 : random intercept + fixed time + Rx + time*Rx #
##################################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
rx1 = ano$Rx1,
rx3 = ano$Rx3,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model3 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- b0j[girl[i]] + g1*time[i] + g2*rx1[i] + g3*rx3[i] + g4*time[i]*rx1[i] + g5*time[i]*rx3[i]
}
for (j in 1:ng) {
b0j[j] ~ dnorm(g0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
g2 ~ dnorm(0,1/(100*sdY^2))
g3 ~ dnorm(0,1/(100*sdY^2))
g4 ~ dnorm(0,1/(100*sdY^2))
g5 ~ dnorm(0,1/(100*sdY^2))
tau ~ dunif(0.0001,200)
ptau <- 1/tau^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model3, con="model3.txt")
start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3), "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),
"g4"=rnorm(1,0,3),"g5"=rnorm(1,0,3), "sigma"=sd(ano$weight), "tau"=.5,
.RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,-1,3), "g5"=rnorm(1,0,3), "sigma"=5, "tau"=1,
.RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,3,4), "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,3,4), "g5"=rnorm(1,0,3), "sigma"=50, "tau"=5,
.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10),"g1"=rnorm(1,-3,10), "g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,2,5),
"g5"=rnorm(1,3,3), "sigma"=10, "tau"=10,
.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
model3.runjags <- run.jags(model=model3,
method="parallel",
monitor=c("g0","g1","g2","g3","g4","g5","sigma","tau"),
data=dataList,
sample=20000,
n.chains=4,
thin=20,
inits=start)
print(model3.runjags)
summary(model3.reml)
plot(model3.runjags)Compute Residuals ICC (intra-class correlation) when using runjags
samps <-combine.mcmc(mcmc.objects = model3.runjags$mcmc)
colnames(samps) # to see which columns I need
sigma2 <- samps[,7]**2
tau2 <- samps[,8]**2
icc3 <- tau2/(sigma2+tau2)
summary(icc3)Model 4: Random slope, fixed time, Rx & time*Rx
Now we’ll switch to random slope and drop the random intercept. This might be OK. This code will run but not run here.
##################################################################################
# Bayesian Model 4 : random slope + fixed time + Rx + time*Rx #
##################################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
rx1 = ano$Rx1,
rx3 = ano$Rx3,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model4 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i] + g2*rx1[i] + g3*rx3[i]
+ g4*time[i]*rx1[i] + g5*time[i]*rx3[i]
}
for (j in 1:ng) {
b1j[j] ~ dnorm(0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
g2 ~ dnorm(0,1/(100*sdY^2))
g3 ~ dnorm(0,1/(100*sdY^2))
g4 ~ dnorm(0,1/(100*sdY^2))
g5 ~ dnorm(0,1/(100*sdY^2))
tau ~ dunif(0.0001,200)
ptau <- 1/tau^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model4, con="model4.txt")
start1 = list("g0"=mean(ano$weight), "g1"=rnorm(1,1,3), "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),
"g4"=rnorm(1,0,3), "g5"=rnorm(1,0,3), "sigma"=sd(ano$weight), "tau"=.5,
.RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,-2,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,-1,3),"g5"=rnorm(1,0,3), "sigma"=5, "tau"=1,
.RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,0,3), "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,3,4), "g5"=rnorm(1,0,3), "sigma"=50, "tau"=5,
.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,5,3),"g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
"g4"=rnorm(1,2,5), "g5"=rnorm(1,3,3), "sigma"=10, "tau"=10,
.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
model4.runjags <- run.jags(model=model4,
method="parallel",
monitor=c("g0","g1","g2","g3","g4","g5", "sigma", "tau"),
data=dataList,
sample=20000,
n.chains=4,
thin=5,
inits=start)
plot(model4.runjags)
print(model4.runjags)
summary(model4reml) Model 5: Random intecept & Slope with 2 time points!
We’ll try a random intercept and slope. With only 2 time points, you can’t fit both random intercept and slope as well as correlation between the random intercept and slope. In this example, \[\tau_01=0\]. This seems to be just enough to get the model to fit. Let’s see what happens when you try to get estimate both.
##################################################################################
# Bayesian Model 5: random intercept & slope + fixed time + Rx1 + Rx3 #
# density of tau_0 doesn't look good #
##################################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
rx1 = ano$Rx1,
rx3 = ano$Rx3,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model5 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- g0 + g1*time[i] + g2*rx1[i] + g3*rx3[i] + U0j[girl[i]] + U1j[girl[i]]*time[i]
}
for (j in 1:ng) {
U0j[j] ~ dnorm(0,ptau0)
U1j[j] ~ dnorm(0,ptau1)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
g2 ~ dnorm(0,1/(100*sdY^2))
g3 ~ dnorm(0,1/(100*sdY^2))
ptau0 ~ dgamma(0.001,0.001)
tau0 <- 1/sqrt(ptau0)
tau1 ~ dunif(0.0001,200)
ptau1 <- 1/tau1^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model5, con="model5.txt")
start1 = list("g0"=mean(ano$weight), "g1"=rnorm(1,1,3), "g2"=rnorm(1,1,3), "g3"=rnorm(1,4,3), "sigma"=sd(ano$weight), "ptau0"=.005, "tau1"=2,
.RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=dnorm(1,0,3), "g1"=rnorm(1,-2,3), "g2"=rnorm(1,0,3), "g3"=rnorm(1,-1,3), "sigma"=5, "ptau0"=.1, "tau1"=3,
.RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=dnorm(1,3,4), "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,4), "sigma"=50, "ptau0"=.045, "tau1"=10,
.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=dnorm(1,-3,10), "g1"=rnorm(1,5,3), "g2"=rnorm(1,5,3), "g3"=rnorm(1,3,4), "sigma"=10, "ptau0"=.10, "tau1"=1,
.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
model5.runjags <- run.jags(model=model5,
method="parallel",
monitor=c("g0","g1","g2","g3","sigma","ptau0","ptau1","tau0", "tau1"),
data=dataList,
sample=20000,
n.chains=4,
thin=20,
inits=start)
print(model5.runjags)
plot(model5.runjags)Refinement of Model 4
So, we’ll take Model 4 as our final model and get some additional statistics from the code (i.e., draws from posterior distribution).
###############################################################################
# Bayesian Model 4 refinement : random slope + fixed time + rx3+ time*Rx3 #
# 1st: drop time*Rx1 #
# 2nd: drop Rx1 #
###############################################################################
dataList <- list(
y = ano$weight,
time = ano$time,
rx3 = ano$Rx3,
sdY = sd(ano$weight),
n = length(ano$weight),
ng= length(ano$weight)/2,
girl= ano$girl
)
model4 <- "model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i],precision)
mu[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i] + g3*rx3[i] + g5*time[i]*rx3[i]
yhat[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i] + g3*rx3[i] + g5*time[i]*rx3[i]
}
for (j in 1:ng) {
b1j[j] ~ dnorm(0,ptau)
}
g0 ~ dnorm(0,1/(100*sdY^2))
g1 ~ dnorm(0,1/(100*sdY^2))
g3 ~ dnorm(0,1/(100*sdY^2))
g5 ~ dnorm(0,1/(100*sdY^2))
tau ~ dunif(0.0001,200)
ptau <- 1/tau^2
sigma ~ dunif(0.0001,2000)
precision <- 1/sigma^2
}"
writeLines(model4, con="model4.txt")
start1 = list("g0"=mean(ano$weight), "g1"=dnorm(1,1,3), "g3"=dnorm(1,0,3),
"g5"=dnorm(1,0,3), "sigma"=sd(ano$weight), "tau"=.5,
.RNG.name="base::Wichmann-Hill", .RNG.seed=523)
start2 = list("g0"=dnorm(1,0,3), "g1"=dnorm(1,-2,3), "g3"=dnorm(1,0,3),
"g5"=dnorm(1,0,3), "sigma"=5, "tau"=1,
.RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=dnorm(1,3,4), "g1"=dnorm(1,0,3), "g3"=dnorm(1,0,3),
"g5"=dnorm(1,0,3), "sigma"=50, "tau"=5,
.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=dnorm(1,-3,10), "g1"=dnorm(1,5,3),"g3"=dnorm(1,0,3),
"g5"=dnorm(1,3,3), "sigma"=10, "tau"=10,
.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
model4.runjags <- run.jags(model=model4,
method="parallel",
monitor=c("g0","g1","g3","g5", "sigma", "tau","yhat"),
data=dataList,
sample=20000,
n.chains=4,
thin=5,
inits=start)## Calling 4 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 Fri Oct 14 10:58:54 2022
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 144
## Unobserved stochastic nodes: 78
## Total graph size: 967
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 4000
## -------------------------------------------------| 4000
## ************************************************** 100%
## . . . . . . . . Updating 100000
## -------------------------------------------------| 100000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Note: Summary statistics were not produced as there are >50 monitored
## variables
## [To override this behaviour see ?add.summary and ?runjags.options]
## FALSEFinished running the simulationprint(model4.runjags)##
## JAGS model with 80000 samples (thin = 5; chains = 4; adapt+burnin = 5000)
##
## Full summary statistics have not been pre-calculated - use either the summary method or add.summary to calculate summary statisticssummary(model4reml) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ time + treatment + time * treatment + (0 + time | girl)
## Data: ano
##
## REML criterion at convergence: 906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59369 -0.46354 -0.00561 0.49016 1.64149
##
## Random effects:
## Groups Name Variance Std.Dev.
## girl time 6.891 2.625
## Residual 22.913 4.787
## Number of obs: 144, groups: girl, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.683 1.988 69.000 40.090 <2e-16 ***
## time 3.007 1.348 104.632 2.230 0.0279 *
## treatment2 2.325 2.891 69.000 0.804 0.4240
## treatment3 -3.718 3.270 69.000 -1.137 0.2594
## time:treatment2 -3.457 1.961 104.632 -1.763 0.0808 .
## time:treatment3 4.258 2.218 104.632 1.920 0.0576 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time trtmn2 trtmn3 tm:tr2
## time -0.885
## treatment2 -0.688 0.608
## treatment3 -0.608 0.538 0.418
## tim:trtmnt2 0.608 -0.688 -0.885 -0.370
## tim:trtmnt3 0.538 -0.608 -0.370 -0.885 0.418Why didn’t print(model4.runjags) give you the table of statistics? There are many of them; in particular, there are an extra 72 X 2 = 143 yhat values, 2 values per girl (weight before and weight after).
plot(model4.runjags)Let’s do some other things with posterior samples. This is a better (more straightforward) than previous version that I had here.
# Takes output and put in format so can use coda plots and diagnostics on it.
samps <- combine.mcmc(model4.runjags)
colnames(samps)## [1] "g0" "g1" "g3" "g5" "sigma" "tau"
## [7] "yhat[1]" "yhat[2]" "yhat[3]" "yhat[4]" "yhat[5]" "yhat[6]"
## [13] "yhat[7]" "yhat[8]" "yhat[9]" "yhat[10]" "yhat[11]" "yhat[12]"
## [19] "yhat[13]" "yhat[14]" "yhat[15]" "yhat[16]" "yhat[17]" "yhat[18]"
## [25] "yhat[19]" "yhat[20]" "yhat[21]" "yhat[22]" "yhat[23]" "yhat[24]"
## [31] "yhat[25]" "yhat[26]" "yhat[27]" "yhat[28]" "yhat[29]" "yhat[30]"
## [37] "yhat[31]" "yhat[32]" "yhat[33]" "yhat[34]" "yhat[35]" "yhat[36]"
## [43] "yhat[37]" "yhat[38]" "yhat[39]" "yhat[40]" "yhat[41]" "yhat[42]"
## [49] "yhat[43]" "yhat[44]" "yhat[45]" "yhat[46]" "yhat[47]" "yhat[48]"
## [55] "yhat[49]" "yhat[50]" "yhat[51]" "yhat[52]" "yhat[53]" "yhat[54]"
## [61] "yhat[55]" "yhat[56]" "yhat[57]" "yhat[58]" "yhat[59]" "yhat[60]"
## [67] "yhat[61]" "yhat[62]" "yhat[63]" "yhat[64]" "yhat[65]" "yhat[66]"
## [73] "yhat[67]" "yhat[68]" "yhat[69]" "yhat[70]" "yhat[71]" "yhat[72]"
## [79] "yhat[73]" "yhat[74]" "yhat[75]" "yhat[76]" "yhat[77]" "yhat[78]"
## [85] "yhat[79]" "yhat[80]" "yhat[81]" "yhat[82]" "yhat[83]" "yhat[84]"
## [91] "yhat[85]" "yhat[86]" "yhat[87]" "yhat[88]" "yhat[89]" "yhat[90]"
## [97] "yhat[91]" "yhat[92]" "yhat[93]" "yhat[94]" "yhat[95]" "yhat[96]"
## [103] "yhat[97]" "yhat[98]" "yhat[99]" "yhat[100]" "yhat[101]" "yhat[102]"
## [109] "yhat[103]" "yhat[104]" "yhat[105]" "yhat[106]" "yhat[107]" "yhat[108]"
## [115] "yhat[109]" "yhat[110]" "yhat[111]" "yhat[112]" "yhat[113]" "yhat[114]"
## [121] "yhat[115]" "yhat[116]" "yhat[117]" "yhat[118]" "yhat[119]" "yhat[120]"
## [127] "yhat[121]" "yhat[122]" "yhat[123]" "yhat[124]" "yhat[125]" "yhat[126]"
## [133] "yhat[127]" "yhat[128]" "yhat[129]" "yhat[130]" "yhat[131]" "yhat[132]"
## [139] "yhat[133]" "yhat[134]" "yhat[135]" "yhat[136]" "yhat[137]" "yhat[138]"
## [145] "yhat[139]" "yhat[140]" "yhat[141]" "yhat[142]" "yhat[143]" "yhat[144]"g0 <- samps[,1]
g1 <- samps[,2]
g3 <- samps[,3]
g5 <- samps[,4]
sigma <- samps[,5]
tau <- samps[,6]
fitted <- samps[,7:150]
ano$yhat <- apply(fitted,2,"mean")
plot(ano$time,ano$yhat,
type="n",
ylim=c(65,104),
xlab="Before/After Treatment",
ylab="Weight in Pounds",
main="Bayesian refined Model 4",
xaxt="n"
)
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))
for (i in 1:n){
g <- subset(ano,ano$girl==i)
lines(g$time, g$yhat, col=i)
}
This graph of fitted values resembles pattern in the data.
More methods for checked convergence should be done.