Question about hierarchical effects

Interfaces rstanarm
1

Operating System: Mac OS X
Interface Version: 2.17.2
Compiler/Toolkit:

I am trying to run a bernoulli model in rstanarm and plot the overall estimates versus each group estimate (see plot attached). I am adapting some code that I have used for rstan output.

What I am thinking are the overall estimates, however, are not lining up with the group-level estimates. I suspect there is something fundamental here I am missing about rstanarm (though it could be an R coding issue I have missed despite my best efforts…).

Code and data attached (I believe). Any help much appreciated.

L

PS: I am new to using rstanarm and discourse so apologies all around for whatever basics I am screwing up here.

plot.pdf (110.4 KB)
Analyzing non-leafouts POST.R (2.1 KB)
dx.csv (265.8 KB)

2

I think you are looking at the posterior distribution of the b parameters — which are deviations from the global parameters (in the non-centered parameterization) — rather than the sum of b and the global parameter, which represents the intercept for a group (in the centered parameterization). If you were just looking at posterior medians, it is the difference between the output of the ranef function and the coef function. But you don’t want to just be looking at posterior medians usually. You can get the raw draws of each parameter by doing as.data.frame or as.matrix, and then you can add things together to get the posterior distribution of the group-level intercepts.

3

Yes, that was it. Thanks! Updated figures and code attached/below …

L

NonBB_sp.pdf (30.2 KB)

iter.m1no <- as.data.frame(m1.no)
iter.m1nl <- as.data.frame(m1.nl)

# manually to get right order, with intercept
mu_params.wi <- c("(Intercept)", "warm20","photo12","chillchill1","chillchill2",
               "siteSH","warm20:photo12",
               "warm20:chillchill1","warm20:chillchill2",
               "photo12:chillchill1","photo12:chillchill2",
               "warm20:siteSH", "photo12:siteSH",
               "chillchill1:siteSH","chillchill2:siteSH")

meanzb.wi <- sumer.m1no[mu_params.wi,col4fig]

rownames(meanzb.wi) = c("Intercept",
                    "Forcing Temperature",
                    "Photoperiod",
                    "Chilling 4°",
                    "Chilling 1.5°C",
                    "Site",
                    "Forcing x Photoperiod",
                    "Forcing x Chilling 4°C",
                    "Forcing x Chilling 1.5°C",
                    "Photoperiod x Chilling 4°C",
                    "Photoperiod x Chilling 1.5°C",
                    "Forcing x Site",
                    "Photoperiod x Site",
                    "Site x Chilling 4°C",
                    "Site x Chilling 1.5°C"
                    )

meanzl.wi <- sumer.m1nl[mu_params.wi,col4fig]
rownames(meanzl.wi) <- rownames(meanzb.wi)

speff.bb <- speff.lo <- vector()
params <- c("(Intercept)", "warm20","photo12","chillchill1",
               "chillchill2","siteSH", "warm20:photo12",
               "warm20:chillchill1","warm20:chillchill2",
               "photo12:chillchill1","photo12:chillchill2",
               "warm20:siteSH", "photo12:siteSH",
               "chillchill1:siteSH","chillchill2:siteSH")

sp.params <- c("(Intercept)", "chillchill1","chillchill2","siteSH",
               "chillchill1:siteSH","chillchill2:siteSH")

params.wsp <- c(1, 4:6, 14:15)
params.nosp <- c(1:15)[-params.wsp]

pdf(file.path(figpath, "NonBB_sp.pdf"), width = 7, height = 8)

par(mfrow=c(1,1), mar = c(2, 10, 2, 1))
# Upper panel: budburst
plot(seq(-4, #min(meanz[,'mean']*1.1),
         5, #max(meanz[,'mean']*1.1),
         length.out = nrow(meanzb.wi)), 
     seq(1, 5*nrow(meanzb.wi), length.out = nrow(meanzb.wi)),
     type="n",
     xlab = "",
     ylab = "",
     yaxt = "n")

legend(x =-4.5, y = 6, bty="n", legend = "a. Budburst", text.font = 2)
# rasterImage(bbpng, -0.25, 1, 0, 4)

axis(2, at = 5*(nrow(meanzb.wi):1), labels = rownames(meanzb.wi), las = 1, cex.axis = 0.8)


# Plot species levels for each predictor
for(i in 1:length(unique(dx$sp))){
  b.params <- iter.m1no[!is.na(match(colnames(iter.m1no), c(paste("b", "[", sp.params, " sp:",
      unique(dx$sp)[i], "]", sep=""))))]

  main.params <- iter.m1no[!is.na(match(colnames(iter.m1no), sp.params))]

  bplusmain <- b.params
  for(c in 1:ncol(main.params)){
      bplusmain[c] <- b.params[c]+main.params[c]
      }

  bplusmain.quant <- sapply(bplusmain, FUN = quantile, probs = c(0.25, 0.50, 0.75))
  
  sp.est <- t(bplusmain.quant)
  
  jt <- jitter(0, factor = 40)

  arrows(sp.est[,"75%"],  jt+(5*(nrow(meanzb.wi):1)-1)[params.wsp], sp.est[,"25%"],  jt+(5*(nrow(meanzb.wi):1)-1)[params.wsp],
         len = 0, col = alpha("firebrick", 0.2)) 
  
  points(sp.est[,'50%'],
         jt+(5*(nrow(meanzb.wi):1)-1)[params.wsp], #[c(3:5,11:12)], # ADJUSTED for just the ranef here
         pch = 16,
         col = alpha("firebrick", 0.5))

  speff.bb = rbind(speff.bb, t(sp.est[,1]))
    }

arrows(meanzb.wi[,"75%"], (5*(nrow(meanzb.wi):1))+1, meanzb.wi[,"25%"], (5*(nrow(meanzb.wi):1))+1,
       len = 0, col = "black", lwd = 3)

points(meanzb.wi[,'mean'],
       (5*(nrow(meanzb.wi):1))+1,
       pch = 16,
       cex = 1,
       col = "midnightblue")
abline(v = 0, lty = 2)


# Lower panel: leafout

par(mfrow=c(1,1), mar = c(2, 10, 2, 1))

plot(seq(-4, #min(meanz[,'mean']*1.1),
         5, #max(meanz[,'mean']*1.1),
         length.out = nrow(meanzl.wi)), 
     seq(1, 5*nrow(meanzl.wi), length.out = nrow(meanzl.wi)),
     type="n",
     xlab = "",
     ylab = "",
     yaxt = "n")

legend(x =-4.5, y = 6, bty="n", legend = "b. Leafout", text.font = 2)

axis(2, at = 5*(nrow(meanzl.wi):1), labels = rownames(meanzl.wi), las = 1, cex.axis = 0.8)


# Plot species levels for each predictor
for(i in 1:length(unique(dx$sp))){
  b.params <- iter.m1nl[!is.na(match(colnames(iter.m1nl), c(paste("b", "[", sp.params, " sp:",
      unique(dx$sp)[i], "]", sep=""))))]

  main.params <- iter.m1nl[!is.na(match(colnames(iter.m1nl), sp.params))]

  bplusmain <- b.params
  for(c in 1:ncol(main.params)){
      bplusmain[c] <- b.params[c]+main.params[c]
      }

  bplusmain.quant <- sapply(bplusmain, FUN = quantile, probs = c(0.25, 0.50, 0.75))
  
  sp.est <- t(bplusmain.quant)
  
  jt <- jitter(0, factor = 40)

  arrows(sp.est[,"75%"],  jt+(5*(nrow(meanzl.wi):1)-1)[params.wsp], sp.est[,"25%"],  jt+(5*(nrow(meanzl.wi):1)-1)[params.wsp],
         len = 0, col = alpha("firebrick", 0.2)) 
  
  points(sp.est[,'50%'],
         jt+(5*(nrow(meanzl.wi):1)-1)[params.wsp], #[c(3:5,11:12)], # ADJUSTED for just the ranef here
         pch = 16,
         col = alpha("firebrick", 0.5))

  speff.lo = rbind(speff.lo, t(sp.est[,1]))
    }

arrows(meanzl.wi[,"75%"], (5*(nrow(meanzl.wi):1))+1, meanzl.wi[,"25%"], (5*(nrow(meanzl.wi):1))+1,
       len = 0, col = "black", lwd = 3)

points(meanzl.wi[,'mean'],
       (5*(nrow(meanzl.wi):1))+1,
       pch = 16,
       cex = 1,
       col = "midnightblue")
abline(v = 0, lty = 2)


dev.off();