--- title: "R Notebook" output: html_notebook --- This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. ```{r} #### Import data library(readxl) Pois <- read_excel("~/*drivename*/Young_Data/Model3.xlsx") attach(Pois) str(Pois) ``` ```{r} ##### Load packages library("car") library("effects") library("lme4") library("brms") library("dplyr") #install.packages("devtools") #devtools::install_github("mvuorre/brmstools") library ("brmstools") ``` ```{r} #######Standardize predictors rescale <- function (x, binary.inputs){ # function to rescale by subtracting the mean and dividing by 2 sd's x.obs <- x[!is.na(x)] if (!is.numeric(x)) x <- as.numeric(factor(x)) if (length(unique(x.obs))==2){ x <- (x-min(x.obs))/(max(x.obs)-min(x.obs)) if (binary.inputs=="0/1") return (x) else if (binary.inputs=="-0.5,0.5") return (x-0.5) else if (binary.inputs=="center") return (x-mean(x.obs)) else if (binary.inputs=="full") return ((x-mean(x.obs))/(2*sd(x.obs))) } else { return ((x-mean(x.obs))/(2*sd(x.obs))) } } Pois$ZTAve =as.vector(rescale(AveTemp)) Pois$ZFood =as.vector(rescale(NDVI)) Pois$ZWater =as.vector(rescale(Water)) ``` ```{r} ###### check VIFs < 4 v1 = vif(glm(Deaths ~ ZTAve+ZFood +ZWater, family = "poisson", data = Pois)) v2 = vif(glm(Deaths ~ ZTAve*ZWater, family = "poisson", data = Pois)) v3 = vif(glm(Deaths ~ ZFood+ZWater, family = "poisson", data = Pois)) v1 v2 v3 ``` ```{r} ##########run models p1 <- brm( Deaths ~ ZWater+ZTAve + (1|Year), family = "poisson", data=Pois , chains = 4, cores = 4, iter = 2000) p2 <- brm( Deaths ~ ZFood+ZTAve +(1|Year), family = "poisson", data=Pois , chains = 4, cores = 4, iter = 2000) ``` ```{r} ###### get summary of models summary(p1) summary(p2) p1.w<-waic(p1) p2.w<-waic(p2) p1.w p2.w ``` ```{r} ####### Postirior checks # first p1 #look at trace plots plot(p1, pars= parnames(p1)[1:10], ask=F) #### Density Check pp_check(p1 ,type = "stat_grouped", group = "ID", stat = "mean") pp_check(p1 ,type = "stat_grouped", group = "ID", stat = "var") pp.plot.dens.1<-bayesplot::pp_check(p1, nsamples = 100) pp.plot.dens.1 # now p2 #look at trace plots plot(p1, pars= parnames(p2)[1:3], ask=F) #### Density Check pp_check(p2 ,type = "stat_grouped", group = "ID", stat = "mean") pp_check(p2 ,type = "stat_grouped", group = "ID", stat = "var") pp.plot.dens.1<-bayesplot::pp_check(p2, nsamples = 100) pp.plot.dens.1 ``` ```{r} ##### look at the r-squareds cond.r2<-bayes_R2(p1, re_formula = NULL, summary = T) marginal.r2<-bayes_R2(p1, re_formula = NA, summary = T) cond.minus.Slopes.r2<-bayes_R2(p1, re_formula = ~(1|Year), summary = T) cond.r2p<-bayes_R2(p2, re_formula = NULL, summary = T) marginal.r2p<-bayes_R2(p2, re_formula = NA, summary = T) cond.minus.Slopes.r2p<-bayes_R2(p2, re_formula = ~(1|Year), summary = T) #output p1 marginal.r2 cond.minus.Slopes.r2 cond.r2 # output p2 marginal.r2p cond.minus.Slopes.r2p cond.r2p ``` Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*. When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.