In brouwern/DRmencia: Data and code to replicate analyses for Latta, Brouwer, et al. 2018. PeerJ.

knitr::opts_chunk$set(echo = TRUE)

Note: this uses lm() b/c of limited sample size

Libraries

library(here)
fi <- here::here("R","load_libraries.R")
source(fi)
load_libraries()

library(wrapr)

Load condition data

load(file = "./data/condition_Mhat.RData")

any(ls(pattern = "condition.Mhat") == "condition.Mhat")

Remove bad entries

Remove bad entries (very rare)

i.mig <- which(condition.Mhat$stat.focals == "mig" &
              is.na(condition.Mhat$sex2) == FALSE)
i.res <- which(condition.Mhat$stat.focals == "res" &
              is.na(condition.Mhat$sex2) == FALSE)

i.all <- c(i.mig,i.res)

i.BANA.bad <- with(condition.Mhat[,],
     which(spp.code == "BANA" & sex2 == "M"))
i.BAWW.bad <- with(condition.Mhat[,],
     which(spp.code == "BAWW" & sex2 == "Unknown"))

i.all <- i.all[-which(i.all %in% c(i.BANA.bad,i.BAWW.bad)) ]


condition.Mhat2 <- condition.Mhat[i.all, ]

No captures of sp (site persist) individuals at Aceitillar for these spp

i.not.in.Aciet <- which(condition.Mhat2$spp.code %in% 
                        qc(BFGR,HLCU,
                           NOMO,RLTH,
                           YFGR))


condition.Mhat3 <- condition.Mhat2[-i.not.in.Aciet,]

Check output

with(condition.Mhat2, 
     table(site, spp.code))
with(condition.Mhat3, 
     table(site, spp.code))

Contrast modeling

Create contrast matrix

#for species that occur at all sites
contrst.mat.comm.5 <- rbind("Linear" = c(-2,-1, 0, 1, 2),
                          "Quad"   = c(+2,-1,-2,-1,+2))

#for species missing from one site
contrst.mat.comm.4 <- rbind("Linear1" = c(-3,-1, 1, 3),
                          "Quad"   = c(1,-1,-1,1))

Loop over species

Species for which there are data

m.hat.spp.to.use <- unique(condition.Mhat3$spp.code)

Fix site names

condition.Mhat3$site <- as.character(condition.Mhat3$site)
condition.Mhat3$site <- gsub("Cueva","La Cueva",condition.Mhat3$site)
condition.Mhat3$site <- gsub("Caoba","La Caoba",condition.Mhat3$site)
condition.Mhat3$site <- gsub("Corral","El Corral",condition.Mhat3$site)
summary(factor(condition.Mhat3$site))

Order factor labels

condition.Mhat3$site <- factor(condition.Mhat3$site,
                               levels = c("La Cueva",
                                           "La Caoba",
                                           "Morelia",  
                                           "El Corral",
                                           "Aceitillar"))

with(condition.Mhat3, 
     table(site, spp.code))

Loop

#store data in list
mhat.multcomp.list <- as.list(m.hat.spp.to.use)
names(mhat.multcomp.list) <- m.hat.spp.to.use
options(warn = -1)

##loop over each focal species
for(i in 1:length(m.hat.spp.to.use)){

  i.spp <- which(condition.Mhat3$spp.code == m.hat.spp.to.use[i])

  print(m.hat.spp.to.use[i])
  length(i.spp)

  #Note: use lm()
  mhat.anova.spp <- lm(M.hat.i ~ -1 + site ,
              data = condition.Mhat3[i.spp, ]
              )

  N.sites <- length(coef(mhat.anova.spp)) 
  coef(mhat.anova.spp)

 if(N.sites == 5) {
mhat.mult.comp.spp <- glht(mhat.anova.spp       
                          ,linfct =contrst.mat.comm.5
                          ,alternative = c("two.sided"))
} 

if(N.sites == 4){
  mhat.mult.comp.spp <- glht(mhat.anova.spp       
                          ,linfct =contrst.mat.comm.4
                          ,alternative = c("two.sided"))
}

 mhat.mult.comp.summary.spp <- summary(mhat.mult.comp.spp
                             ,test = adjusted(type = "none")) 

#print(mhat.mult.comp.summary.spp)

mhat.multcomp.list[[i]] <- mhat.mult.comp.summary.spp

}

Load function to extaction multocmp output

source('./R/getmultcompstats.R')

getmultcompstats(mhat.multcomp.list[[3]],dist = "")

extract

mhat.mcomp.stats.tab <- data.frame(Species = names(mhat.multcomp.list),
           Linear.B = NA,
           Quad.B = NA,
           Linear.t = NA,
           Quad.t = NA,
           Linear.p = NA,
           Quad.p = NA,

           Linear.SE = NA,
           Quad.SE = NA,

           siteLaCueva=NA,
           siteLaCaoba =NA,
           siteMorelia  =NA,
           siteEl.Corral=NA,
           siteAceitillar = NA,
                     #CI lo
           CI.lo.LaCueva=NA,
           CI.lo.LaCaoba =NA,
           CI.lo.Morelia  =NA,
           CI.lo.El.Corral=NA,
           CI.lo.Aceitillar = NA,

           #CI hi
           CI.hi.LaCueva=NA,
           CI.hi.LaCaoba =NA,
           CI.hi.Morelia  =NA,
           CI.hi.El.Corral=NA,
           CI.hi.Aceitillar = NA)

for(i in 1:length(mhat.multcomp.list)){
  print(i)
  mhat.mcomp.stats.tab[i,-1] <- getmultcompstats(mhat.multcomp.list[[i]], dist = "")
}

Set order

library(wrapr)
j.names <- qc(Species,
              Linear.B, Linear.SE, Linear.t, Linear.p,
              Quad.B,   Quad.SE,   Quad.t,   Quad.p,
              siteLaCueva,    CI.lo.LaCueva, CI.hi.LaCueva,
              siteLaCaoba,    CI.lo.LaCaoba, CI.hi.LaCaoba,
              siteMorelia,    CI.lo.Morelia, CI.hi.Morelia,
              siteEl.Corral,  CI.lo.El.Corral, CI.hi.El.Corral,
              siteAceitillar, CI.lo.Aceitillar, CI.hi.Aceitillar)

mhat.mcomp.stats.tab2 <- mhat.mcomp.stats.tab[,j.names]

Clean up

names(mhat.mcomp.stats.tab2) <- gsub("site","",names(mhat.mcomp.stats.tab2))
names(mhat.mcomp.stats.tab2) <- gsub("^La","La ",names(mhat.mcomp.stats.tab2))
names(mhat.mcomp.stats.tab2) <- gsub("^El.","El ",names(mhat.mcomp.stats.tab2))

write.csv(mhat.mcomp.stats.tab,file = "./tables/ind_spp_condition_trends_with_Aceitillar.csv")

Make plot

j.names <- qc(
   LaCueva, CI.lo.LaCueva, CI.hi.LaCueva,
   LaCaoba, CI.lo.LaCaoba, CI.hi.LaCaoba,
   Morelia, CI.lo.Morelia, CI.hi.Morelia,
   ElCorral, CI.lo.El.Corral, CI.hi.El.Corral,
   Aceitillar, CI.lo.Aceitillar, CI.hi.Aceitillar)

#the qc() function doesnt like spaces
j.names <- gsub("^La", "La ",j.names)
j.names <- gsub("^El", "El ",j.names)

mcomp.melt <- melt(data = mhat.mcomp.stats.tab2,
                   id.vars = "Species"#,
                   measure.vars = j.names)

#create site variable
mcomp.melt$Site <- mcomp.melt$variable
mcomp.melt$Site <- gsub("CI.lo.","",mcomp.melt$Site)
mcomp.melt$Site <- gsub("CI.hi.","",mcomp.melt$Site)
mcomp.melt$Site <- gsub("^La", "La ",mcomp.melt$Site)
mcomp.melt$Site <- gsub("^El", "El ",mcomp.melt$Site)
mcomp.melt$Site <- gsub("^El .Corral", "El Corral ",mcomp.melt$Site)

mcomp.melt$Site <- gsub("  ", " ",mcomp.melt$Site)


#create type of var column
mcomp.melt$variable2 <- mcomp.melt$variable
mcomp.melt$variable2 <- gsub("(CI.lo.)(.*)","\\1",mcomp.melt$variable2)
mcomp.melt$variable2 <- gsub("(CI.hi.)(.*)","\\1",mcomp.melt$variable2)
mcomp.melt$variable2[-grep("CI",mcomp.melt$variable2)] <- "mean"

Reshape

tab.forplotting <- dcast(data = mcomp.melt,
      formula = Species + Site ~ variable2,
      value.var = "value",
      fun.aggregate = mean)

mcomp.melt[which(mcomp.melt$Species == "PAWA"), ]
tab.forplotting[which(tab.forplotting$Species == "PAWA"), ]

Plot

tab.forplotting$Site <- gsub(" $","",tab.forplotting$Site)
tab.forplotting$Site2 <- gsub("La ","",tab.forplotting$Site)
tab.forplotting$Site2 <- gsub("El ","",tab.forplotting$Site2)
summary(factor(tab.forplotting$Site2))

tab.forplotting$Site2 <- factor(tab.forplotting$Site2,
        levels = c("Cueva",
                   "Caoba",
                   "Morelia",
                   "Corral",
                   "Aceitillar"))

ggplot(data = tab.forplotting,
       aes(y = mean,
           x = Site2,
           group = Species)) +
  geom_point() +
  geom_line() +
  xlab("Site") +
  ylab("Scaled Mass Index (SMI)") +
  geom_errorbar(aes(ymax = CI.hi.,
                    ymin = CI.lo.),
                width = 0) +
  facet_wrap(~Species,
             scale = "free")+ 
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

save.image("WRKSPC_condition_ANOVA.RData")

brouwern/DRmencia documentation built on May 6, 2019, 12:24 p.m.