11: Childhood Undernutrition in India

Since the data set used in Chapter 11 is proprietary (see https://dhsprogram.com/ for details), we are not able to provide it. In Section 1 we give the coding and output of the analysis in Chapter 11; in Section 2 we provide reproducible code, using a data set that is in the public domain, that mimics the analysis in Chapter 11.

Data analysis in Chapter 11

library(parallel)
library(BayesX)
library(bamlss)
library(modelr)
library(latex2exp)

wd <- paste0(getwd(), "/") 
setwd(wd)
pathdata <- paste(wd, "data/", sep = "")

# load preprocessed data

data <- read.table("data/india.raw", header = TRUE)
data$state <- as.factor(data$state)

K <- read.gra("data/K_penalty_connected.gra")

Note: map consists of 36 regions
Reading map ... finished

11.2 Univariate Analysis

11.2.1 Model Estimation

rerun <- FALSE

if (!rerun)  {

  load("data/univnormal_models.RData")

} else {

  set.seed(42)

  # model with constant variance ----------

  f1 <- wasting2 ~ s(state, bs = "mrf", xt = list(penalty = K, prior = "ig",
                                                  a = 0.001, b = 0.001)) +
                   csex +
                   s(cage, xt = list(prior = "ig", a = 0.001, b = 0.001)) +
                   s(cage, by = csex, xt = list(prior = "ig", a = 0.001, b = 0.001))

  # Estimate model

  m1 <- bamlss(f1, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

  # model with varying variance ----------

  f2 <- list(wasting2 ~ s(state, bs = "mrf", xt = list(penalty = K, prior = "ig",
                                                       a = 0.001, b = 0.001)) +
                        csex +
                        s(cage, xt = list(prior = "ig", a = 0.001, b = 0.001)) +
                        s(cage, by = csex, xt = list(prior = "ig",
                                                     a = 0.001, b = 0.001)),
             sigma ~ s(state, bs = "mrf", xt = list(penalty = K)) +
                     csex +
                     s(cage) +
                     s(cage, by = csex))

  # Estimate model

  m2 <- bamlss(f2, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)


### Estimate models with different hyperparameters for tau^2


  f_sd <- list(wasting2 ~ s(state, bs = "mrf", xt = list(penalty = K, prior = "sd")) +
                           csex +
                           s(cage, xt = list(prior = "sd")) +
                           s(cage, by = csex, xt = list(prior = "sd")),
                sigma ~ s(state, bs = "mrf", xt = list(penalty = K)) +
                        csex +
                        s(cage) +
                        s(cage, by = csex))

  f_hc <- list(wasting2 ~ s(state, bs = "mrf", xt = list(penalty = K, prior = "hc")) +
                          csex +
                          s(cage, xt = list(prior = "hc")) +
                          s(cage, by = csex, xt = list(prior = "hc")),
               sigma ~ s(state, bs = "mrf", xt = list(penalty = K)) +
                       csex +
                       s(cage) +
                       s(cage, by = csex))

  f_u <- list(wasting2 ~ s(state, bs = "mrf", xt = list(penalty = K, prior = "u")) +
                         csex +
                         s(cage, xt = list(prior = "u")) +
                         s(cage, by = csex, xt = list(prior = "u")),
              sigma ~ s(state, bs = "mrf", xt = list(penalty = K)) +
                      csex +
                      s(cage) +
                      s(cage, by = csex))

  # Estimate models

  m_ig <- m2

  m_sd <- bamlss(f_sd, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

  m_hc <- bamlss(f_hc, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

  m_u <- bamlss(f_u, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

  save(list = setdiff(ls(), c("wd", "pathdata", "pathplots")),
       file = paste(pathdata,"univnormal_models.RData", sep = ""))
}

11.2.2 MCMC Diagnostics

Effective sample size

samples <- data.frame(m1[["samples"]][[1]])$mu.s.s.state..x33
samples <- as.mcmc(samples)
effectiveSize(samples)

var1 
1001 

Individual Plots

# trajectory ######

par(mar=c(4.1, 6.3, 4.1, 1.1))
plot(1:1001,data.frame(m1[["samples"]][[1]])$mu.s.s.state..x33, type = "l",
     main = TeX(r"((a) MCMC sampling path of $\beta^{\mu}_{{spat}_{33}}$)"),
     ylab = TeX(r"($\beta^{\mu}_{{spat}_{33}}$)"),
     xlab = "Iterations", cex.main = 1, cex.lab = 1, cex.axis = 0.75)

# acf ######

par(mar=c(4.5, 5, 4.1, 1.1))
autocorr.plot(samples, auto.layout = T, lwd = 3,
              main = TeX(r"((b) ACF plot of $\beta^{\mu}_{{spat}_{33}}$)"),
              cex.main = 1, cex.lab = 1, cex.axis = 0.75)

11.2.3 Model Comparison

Normalized Quantile Residuals

pred <- predict(m1,newdata=data,what="samples",type = "parameter", FUN=identity)

param_sigma <- rowMeans(data.frame(pred["sigma"]))
param_mu <- rowMeans(data.frame(pred["mu"]) )

y <- data$wasting2

pn <- function (y,  param)
{
  mean <-  param[, 1]
  sd <- (param[, 2])
  pnorm(y, mean = mean, sd = sd)
}

resid_m1=qnorm(pn(y, cbind(param_mu, param_sigma)))


##### m2 --------------

pred <- predict(m2,newdata=data,what="samples",type = "parameter", FUN=identity)

param_sigma <- rowMeans(data.frame(pred["sigma"]))
param_mu <- rowMeans(data.frame(pred["mu"]) )

y <- data$wasting2

pn <- function (y,  param)
{
  mean <-  param[, 1]
  sd <- (param[, 2])
  pnorm(y, mean = mean, sd = sd)
}

resid_m2=qnorm(pn(y, cbind(param_mu, param_sigma)))


par(mar=c(4.1, 5, 4.1, 1.1))

qqnorm(resid_m1, main = expression('(a) Normalized quantile residuals (M1)'),
       cex.main=1,cex.lab=1, ylim = c(-4.3, 6.8), xlim = c(-4.3, 6.8),
       cex.axis=0.75)
abline(0,1)

qqnorm(resid_m2, main = expression('(b) Normalized quantile residuals (M2)'),
       cex.main=1,cex.lab=1, ylim = c(-4.3, 6.8), xlim = c(-4.3, 6.8),
       cex.axis=0.75)
abline(0,1)

Information Criteria and Log-Score

cv_log_score <- function(model, rerun = FALSE, n_folds = 10,
                         n_cores = ifelse(.Platform$OS.type == "windows", 1, n_folds)) {

  RNGkind("L'Ecuyer-CMRG")
  set.seed(42)

  model_name <- deparse(substitute(model))
  d <- model$model.frame
  n_rows <- nrow(d)
  rand_ind <- sample(n_rows)
  folds <- cut(seq_len(n_rows), breaks = n_folds, labels = FALSE)

  log_score_k <- function(k) {

    test_ind <- rand_ind[folds == k]

    if (rerun) {
      m_train <- update(model, data = d[-test_ind, ])
      saveRDS(m_train, sprintf("%s%s_train_fold_%d.rds", pathdata, model_name, k))
    } else {
      m_train <- readRDS(sprintf("%s%s_train_fold_%d.rds", pathdata, model_name, k))
    }

    d_test <- d[test_ind, ]
    pred <- predict(m_train, newdata = d_test[, -1, drop = FALSE], what = "samples",
                    type = "parameter", FUN = identity)
    pred_mu <- as.matrix(pred[["mu"]])
    pred_sigma <- as.matrix(pred[["sigma"]])

    ptw_pred_lik <-
      sapply(seq_along(test_ind),
             function(row) {
               mean(dnorm(d_test[row, 1], pred_mu[row, ], pred_sigma[row, ]))
             })

    sum(log(ptw_pred_lik))
  }

 log_score_folds <- lapply(seq_len(n_folds), log_score_k)
 
 sum(unlist(log_score_folds))
}

bamlss::DIC(m1, m2)

        DIC       pd
m1 82156.59 43.41082
m2 81128.05 84.69140

bamlss::WAIC(m1, m2)

      WAIC1    WAIC2        p1       p2
m1 82159.06 82159.46  45.87964  46.0834
m2 81151.72 81156.62 108.36408 110.8171

log_score_m1 <- cv_log_score(m1)
log_score_m2 <- cv_log_score(m2)
log_score_m1

[1] -41082.62

log_score_m2

[1] -40589.01

11.2.4 Regression Effects and Hyperprior Sensitivity

Hyperprior Sensitivity

Predict mu over age

pred_mu_func <- function(state = 33, csex = 0, model){
  nd <- data.frame("cage" = seq_range(data$cage,100))

  nd$state <- state
  nd$csex <- csex # men

  pred <- predict(model,newdata=nd,what="samples",FUN=function(x){x},intercept = TRUE)
  pred <- data.frame(pred["mu"])

  pred_mu <- data.frame("cage" = seq_range(data$cage,100))
  pred_mu$mean <- rowMeans(data.frame(pred))
  pred_mu$lq <- apply(data.frame(pred),  FUN="quantile",prob=0.025,MAR=1)
  pred_mu$uq <- apply(data.frame(pred),  FUN="quantile",prob=0.975,MAR=1)
  return(list(pred,pred_mu))
}

### CI; men

pred_mu_men <- pred_mu_func(model = m2)[[2]]

### CI; women

pred_mu_women <- pred_mu_func(model = m2, csex = 1)[[2]]

mu plot over age for all hyperparameters

plot_mu <- function(model , colour = TRUE, main, both = TRUE) {
  col <- ifelse(colour,"_col", "")
  col_ <- ifelse(colour,1, 0)

  pred_mu_women <- pred_mu_func(csex = 1, model = model)[[2]]
  pred_mu_men <- pred_mu_func(csex = 0, model = model)[[2]]

  main <- ifelse(substitute(model) == "m_ig",
                 TeX(r"((a) Predicted $\mu$ with $IG(0.001,0.001)$)"),
                 ifelse(substitute(model) == "m_sd",
                        TeX(r"((b) Predicted $\mu$ with $SD(0.009)$)"),
                        ifelse(substitute(model) == "m_hc",
                               TeX(r"((c) Predicted $\mu$ with $\it{C}^{+}(0.01)$)"),
                               TeX(r"((d) Predicted $\mu$ with $\it{U}(0.27)$)"))))

  plot(pred_mu_women$cage,pred_mu_women$mean,type="l",
       col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1, main = main,
       ylab=TeX(r"($\widehat{\mu}(Z_{wa}\,|\,x)$)"),xlab="child's age", ylim = c(-1.5, 0))
  lines(pred_mu_women$cage,pred_mu_women$uq,type="l",
        col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
  lines(pred_mu_women$cage,pred_mu_women$lq,type="l",
        col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
  polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
          c(pred_mu_women$mean,rev(pred_mu_women$uq)),border=NA,
          col=ifelse(colour,rgb(red = col_, green = 0, blue = 0,alpha = 0.5),
                     rgb(127,127,127,  maxColorValue = 255)))
  polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
          c(pred_mu_women$mean,rev(pred_mu_women$lq)),border=NA,
          col=ifelse(colour,rgb(red = col_, green = 0, blue = 0,alpha = 0.5),
                     rgb(127,127,127,  maxColorValue = 255)))
if (both) {
  lines(pred_mu_men$cage,pred_mu_men$uq,type="l",
        col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
  lines(pred_mu_men$cage,pred_mu_men$lq,type="l",
        col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
  polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
          c(pred_mu_men$mean,rev(pred_mu_men$uq)),border=NA,
          col=ifelse(colour, rgb(red = 0, green = 0, blue = col_,alpha = 0.3),
                     rgb(195,195,195,  maxColorValue = 255)))
  polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
          c(pred_mu_men$mean,rev(pred_mu_men$lq)),border=NA,
          col=ifelse(colour, rgb(red = 0, green = 0, blue = col_,alpha = 0.3),
                     rgb(195,195,195, maxColorValue = 255)))


  lines(pred_mu_men$cage,pred_mu_men$mean,type="l",
        col=rgb(red = 0, green = 0, blue = 0,alpha = 0.5))


  lines(pred_mu_women$cage,pred_mu_women$mean,type="l",
        col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
  legend("topright",cex = 1, legend=c("Girls", "Boys"),
         fill= c(rgb(red = col_, green = 0, blue = 0,alpha = 0.5),
                 rgb(red = 0, green = 0, blue = col_,alpha = 0.3)), bty="n")
}  else {
  lines(pred_mu_women$cage,pred_mu_women$mean,type="l",
          col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
}

}

par(mar=c(4.4, 5.6, 3.1, 1.1), mfrow = c(2,2))

plot_mu(model = m_ig , colour = T, main = "(a)")
plot_mu(model = m_sd , colour = T, main = "(b)")
plot_mu(model = m_hc , colour = T, main = "(c)")
plot_mu(model = m_u , colour = T, main = "(d)")

Regression Effects

Compute and plot simultaneous CI

#### simultaneous confidence bands; men -------------

pred <- pred_mu_func(model = m2)[[1]]

pred_mu_men$lq_sim <- pred_mu_men$lq
pred_mu_men$uq_sim <- pred_mu_men$uq

temp <- pred

n = 0

while((n/ncol(temp))<0.95){
  pred_mu_men$uq_sim <-
    pred_mu_men$mean + 1.001*(abs(pred_mu_men$uq_sim - pred_mu_men$mean))
  pred_mu_men$lq_sim <-
    pred_mu_men$mean - 1.001*(abs(pred_mu_men$lq_sim - pred_mu_men$mean))
  n = 0
  for (x in 1:ncol(temp)){
    if (all(temp[,x] <= pred_mu_men$uq_sim) &   all(temp[,x] >= pred_mu_men$lq_sim)){
      n <- n+1
    }
  }
}

#### simultaneous confidence bands; women -------------

pred <- pred_mu_func(model = m2,csex = 1)[[1]]

pred_mu_women$lq_sim <- pred_mu_women$lq
pred_mu_women$uq_sim <- pred_mu_women$uq

temp <- pred

n = 0

while((n/ncol(temp))<0.95){
  pred_mu_women$uq_sim <-
    pred_mu_women$mean + 1.001*(abs(pred_mu_women$uq_sim - pred_mu_women$mean))
  pred_mu_women$lq_sim <-
    pred_mu_women$mean - 1.001*(abs(pred_mu_women$lq_sim - pred_mu_women$mean))
  n = 0
  for (x in 1:ncol(temp)){
    if (all(temp[,x] <= pred_mu_women$uq_sim) &   all(temp[,x] >= pred_mu_women$lq_sim)){
      n <- n+1
    }
  }
}

par(mar=c(4.4, 5.6, 3.1, 1.1), mfrow = c(1,2))

plot(pred_mu_women$cage,pred_mu_women$mean,type="l",
     col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1,
     main = expression("(a) Predicted" ~ mu ~ "- Girls") ,
     ylab=TeX(r"($\widehat{\mu}(Z_{wa}\,|\,x)$)"),xlab="child's age", ylim = c(-1.5, 0))
lines(pred_mu_women$cage,pred_mu_women$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_women$cage,pred_mu_women$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_women$cage,pred_mu_women$uq_sim,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_women$cage,pred_mu_women$lq_sim,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
        c(pred_mu_women$mean,rev(pred_mu_women$uq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.5))
polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
        c(pred_mu_women$mean,rev(pred_mu_women$lq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.5))
polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
        c(pred_mu_women$uq_sim,rev(pred_mu_women$uq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.2))
polygon(c(pred_mu_women$cage,rev(pred_mu_women$cage)),
        c(pred_mu_women$lq_sim,rev(pred_mu_women$lq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.2))
legend("topright",cex = 1, legend=c("Pointwise CI", "Simultaneous CI"),
       fill= c(rgb(red = 1, green = 0, blue = 0,alpha = 0.5),
               rgb(red = 1, green = 0, blue = 0,alpha = 0.2)), bty="n")

plot(pred_mu_men$cage,pred_mu_men$mean,type="l",
     col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1,
     main = expression("(b) Predicted" ~ mu ~ "- Boys") ,
     ylab="",xlab="child's age", ylim = c(-1.5, 0))
lines(pred_mu_men$cage,pred_mu_men$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_men$cage,pred_mu_men$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_men$cage,pred_mu_men$uq_sim,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu_men$cage,pred_mu_men$lq_sim,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
        c(pred_mu_men$mean,rev(pred_mu_men$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.3))
polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
        c(pred_mu_men$mean,rev(pred_mu_men$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.3))
polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
        c(pred_mu_men$uq_sim,rev(pred_mu_men$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.1))
polygon(c(pred_mu_men$cage,rev(pred_mu_men$cage)),
        c(pred_mu_men$lq_sim,rev(pred_mu_men$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.1))
legend("topright",cex = 1, legend=c("Pointwise CI", "Simultaneous CI"),
       fill= c(rgb(red = 0, green = 0, blue = 1,alpha = 0.3),
               rgb(red = 0, green = 0, blue = 1,alpha = 0.1)), bty="n")

Predict sigma over age

pred_sigma_func <- function(state = 33, csex = 0, model){
  nd <- data.frame("cage" = seq_range(data$cage,100))

  nd$state <- state
  nd$csex <- csex # men

  pred <- predict(model,newdata=nd,what="samples",FUN=function(x){x},intercept = TRUE)
  pred <- data.frame(pred["sigma"])

  pred_sigma <- data.frame("cage" = seq_range(data$cage,100))
  pred_sigma$mean <- rowMeans(exp(data.frame(pred)))
  pred_sigma$lq <- apply(exp(pred),  FUN="quantile",prob=0.025,MAR=1)
  pred_sigma$uq <- apply(exp(pred),  FUN="quantile",prob=0.975,MAR=1)
  return(list(pred,pred_sigma))
}


######## men

pred_sigma_men <- pred_sigma_func(model = m2)[[2]]

######## women


pred_sigma_women <- pred_sigma_func(model = m2, csex = 1)[[2]]

Plot sigma over age

par(mar=c(4.4, 5.6, 3.1, 1.1))

plot(pred_sigma_women$cage,pred_sigma_women$mean,type="l",
     col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1,
     main =  expression("Predicted" ~ sigma),
     ylab=TeX(r"($\widehat{\sigma}(Z_{wa}\,|\,x)$)"),xlab="child's age",
     ylim = c(0.8,1.4))
lines(pred_sigma_women$cage,pred_sigma_women$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_sigma_women$cage,pred_sigma_women$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_sigma_women$cage,rev(pred_sigma_women$cage)),
        c(pred_sigma_women$mean,rev(pred_sigma_women$uq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.5))
polygon(c(pred_sigma_women$cage,rev(pred_sigma_women$cage)),
        c(pred_sigma_women$mean,rev(pred_sigma_women$lq)),border=NA,
        col=rgb(red = 1, green = 0, blue = 0,alpha = 0.5))

lines(pred_sigma_men$cage,pred_sigma_men$mean,type="l",
      col=rgb(red = 0, green = 0, blue = 0))
lines(pred_sigma_men$cage,pred_sigma_men$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_sigma_men$cage,pred_sigma_men$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_sigma_men$cage,rev(pred_sigma_men$cage)),
        c(pred_sigma_men$mean,rev(pred_sigma_men$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.5))
polygon(c(pred_sigma_men$cage,rev(pred_sigma_men$cage)),
        c(pred_sigma_men$mean,rev(pred_sigma_men$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.5))
legend("topright",cex = 0.75, legend=c("Girls", "Boys"),
       fill= c(rgb(red = 1, green = 0, blue = 0,alpha = 0.5),
               rgb(red = 0, green = 0, blue = 1,alpha = 0.5)), bty="n")

Reproducible data analysis

library(parallel)
library(BayesX)
library(bamlss)
library(modelr)
library(latex2exp)
library(gamboostLSS)

# load data
data <- gamboostLSS::india

11.2 Univariate Analysis

11.2.1 Model Estimation

set.seed(42)

# model with constant variance ----------

f1 <- stunting ~ s(cage, xt = list(prior = "ig", a = 0.001, b = 0.001)) 

# Estimate model 

m1 <- bamlss(f1, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

# model with varying variance ----------

f2 <- list(stunting ~  s(cage, xt = list(prior = "ig", a = 0.001, b = 0.001)),
           sigma ~  s(cage))

# Estimate model 

m2 <- bamlss(f2, gaussian, data = data, n.iter = 12000, burnin = 2000, thin = 10)

11.2.2 MCMC Diagnostics

Effective sample size

samples <- data.frame(m1[["samples"]][[1]])$mu.s.s.cage..b2
samples <- as.mcmc(samples)
effectiveSize(samples) 

var1 
1001 

Trajectory

par(mar=c(4.1, 6.3, 4.1, 1.1))
plot(1:1001,data.frame(m1[["samples"]][[1]])$mu.s.s.cage..b2, type = "l",
     main = TeX(r"((a) MCMC sampling path of $\beta^{\mu}_{{age}_{2}}$)"), 
     ylab = TeX(r"($\beta^{\mu}_{{spat}_{33}}$)"), 
     xlab = "Iterations", cex.main = 1, cex.lab = 1, cex.axis = 0.75)

ACF

par(mar=c(4.5, 5, 4.1, 1.1))
autocorr.plot(samples, auto.layout = T, lwd = 3, 
              main = TeX(r"((b) ACF plot of $\beta^{\mu}_{{age}_{2}}$)"),
              cex.main = 1, cex.lab = 1, cex.axis = 0.75)

11.2.3 Model Comparison

Normalized Quantile Residuals

##### m1 --------------

pred <- predict(m1,newdata=data,what="samples",type = "parameter", FUN=identity)

param_sigma <- rowMeans(data.frame(pred["sigma"]))
param_mu <- rowMeans(data.frame(pred["mu"]) )

y <- data$stunting

pn <- function (y,  param)
{
  mean <-  param[, 1]
  sd <- (param[, 2])
  pnorm(y, mean = mean, sd = sd)
}


resid_m1=qnorm(pn(y, cbind(param_mu, param_sigma)))

##### m2 --------------
pred <- predict(m2,newdata=data,what="samples",type = "parameter", FUN=identity)

param_sigma <- rowMeans(data.frame(pred["sigma"]))
param_mu <- rowMeans(data.frame(pred["mu"]) )

y <- data$stunting

pn <- function (y,  param)
{
  mean <-  param[, 1]
  sd <- (param[, 2])
  pnorm(y, mean = mean, sd = sd)
}

resid_m2=qnorm(pn(y, cbind(param_mu, param_sigma)))

par(mar=c(4.1, 5, 4.1, 1.1))

qqnorm(resid_m1, main = expression('(a) Normalized quantile residuals (M1)'),
       cex.main=1,cex.lab=1, ylim = c(-4.3, 5.5), xlim = c(-4.3, 5.5),
       cex.axis=0.75)
abline(0,1)

qqnorm(resid_m2, main = expression('(b) Normalized quantile residuals (M2)'),
       cex.main=1,cex.lab=1, ylim = c(-4.3, 5.5), xlim = c(-4.3, 5.5),
       cex.axis=0.75)
abline(0,1)

Information Criteria

bamlss::DIC(m1, m2)

        DIC        pd
m1 14820.27  8.731034
m2 14801.21 10.902704

bamlss::WAIC(m1, m2)

      WAIC1    WAIC2       p1        p2
m1 14820.55 14820.62  9.01215  9.048768
m2 14802.22 14802.47 11.91700 12.042789

11.2.4 Regression Effects

pred_mu_func <- function(model){
  nd <- data.frame("cage" = seq_range(data$cage,100))

  pred <- predict(model,newdata=nd,what="samples",FUN=function(x){x},intercept = TRUE)
  pred <- data.frame(pred["mu"])

  pred_mu <- data.frame("cage" = seq_range(data$cage,100))
  pred_mu$mean <- rowMeans(data.frame(pred))
  pred_mu$lq <- apply(data.frame(pred),  FUN="quantile",prob=0.025,MAR=1)
  pred_mu$uq <- apply(data.frame(pred),  FUN="quantile",prob=0.975,MAR=1)
  return(list(pred,pred_mu))
}

Compute and plot simultaneous CI

pred_mu <- pred_mu_func(model = m2)[[2]]
pred <- pred_mu_func(model = m2)[[1]]

pred_mu$lq_sim <- pred_mu$lq
pred_mu$uq_sim <- pred_mu$uq

temp <- pred

n = 0

while((n/ncol(temp))<0.95){
  pred_mu$uq_sim <-
    pred_mu$mean + 1.001*(abs(pred_mu$uq_sim - pred_mu$mean))
  pred_mu$lq_sim <-
    pred_mu$mean - 1.001*(abs(pred_mu$lq_sim - pred_mu$mean))
  n = 0
  for (x in 1:ncol(temp)){
    if (all(temp[,x] <= pred_mu$uq_sim) &   all(temp[,x] >= pred_mu$lq_sim)){
      n <- n+1
    }
  }
}

par(mar=c(4.4, 5.6, 3.1, 1.1))

plot(pred_mu$cage,pred_mu$mean,type="l",
     col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1,
     main = expression("Predicted" ~ mu) ,
     ylab=TeX(r"($\widehat{\mu}(Z_{wa}\,|\,x)$)"),xlab="child's age", ylim = c(-3, 1))
lines(pred_mu$cage,pred_mu$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu$cage,pred_mu$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu$cage,pred_mu$uq_sim,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_mu$cage,pred_mu$lq_sim,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_mu$cage,rev(pred_mu$cage)),
        c(pred_mu$mean,rev(pred_mu$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.3))
polygon(c(pred_mu$cage,rev(pred_mu$cage)),
        c(pred_mu$mean,rev(pred_mu$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.3))
polygon(c(pred_mu$cage,rev(pred_mu$cage)),
        c(pred_mu$uq_sim,rev(pred_mu$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.1))
polygon(c(pred_mu$cage,rev(pred_mu$cage)),
        c(pred_mu$lq_sim,rev(pred_mu$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.1))
legend("topright",cex = 0.75, legend=c("Pointwise CI", "Simultaneous CI"),
       fill= c(rgb(red = 0, green = 0, blue = 1,alpha = 0.3),
               rgb(red = 0, green = 0, blue = 1,alpha = 0.1)), bty="n")

Predict sigma over age

pred_sigma_func <- function(model){
  nd <- data.frame("cage" = seq_range(data$cage,100))

  pred <- predict(model,newdata=nd,what="samples",FUN=function(x){x},intercept = TRUE)
  pred <- data.frame(pred["sigma"])

  pred_sigma <- data.frame("cage" = seq_range(data$cage,100))
  pred_sigma$mean <- rowMeans(exp(data.frame(pred)))
  pred_sigma$lq <- apply(exp(pred),  FUN="quantile",prob=0.025,MAR=1)
  pred_sigma$uq <- apply(exp(pred),  FUN="quantile",prob=0.975,MAR=1)
  return(list(pred,pred_sigma))
}

pred_sigma <- pred_sigma_func(model = m2)[[2]]

Plot sigma over age

par(mar=c(4.4, 5.6, 3.1, 1.1))

plot(pred_sigma$cage,pred_sigma$mean,type="l",
     col=rgb(red = 0, green = 0, blue = 0),cex.lab=1,cex.main=1,cex.axis=1,
     main =  expression("Predicted" ~ sigma),
     ylab=TeX(r"($\widehat{\sigma}(Z_{wa}\,|\,x)$)"),xlab="child's age",
     ylim = c(1.2, 1.8))
lines(pred_sigma$cage,pred_sigma$uq,type="l",
      col=rgb(red =0, green = 0, blue = 0,alpha = 0.7))
lines(pred_sigma$cage,pred_sigma$lq,type="l",
      col=rgb(red = 0, green = 0, blue = 0,alpha = 0.7))
polygon(c(pred_sigma$cage,rev(pred_sigma$cage)),
        c(pred_sigma$mean,rev(pred_sigma$uq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.5))
polygon(c(pred_sigma$cage,rev(pred_sigma$cage)),
        c(pred_sigma$mean,rev(pred_sigma$lq)),border=NA,
        col=rgb(red = 0, green = 0, blue = 1,alpha = 0.5))