knitr::opts_chunk$set(echo = TRUE)

Using lmer to model condition in all spp at same time

Use lmList() to fit seperate models to each species

Look at models fit to each subset of the data

lmList1 <- lmList(mass.log ~ wing.log|spp.code,
           data = condition[i.use,])

str(lmList1@.Data,1)


par(mfrow = c(4,4), mar = c(1,1,1,1))
for(i in 1:16){
  plot(lmList1@.Data[[i]], which = 2, main = "")
}

Model

i.use <- which(is.na(condition$stat.focals) == FALSE)
i.dont.use <- which(condition$spp.code %in% c("BFGR","HLCU","NOMO","RLTH","YFGR"))


condition$wing.log.cent <- scale(condition$wing.log, scale = F)
#condition$spp.site <- with()
lmer4blups <- blmer(mass.log~ 1 + 
                      (1|site.age) +
                (wing.log.cent|spp.code) +
                       (1|spp.code:site.age),
           data = condition[-i.dont.use,],
           control = lmerControl(optimizer="Nelder_Mead"))

summary(lmer4blups)

summary(factor(condition[-i.dont.use,"spp.code"]))
par(mfrow = c(1,1))
hist(resid(lmer4blups))

plot(resid(lmer4blups) ~lmer4blups@frame$spp.code)

plot(lmer4blups)

plot_model(lmer4blups, type = "diag")

Create new data frame for prediction

newdat <- lmer4blups@frame

#isolate unique values
i <- with(newdat, paste(spp.code, site.age))
i <- match(unique(i),i)
newdat[i,]
x <- predictInterval(lmer4blups,
                     newdata = newdat[i,],
                     which = "random",
                     level = 0.95,
                     #include.resid.var = FALSE,
                     type = "linear.prediction")

x <- cbind( newdat[i,],x)


x$site.age <- factor(x$site.age,
                                  levels = c("2","5","10","20",
                                             "mature"))
i.bana <- which(x$spp.code == "BANA")
i.gabu <- which(x$spp.code == "GABU")
i.bcpt <- which(x$spp.code == "BCPT")
i.grwa <- which(x$spp.code == "GRWA")

ggplot(data = x[i.bana,],
        aes(y = fit,
            x = site.age,
            color = site.age)) +
  geom_point() +
  #facet_wrap(~spp.code, scales = "free") +
  #geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lwr,
                    ymax = upr),
                width = 0) +
  coord_flip()
library(nlme)
#constant variance(s), generally used to allow different variances according to the levels of a classification factor

gls1 <- lme(log(mass) ~ log(wing), 
            random = list(~1 +log(wing)|spp.code),
           data = condition[i.use,],
           weights = varFixed(~log(wing)))
plot_model(gls1, type = "diag")

Run models

lmer.site <- update(lmer1, . ~ . + site)
lmer.site.status <- update(lmer1, . ~ . + site + stat.focals)

lmer.site.X.status <- update(lmer1, . ~ . + site*stat.focals)
lmer.means <- update(lmer1, . ~ . + -1 + site:stat.focals)
AICtab(lmer1,
       lmer.site,
       lmer.site.status,
       lmer.site.X.status,
       lmer.means)
plot_model(lmer.site.X.status,
           type = "est")

FEsim.out <- FEsim(lmer.means)

ggplot(data = FEsim.out[FEsim.out$term != "wing.log",],
        aes(y = mean,
            x = term)) +
  geom_point()
newdat <- lmer1@frame

i <- with(newdat, paste(spp.code, site))
i <- match(unique(i),i)
newdat[i,]

x <- predictInterval(lmer1,newdata = newdat[i,],
                     which = "random")

x <- cbind( newdat[i,],x)
ggplot(data = x,
        aes(y = fit,
            x = site,
            color = site)) +
  geom_point() +
  facet_wrap(~spp.code, scales = "free") +
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lwr,
                    ymax = upr),
                width = 0) +
  coord_flip()

predictInterval

lmer4blups <- lmer(mass.log~ 1 + 
                (1|spp.code) +    #slope for each spp
                (1 + wing.log.cent|spp.code:site),#intercept for each spp
           data = condition[i.use,])

summary(lmer4blups)


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