--- 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) ```