Linear mixed model MCMC
2024-09-05
Linear mixed model in JAGS for random intercept
library(R2jags)## Loading required package: rjags## Loading required package: coda## Linked to JAGS 4.3.1## Loaded modules: basemod,bugs##
## Attaching package: 'R2jags'## The following object is masked from 'package:coda':
##
## traceplot# create dataset
library(nlme)
snakes3=Orthodont
snakes3$bw = rep(c(1:27),each=4)
names(snakes3)=c( "b_mass" ,"b_length","site" ,"sex","bw" )
snakes3=snakes3[sample(nrow(snakes3)*0.8 ), ]
head(snakes3)## b_mass b_length site sex bw
## 37 27.5 8 M10 Male 10
## 48 28.0 14 M12 Male 12
## 27 24.5 12 M07 Male 7
## 12 27.5 14 M03 Male 3
## 74 24.0 10 F03 Female 19
## 59 26.0 12 M15 Male 15 table(snakes3$site)##
## M16 M05 M02 M11 M07 M08 M03 M12 M13 M14 M09 M15 M06 M04 M01 M10 F10 F09 F06 F01
## 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 0 0 2 4
## F05 F07 F02 F08 F03 F04 F11
## 4 0 4 0 4 4 0 nrow(snakes3)## [1] 86 plot(b_mass ~ b_length, col = site, data = snakes3)
# snakes3 <- within(snakes3, sex <- relevel(sex, ref = "Female"))
summary(lme(b_mass ~ b_length + (sex), data = snakes3, random = ~ 1|site))## Linear mixed-effects model fit by REML
## Data: snakes3
## AIC BIC logLik
## 360.7635 372.8577 -175.3817
##
## Random effects:
## Formula: ~1 | site
## (Intercept) Residual
## StdDev: 1.52578 1.52631
##
## Fixed effects: b_mass ~ b_length + (sex)
## Value Std.Error DF t-value p-value
## (Intercept) 16.918734 0.9199308 63 18.391311 0.0000
## b_length 0.731820 0.0740991 63 9.876227 0.0000
## sexFemale -2.009645 0.8253985 20 -2.434758 0.0244
## Correlation:
## (Intr) b_lngt
## b_length -0.886
## sexFemale -0.262 0.026
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.40584199 -0.58577184 -0.03230801 0.54900703 3.39016809
##
## Number of Observations: 86
## Number of Groups: 22 # summary(lme(b_mass ~ b_length , data = snakes3, random = ~ 1|site)) Nsites <- length(levels(as.factor(snakes3$site)))
jagsdata_s3 <- with(snakes3, list(b_mass = b_mass, b_length = b_length, site = site, sex=sex,
N = length(b_mass), Nsites = Nsites))
# The JAGS model
lm3_jags <- function(){
# Likelihood:
for (i in 1:N){
b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
mu[i] <- alpha + a[site[i]] + beta * b_length[i]+ gama*sex[i] # Random intercept for site
}
# Priors:
alpha ~ dnorm(0, 0.01) # intercept
sigma_a ~ dunif(0, 100) # standard deviation of random effect (variance between sites)
tau_a <- 1 / (sigma_a * sigma_a) # convert to precision
for (j in 1:Nsites){
a[j] ~ dnorm(0, tau_a) # random intercept for each site
}
beta ~ dnorm(0, 0.01) # slope
gama ~ dnorm(0, 0.01) # slope
sigma ~ dunif(0, 100) # standard deviation of fixed effect (variance within sites)
tau <- 1 / (sigma * sigma) # convert to precision
}
init_values <- function(){
list(alpha = rnorm(1), sigma_a = runif(1), beta = rnorm(1), sigma = runif(1),gama= rnorm(1) )
}
params <- c("alpha", "beta", "sigma", "sigma_a", "gama")
fit_lm3 <- jags(data = jagsdata_s3, inits = init_values, parameters.to.save = params, model.file = lm3_jags,
n.chains = 3, n.iter = 20000, n.burnin = 5000, n.thin = 10, DIC = F)## module glm loaded## module dic loaded## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 86
## Unobserved stochastic nodes: 32
## Total graph size: 478
##
## Initializing model fit_lm3## Inference for Bugs model at "C:/Users/hed2/AppData/Local/Temp/RtmpYBOBop/model39446ec37b3f", fit using jags,
## 3 chains, each with 20000 iterations (first 5000 discarded), n.thin = 10
## n.sims = 4500 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## alpha 18.455 1.460 15.570 17.499 18.484 19.429 21.276 1.002 2000
## beta 0.746 0.076 0.595 0.695 0.746 0.796 0.898 1.001 3100
## gama -1.781 0.886 -3.489 -2.367 -1.787 -1.216 -0.021 1.001 3000
## sigma 1.561 0.147 1.311 1.459 1.550 1.651 1.883 1.001 4500
## sigma_a 1.626 0.350 1.052 1.383 1.590 1.825 2.411 1.001 2800
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).Linear mixed model in JAGS for random slope
we suppose they are independent of random intercept and slope
(lme(b_mass ~ b_length + (sex), data = snakes3, random = ~ 1+sex|site))## Linear mixed-effects model fit by REML
## Data: snakes3
## Log-restricted-likelihood: -174.9241
## Fixed: b_mass ~ b_length + (sex)
## (Intercept) b_length sexFemale
## 16.896158 0.733872 -1.990918
##
## Random effects:
## Formula: ~1 + sex | site
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.662072 (Intr)
## sexFemale 1.705354 -0.827
## Residual 1.526174
##
## Number of Observations: 86
## Number of Groups: 22 # summary(lme(b_mass ~ b_length , data = snakes3, random = ~ 1|site))
Nsites <- length(levels(as.factor(snakes3$site)))
jagsdata_s3 <- with(snakes3, list(b_mass = b_mass, b_length = b_length, site = site, sex=sex,
N = length(b_mass), Nsites = Nsites))
# The JAGS model
lm3_jags <- function(){
# Likelihood:
for (i in 1:N){
b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
mu[i] <- (alpha + a[site[i]] )+ beta * b_length[i]+ (gama+g[site[i]]) *sex[i] # Random intercept for site
}
# Priors:
alpha ~ dnorm(0, 0.01) # intercept
sigma_a ~ dunif(0, 100) # standard deviation of random effect (variance between sites)
tau_a <- 1 / (sigma_a * sigma_a) # convert to precision
for (j in 1:Nsites){
a[j] ~ dnorm(0, tau_a) # random intercept for each site
}
sigma_g ~ dunif(0, 100) # standard deviation of random effect (variance between sites)
tau_g <- 1 / (sigma_g * sigma_g) # convert to precision
for (j in 1:Nsites){
g[j] ~ dnorm(0, tau_g) # random intercept for each site
}
beta ~ dnorm(0, 0.01) # slope
gama ~ dnorm(0, 0.01) # slope
sigma ~ dunif(0, 100) # standard deviation of fixed effect (variance within sites)
tau <- 1 / (sigma * sigma) # convert to precision
}
init_values <- function(){
list(alpha = rnorm(1), sigma_a = runif(1), sigma_g = runif(1),beta = rnorm(1), sigma = runif(1),gama= rnorm(1) )
}
params <- c("alpha", "beta", "sigma", "sigma_a","sigma_g", "gama")
fit_lm34 <- jags(data = jagsdata_s3, inits = init_values, parameters.to.save = params, model.file = lm3_jags,
n.chains = 3, n.iter = 20000, n.burnin = 5000, n.thin = 10, DIC = F)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 86
## Unobserved stochastic nodes: 60
## Total graph size: 572
##
## Initializing model fit_lm34## Inference for Bugs model at "C:/Users/hed2/AppData/Local/Temp/RtmpYBOBop/model39443b7d1b3c", fit using jags,
## 3 chains, each with 20000 iterations (first 5000 discarded), n.thin = 10
## n.sims = 4500 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## alpha 18.400 1.581 15.234 17.365 18.425 19.467 21.398 1.001 4500
## beta 0.745 0.075 0.597 0.696 0.744 0.794 0.891 1.001 4500
## gama -1.701 1.069 -3.685 -2.403 -1.746 -1.044 0.510 1.001 4500
## sigma 1.563 0.146 1.309 1.461 1.554 1.654 1.881 1.001 4500
## sigma_a 1.330 0.572 0.056 1.031 1.386 1.690 2.417 1.038 200
## sigma_g 0.718 0.509 0.044 0.297 0.618 1.062 1.824 1.011 290
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).The results are a little different
we suppose they are not independent of random intercept and slope
Nsites <- length(levels(as.factor(snakes3$site)))
jagsdata_s3 <- with(snakes3, list(b_mass = b_mass, b_length = b_length, site = site, sex=sex,
N = length(b_mass), Nsites = Nsites))
# The JAGS model
lm3_jags <- function(){
# Likelihood:
for (i in 1:N){
b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
mu[i] <- (alpha + a[site[i]] )+ beta * b_length[i]+ (gama+g[site[i]]) *sex[i] # Random intercept for site
}
# Priors:
alpha ~ dnorm(0, 0.01) # intercept
for (j in 1:Nsites){
## random effects
## mean is 0
b0[j,1] = 0
b0[j,2] = 0
b[j,1:2] ~ dmnorm(b0[j,1:2],ISigma[,]) ## MV normal
a[j] =b[j,1]
g[j]=b[j,2]
}
#(3) Variance-covariance matrix
ISigma[1:2,1:2] ~dwish(R[,],3)
Sigma[1:2,1:2] = inverse(ISigma[,])
R[1,1] =1
R[1,2] =0
R[2,2] =1
R[2,1] =0
beta ~ dnorm(0, 0.01) # slope
gama ~ dnorm(0, 0.01) # slope
sigma ~ dunif(0, 100) # standard deviation of fixed effect (variance within sites)
tau <- 1 / (sigma * sigma) # convert to precision
}
init_values <- function(){
list(alpha = rnorm(1), sigma_a = runif(1), sigma_g = runif(1),beta = rnorm(1), sigma = runif(1),gama= rnorm(1) )
}
params <- c("alpha", "beta", "sigma", "Sigma", "gama")
fit_lm35 <- jags(data = jagsdata_s3, inits = init_values, parameters.to.save = params, model.file = lm3_jags,
n.chains = 3, n.iter = 20000, n.burnin = 5000, n.thin = 10, DIC = F)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 86
## Unobserved stochastic nodes: 32
## Total graph size: 624## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 1## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 1## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 2## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 2## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 3## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 3## Initializing model fit_lm35## Inference for Bugs model at "C:/Users/hed2/AppData/Local/Temp/RtmpYBOBop/model3944234e60d0", fit using jags,
## 3 chains, each with 20000 iterations (first 5000 discarded), n.thin = 10
## n.sims = 4500 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## Sigma[1,1] 2.054 2.802 0.193 0.625 1.199 2.343 9.511 1.002 2000
## Sigma[2,1] -0.361 1.695 -4.439 -0.372 0.080 0.286 0.722 1.002 3000
## Sigma[1,2] -0.361 1.695 -4.439 -0.372 0.080 0.286 0.722 1.002 3000
## Sigma[2,2] 0.854 1.133 0.135 0.326 0.557 0.994 3.358 1.001 4300
## alpha 18.547 1.483 15.534 17.552 18.578 19.519 21.408 1.001 4500
## beta 0.743 0.076 0.594 0.692 0.744 0.792 0.893 1.001 4500
## gama -1.806 0.973 -3.687 -2.441 -1.823 -1.193 0.195 1.001 4500
## sigma 1.573 0.147 1.322 1.466 1.564 1.664 1.891 1.001 4500
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).similar results for the above scenarios
Join model for missing data
Nsites <- length(levels(as.factor(snakes3$site)))
jagsdata_s3 <- with(snakes3, list(b_mass = b_mass, b_length = b_length, site = site, sex=sex,
N = length(b_mass), Nsites = Nsites,bw=bw))
# The JAGS model
lm3_jags <- function(){
# Likelihood:
for (i in 1:N){
b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
mu[i] <- (alpha + a[site[i]] )+ beta * b_length[i]+ (gama+g[site[i]]) *sex[i] # Random intercept for site
## linear regression model
bw[i] ~ dnorm(m[i],1)
m[i] = bb[1] + bb[2]*a[site[i]] + bb[3]*g[site[i]]
}
# Priors:
alpha ~ dnorm(0, 0.01) # intercept
#(1) Coefficients
for (ll in 1:3) { bb[ll] ~ dnorm(0, .01) }
for (j in 1:Nsites){
## random effects
## mean is 0
b0[j,1] = 0
b0[j,2] = 0
b[j,1:2] ~ dmnorm(b0[j,1:2],ISigma[,]) ## MV normal
a[j] =b[j,1]
g[j]=b[j,2]
}
#(3) Variance-covariance matrix
ISigma[1:2,1:2] ~dwish(R[,],3)
Sigma[1:2,1:2] = inverse(ISigma[,])
R[1,1] =1
R[1,2] =0
R[2,2] =1
R[2,1] =0
beta ~ dnorm(0, 0.01) # slope
gama ~ dnorm(0, 0.01) # slope
sigma ~ dunif(0, 100) # standard deviation of fixed effect (variance within sites)
tau <- 1 / (sigma * sigma) # convert to precision
}
init_values <- function(){
list(alpha = rnorm(1), sigma_a = runif(1), sigma_g = runif(1),beta = rnorm(1), sigma = runif(1),gama= rnorm(1) ,bb= rnorm(3))
}
params <- c("alpha", "beta", "sigma", "Sigma", "gama","bb","a","g")
fit_lm35 <- jags(data = jagsdata_s3, inits = init_values, parameters.to.save = params, model.file = lm3_jags,
n.chains = 3, n.iter = 20000, n.burnin = 5000, n.thin = 10, DIC = F)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 172
## Unobserved stochastic nodes: 35
## Total graph size: 779## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 1## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 1## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 2## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 2## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_a" in chain 3## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : Unused initial value for "sigma_g" in chain 3## Initializing model fit_lm35## Inference for Bugs model at "C:/Users/hed2/AppData/Local/Temp/RtmpYBOBop/model3944117f7deb", fit using jags,
## 3 chains, each with 20000 iterations (first 5000 discarded), n.thin = 10
## n.sims = 4500 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## Sigma[1,1] 1.902 3.051 0.248 0.600 0.994 1.771 11.361 1.178 20
## Sigma[2,1] -0.250 1.699 -5.721 -0.094 0.172 0.366 0.896 1.353 22
## Sigma[1,2] -0.250 1.699 -5.721 -0.094 0.172 0.366 0.896 1.353 22
## Sigma[2,2] 0.895 1.007 0.170 0.363 0.571 0.995 3.726 1.063 42
## a[1] -1.023 0.958 -3.833 -1.311 -0.852 -0.494 0.309 1.286 15
## a[2] -0.652 0.997 -2.990 -1.206 -0.511 0.026 0.956 1.497 8
## a[3] -0.359 1.072 -2.705 -1.051 -0.174 0.388 1.370 1.888 5
## a[4] -0.549 0.720 -2.397 -0.844 -0.423 -0.117 0.485 1.125 41
## a[5] -0.354 0.792 -2.143 -0.748 -0.271 0.157 0.941 1.305 11
## a[6] -0.330 0.746 -2.070 -0.700 -0.236 0.139 0.846 1.196 16
## a[7] 0.014 0.934 -1.959 -0.607 0.147 0.623 1.644 1.807 6
## a[8] -0.288 0.659 -1.788 -0.575 -0.222 0.068 0.832 1.082 130
## a[9] -0.323 0.675 -2.049 -0.607 -0.259 0.030 0.891 1.119 42
## a[10] -0.039 0.654 -1.283 -0.394 -0.072 0.299 1.319 1.187 18
## a[11] 0.243 0.647 -0.964 -0.094 0.210 0.542 1.693 1.101 51
## a[12] 0.404 0.794 -0.911 -0.104 0.310 0.805 2.242 1.321 10
## a[13] 0.964 0.894 -0.392 0.458 0.836 1.267 3.484 1.283 14
## a[14] 1.157 1.038 -0.469 0.584 1.017 1.512 4.111 1.388 10
## a[15] 1.819 1.438 -0.360 1.051 1.593 2.208 6.212 1.462 9
## a[16] 2.349 1.676 -0.192 1.366 2.026 2.922 7.214 1.251 16
## a[17] 0.038 1.348 -2.547 -0.657 0.025 0.680 2.871 1.051 4500
## a[18] -0.021 1.408 -2.727 -0.694 -0.018 0.644 2.753 1.068 4500
## a[19] -0.917 1.373 -4.109 -1.470 -0.924 -0.265 1.510 1.468 8
## a[20] -1.048 1.112 -3.796 -1.423 -0.912 -0.494 0.580 1.133 69
## a[21] -0.697 1.272 -3.516 -1.224 -0.739 -0.081 1.645 1.459 9
## a[22] -0.013 1.372 -2.753 -0.667 -0.013 0.662 2.638 1.059 2900
## a[23] -0.426 1.027 -2.503 -0.864 -0.449 0.057 1.443 1.220 19
## a[24] -0.010 1.396 -2.746 -0.667 -0.012 0.676 2.748 1.062 2000
## a[25] -0.157 1.164 -2.044 -0.719 -0.290 0.413 2.084 1.347 11
## a[26] 0.269 1.338 -1.744 -0.520 0.073 0.966 3.211 1.467 8
## a[27] -0.034 1.339 -2.689 -0.687 0.012 0.644 2.473 1.065 3200
## alpha 17.899 1.954 13.969 16.629 17.926 19.218 21.694 1.106 24
## bb[1] 11.267 1.487 8.343 10.303 11.271 12.235 14.115 1.004 590
## bb[2] -3.166 6.909 -14.473 -8.617 -4.712 2.863 9.321 2.204 5
## bb[3] 1.729 10.079 -17.134 -8.334 5.366 9.511 16.188 3.853 3
## beta 0.743 0.076 0.596 0.691 0.744 0.794 0.892 1.003 1000
## g[1] -0.308 0.799 -1.629 -0.727 -0.402 -0.063 2.039 1.361 11
## g[2] -0.544 0.826 -2.097 -1.075 -0.623 0.080 1.027 1.834 6
## g[3] -0.531 0.942 -2.213 -1.159 -0.696 0.256 1.236 2.470 4
## g[4] -0.263 0.539 -1.367 -0.542 -0.236 0.013 0.897 1.126 29
## g[5] -0.331 0.597 -1.549 -0.713 -0.349 0.109 0.731 1.594 7
## g[6] -0.268 0.548 -1.420 -0.596 -0.269 0.107 0.702 1.373 9
## g[7] -0.294 0.791 -1.771 -0.844 -0.402 0.363 1.160 2.344 4
## g[8] -0.112 0.460 -1.062 -0.344 -0.112 0.106 0.838 1.064 82
## g[9] -0.087 0.475 -1.006 -0.338 -0.106 0.146 1.024 1.115 29
## g[10] 0.086 0.473 -0.748 -0.221 0.051 0.349 1.121 1.264 12
## g[11] 0.040 0.473 -0.954 -0.220 0.065 0.311 0.931 1.139 21
## g[12] 0.344 0.610 -0.738 -0.081 0.336 0.706 1.629 1.518 7
## g[13] 0.275 0.786 -1.874 -0.024 0.364 0.692 1.707 1.395 10
## g[14] 0.301 0.935 -2.288 -0.034 0.432 0.825 1.840 1.517 8
## g[15] 0.492 1.342 -3.505 0.097 0.722 1.225 2.567 1.540 8
## g[16] 1.092 1.520 -3.253 0.622 1.242 1.860 3.711 1.259 16
## g[17] 0.007 0.927 -1.878 -0.518 0.001 0.527 1.878 1.010 4500
## g[18] -0.009 0.951 -1.963 -0.517 0.000 0.506 1.855 1.016 2600
## g[19] -0.106 1.061 -1.893 -0.852 -0.145 0.480 2.415 1.840 6
## g[20] -0.336 0.837 -1.821 -0.794 -0.404 0.027 1.760 1.354 11
## g[21] -0.030 0.963 -1.680 -0.715 -0.017 0.523 2.169 1.889 5
## g[22] -0.004 0.945 -1.928 -0.524 0.009 0.517 1.892 1.016 4500
## g[23] -0.012 0.752 -1.349 -0.529 -0.038 0.445 1.547 1.657 6
## g[24] -0.027 0.969 -1.919 -0.536 -0.019 0.506 1.904 1.017 3600
## g[25] 0.148 0.815 -1.282 -0.461 0.196 0.677 1.654 1.866 5
## g[26] 0.373 0.900 -1.329 -0.316 0.498 0.979 1.994 2.058 5
## g[27] 0.020 0.941 -1.852 -0.486 0.015 0.527 1.918 1.010 1500
## gama -1.302 1.523 -4.154 -2.380 -1.313 -0.277 1.806 1.217 13
## sigma 1.586 0.152 1.319 1.483 1.577 1.680 1.921 1.002 1700
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).table(snakes3$site)##
## M16 M05 M02 M11 M07 M08 M03 M12 M13 M14 M09 M15 M06 M04 M01 M10 F10 F09 F06 F01
## 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 0 0 2 4
## F05 F07 F02 F08 F03 F04 F11
## 4 0 4 0 4 4 0