########################################################### ######## STAN models and plotting for Pierson et al. ###### ########################################################### ########################################################### ###################### Manipulate data #################### ########################################################### Sys.setenv(TZ='EST') # Load packages. library(gdata) library(shinystan) library(rstan) library(loo) library(ggplot2) library(gridExtra) library(smoothr) library(rstanarm) library(bayesplot) # Read in data. plethDat <- read.xls("Plethodon_Data.xlsx") # Pull out all weights. weight1 <- plethDat$Weight1 weight2 <- plethDat$Weight2 weight3 <- plethDat$Weight3 weight4 <- plethDat$Weight4 weight5 <- plethDat$Weight5 weight6 <- plethDat$Weight6 weight7 <- plethDat$Weight7 # Average weights across measurements 4-7. newWeight <- NA for (i in 1:length(weight4)) { newWeight[i] <- mean(weight4[i],weight5[i], weight6[i],weight7[i]) } # Pull out and average all SVLs. svl1 <- vector(length=length(plethDat$Tag)) svl2 <- vector(length=length(plethDat$Tag)) svl3 <- vector(length=length(plethDat$Tag)) svl4 <- vector(length=length(plethDat$Tag)) svl5 <- vector(length=length(plethDat$Tag)) svl6 <- vector(length=length(plethDat$Tag)) svl7 <- vector(length=length(plethDat$Tag)) for (i in 1:length(plethDat$Tag)){ svl1[i] <- mean(plethDat$SVL1A[i],plethDat$SVL1B[i]) svl2[i] <- mean(plethDat$SVL2A[i],plethDat$SVL2B[i]) svl3[i] <- mean(plethDat$SVL3A[i],plethDat$SVL3B[i]) svl4[i] <- mean(plethDat$SVL4A[i],plethDat$SVL4B[i]) svl5[i] <- mean(plethDat$SVL5A[i],plethDat$SVL5B[i]) svl6[i] <- mean(plethDat$SVL6A[i],plethDat$SVL6B[i]) svl7[i] <- mean(plethDat$SVL7A[i],plethDat$SVL7B[i]) } # Average SVL across measurements 4-7. newSVL <- NA for (i in 1:length(svl4)) { newSVL[i] <- mean(svl4[i],svl5[i], svl6[i],svl7[i]) } # Pull out and average all HWs. hw1 <- vector(length=length(plethDat$Tag)) hw2 <- vector(length=length(plethDat$Tag)) hw3 <- vector(length=length(plethDat$Tag)) hw4 <- vector(length=length(plethDat$Tag)) hw5 <- vector(length=length(plethDat$Tag)) hw6 <- vector(length=length(plethDat$Tag)) hw7 <- vector(length=length(plethDat$Tag)) for (i in 1:length(plethDat$Tag)){ hw1[i] <- mean(plethDat$HW1A[i],plethDat$HW1B[i]) hw2[i] <- mean(plethDat$HW2A[i],plethDat$HW2B[i]) hw3[i] <- mean(plethDat$HW3A[i],plethDat$HW3B[i]) hw4[i] <- mean(plethDat$HW4A[i],plethDat$HW4B[i]) hw5[i] <- mean(plethDat$HW5A[i],plethDat$HW5B[i]) hw6[i] <- mean(plethDat$HW6A[i],plethDat$HW6B[i]) hw7[i] <- mean(plethDat$HW7A[i],plethDat$HW7B[i]) } # Average HW across measurements 3-7. newHW <- NA for (i in 1:length(hw4)) { newHW[i] <- mean(hw4[i],hw5[i], hw6[i],hw7[i]) } # Pull out treatments and calculate changes. treat <- plethDat$Treatment weightChange <- (weight7-weight1)/weight1 SVLchange <- (svl7-svl1)/svl1 HWchange <- (hw7-hw1)/hw1 ratioChange <- ((svl7/hw7)-(svl1/hw1))/(svl1/hw1) ########################################################### ############ Assess interobserver reliability. ############ ########################################################### # For SVL. allSVLa <- c(plethDat$SVL1A,plethDat$SVL2A,plethDat$SVL3A, plethDat$SVL4A,plethDat$SVL5A,plethDat$SVL6A, plethDat$SVL7A) allSVLb <- c(plethDat$SVL1B,plethDat$SVL2B,plethDat$SVL3B, plethDat$SVL4B,plethDat$SVL5B,plethDat$SVL6B, plethDat$SVL7A) # Measure correlation. cor(allSVLa,allSVLb) # For head width. allHWa <- c(plethDat$HW1A,plethDat$HW2A,plethDat$HW3A, plethDat$HW4A,plethDat$HW5A,plethDat$HW6A, plethDat$HW7A) allHWb <- c(plethDat$HW1B,plethDat$HW2B,plethDat$HW3B, plethDat$HW4B,plethDat$HW5B,plethDat$HW6B, plethDat$HW7A) # Measure correlation. cor(allHWa,allHWb) ########################################################### ################ Build and run STAN models ################ ########################################################### # Settings for STAN. rstan_options(auto_write=TRUE) options(mc.cores = parallel::detectCores()) # Set colors and symbols. col <- alpha(c("red","blue","green","orange","purple"),0.01) col2 <- c("red","blue","green","orange","purple") symbs <- c(15,17,18,19,20) # Full model. newSVL <- runif(1000,min=25,max=75) fullMod <- list(nObs=length(weightChange), nTreats=length(unique(treat)), weightChange=weightChange, SVL=svl1, treat=treat, HWchange = HWchange, SVLchange = SVLchange, ratioChange = ratioChange ) fullMod <- stan(file="fullMod.stan", dat=fullMod, iter=1000, chains=2, seed=4) plot(fullMod, pars="alpha[1]") plot(fullMod) fullModSummary <- summary(fullMod) postFullMod <- as.data.frame(fullMod, pars=c("alpha","beta1","beta2","beta3", "rho1", "rho2") ) postFullMod <-as.matrix(postFullMod) posterior_interval(postFullMod, prob = 0.80, type = "central", pars = c("alpha","beta1","beta2","beta3", "rho1", "rho2") ) posterior_interval(postFullMod, prob = 0.95, type = "central", pars = c("alpha","beta1","beta2","beta3", "rho1", "rho2") ) # Calculate changes. posterior <- as.matrix(fullMod) SVLchange1 <- posterior[,"alpha[1]"]*posterior[,"rho1"] SVLchange2 <- posterior[,"alpha[2]"]*posterior[,"rho1"] SVLchange3 <- posterior[,"alpha[3]"]*posterior[,"rho1"] SVLchange4 <- posterior[,"alpha[4]"]*posterior[,"rho1"] SVLchange5 <- posterior[,"alpha[5]"]*posterior[,"rho1"] HWchange1 <- posterior[,"alpha[1]"]*posterior[,"rho2"] HWchange2 <- posterior[,"alpha[2]"]*posterior[,"rho2"] HWchange3 <- posterior[,"alpha[3]"]*posterior[,"rho2"] HWchange4 <- posterior[,"alpha[4]"]*posterior[,"rho2"] HWchange5 <- posterior[,"alpha[5]"]*posterior[,"rho2"] ratioChange1 <- (1+SVLchange1)/(1+HWchange1)-1 ratioChange2 <- (1+SVLchange2)/(1+HWchange2)-1 ratioChange3 <- (1+SVLchange3)/(1+HWchange3)-1 ratioChange4 <- (1+SVLchange4)/(1+HWchange4)-1 ratioChange5 <- (1+SVLchange5)/(1+HWchange5)-1