---
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.