Question about hierarchical effects

lizzieinvancouver · February 8, 2018, 1:33am

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)

bgoodri · February 8, 2018, 2:13am

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.

lizzieinvancouver · February 15, 2018, 4:17am

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();

Topic		Replies	Views
New group level in RStanarm output rstanarm	2	191	March 6, 2024
Intercept only model rstanarm rstan	7	1630	August 26, 2020
Rstanarm: ranef vs b parameters rstanarm hierarchical-model , interpret-results	1	630	March 23, 2021
How correlation is computed in rstanarm output rstanarm	1	528	February 24, 2020
Can I obtain bayes_R2 for reach level using rstanarm? rstanarm techniques	7	1317	January 26, 2018

Question about hierarchical effects

Related topics