---
title: "R Notebook"
output: html_notebook
---
```{r}
rm(list = ls())
require(glmulti)
require(arm)
require(rsq)
df <- read.csv("~/Dropbox/Active/Arizona/1) CA Smith Fellowship/Manuscript/Datasets/DataS2.csv", skip=1)
df <- df[which(df$Invertebrate.richness>0),]
```
###Wrapper function from developer to only allow by basin interactions: https://vcalcagnoresearch.wordpress.com/package-glmulti/
```{r}
nullos <- glm(Invertebrate.richness~1,data=df)
myglm=function(y, data) {
termz = terms(y,data=data)
index=which(dimnames(attr(termz,"factors"))[[2]] == "Days.flowing:Drying.frequency" | dimnames(attr(termz,"factors"))[[2]] == "Drying.frequency:Days.flowing" | dimnames(attr(termz,"factors"))[[2]] == "Days.flowing:Flow.permanence" | dimnames(attr(termz,"factors"))[[2]] == "Flow.permanence:Days.flowing" | dimnames(attr(termz,"factors"))[[2]] == "Days.flowing:Stream.distance.to.perennial.refuge..m." | dimnames(attr(termz,"factors"))[[2]] == "Stream.distance.to.perennial.refuge..m.:Days.flowing" | dimnames(attr(termz,"factors"))[[2]] == "Days.flowing:Direction.to.perennial.water" | dimnames(attr(termz,"factors"))[[2]] == "Direction.to.perennial.water:Days.flowing" | dimnames(attr(termz,"factors"))[[2]] == "Drying.frequency:Flow.permanence" | dimnames(attr(termz,"factors"))[[2]] == "Flow.permanence:Drying.frequency" | dimnames(attr(termz,"factors"))[[2]] == "Drying.frequency:Stream.distance.to.perennial.refuge..m." | dimnames(attr(termz,"factors"))[[2]] == "Stream.distance.to.perennial.refuge..m.:Drying.frequency" | dimnames(attr(termz,"factors"))[[2]] == "Drying.frequency:Direction.to.perennial.water" | dimnames(attr(termz,"factors"))[[2]] == "Direction.to.perennial.water:Drying.frequency" | dimnames(attr(termz,"factors"))[[2]] == "Flow.permanence:Stream.distance.to.perennial.refuge..m." | dimnames(attr(termz,"factors"))[[2]] == "Stream.distance.to.perennial.refuge..m.:Flow.permanence" | dimnames(attr(termz,"factors"))[[2]] == "Flow.permanence:Direction.to.perennial.water" | dimnames(attr(termz,"factors"))[[2]] == "Direction.to.perennial.water:Flow.permanence" | dimnames(attr(termz,"factors"))[[2]] == "Stream.distance.to.perennial.refuge..m.:Direction.to.perennial.water" | dimnames(attr(termz,"factors"))[[2]] == "Direction.to.perennial.water:Stream.distance.to.perennial.refuge..m.")
if (length(index)==0) { # the desired effect must be present
return(glm(formula=y, data=data))
} else return(nullos);
}
```
###Scale all variables
```{r}
df$Days.flowing <- rescale(df$Days.flowing)
df$Drying.frequency <- rescale(df$Drying.frequency)
df$Flow.permanence <- rescale(df$Flow.permanence)
df$Stream.distance.to.perennial.refuge..m. <- rescale(df$Stream.distance.to.perennial.refuge..m.)
df$Protected.area <- as.factor(df$Basin)
df$Direction.to.perennial.water <- as.factor(df$Direction.to.perennial.water)
```
###Get all models and see which are set to null by wrapper function
```{r}
nmods <- glmulti(Invertebrate.richness ~ Drying.frequency + Days.flowing + Stream.distance.to.perennial.refuge..m. + Direction.to.perennial.water + Protected.area, data=df, fitfunc=myglm, marginality=T, crit = aicc, method = "d", level=2, report =F, plotty=F)
res <- glmulti(Invertebrate.richness ~ Drying.frequency + Days.flowing + Stream.distance.to.perennial.refuge..m. + Direction.to.perennial.water + Protected.area, data=df, fitfunc=myglm, marginality=T, crit = aicc, method = "h", level=2, report =F, plotty=F, confsetsize=nmods)
write(res, file = "AllModels.csv", ncolumns = if (is.character(x)) 1 else 5, append = FALSE, sep = ",")
```
###Take top #97 models (full set with interactions including basin)
```{r}
res2 <- glmulti(Invertebrate.richness ~ Drying.frequency + Days.flowing + Stream.distance.to.perennial.refuge..m. + Direction.to.perennial.water + Protected.area, data=df, fitfunc=myglm, marginality=T, crit = aicc, method = "h", level=2, report =F, plotty=F, confsetsize=97)
write(res2, file = "ModelSet.csv", append = FALSE, sep = ",")
print(res2)
```
```{r}
top <- weightable(res2)
write.csv(top, "ModelWeights.csv")
```
```{r}
coefs <- coef.glmulti(res2, select=2)
write.csv(coefs, "Coefficients2ICAverage.csv")
```
```{r}
pdf("ModelICs.pdf")
plot(res2)
dev.off()
```
```{r}
res13 <- glmulti(Invertebrate.richness ~ Drying.frequency + Days.flowing + Stream.distance.to.perennial.refuge..m. + Direction.to.perennial.water + Protected.area, data=df, fitfunc=myglm, marginality=T, crit = aicc, method = "h", level=2, report =F, plotty=F, confsetsize=13)
```
```{r}
pdf("Importance.pdf")
plot(res13, type="s")
dev.off()
```
```{r}
res@formulas[[1]]
best <- glm(Invertebrate.richness ~ 1 + Basin + Drying.frequency + Days.flowing +
Stream.distance.to.perennial.refuge..m. + Protected.area:Days.flowing, df, family = gaussian)
rsq(best, adj=TRUE)
```