lme4_GH858.R

library(lme4)
load("areaRep26_glmer.Rdata")
n <- nrow(areaRep26_glmer)
nsp <- length(unique(areaRep26_glmer$especie))
## experimental design: each species appears in exactly 1 'sistema'
##  each species appears once in each 'tempo'
with(areaRep26_glmer,table(especie,sistema))
with(areaRep26_glmer,table(especie,tempo))
summary(areaRep26_glmer)

library(ggplot2)

## exploratory graphs
## rearrange data for graphing purposes
dd <- areaRep26_glmer
dd$especie <- reorder(dd$especie,dd$area)
(ggplot(dd, aes(especie,area))
    + geom_boxplot()
    + coord_flip()
    + scale_y_sqrt()
)
(ggplot(dd, aes(tempo,area))
    + geom_boxplot()
    + coord_flip()
    + geom_hline(yintercept=0,lty=2)
    + facet_wrap(~sistema)
    + scale_y_sqrt()
)

## fit linear model to get means by species
m1 <- lm(area ~ especie-1, data=areaRep26_glmer)
s1 <- log(coef(m1))[factor(areaRep26_glmer$especie)]

## try a GLM (this failed with various choices of mustart/start,
## default args: "cannot find valid starting values"
## start=rep(1,nsp):
##    Error: inner loop 1; cannot correct step size
##  but etastart=s1 makes it work
m2 <- glm(area ~ especie-1, data=areaRep26_glmer,
          family = gaussian(link= "log"),
          etastart=s1)

m0 <- glmer(area ~ (1|especie), data=areaRep26_glmer,
            family = gaussian(link= "log"),
            etastart=s1)

update(m0, . ~ tempo + (1|especie))
## Warning message:
## In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge with max|grad| = 0.0403883 (tol = 0.002, component 

## works if we rescale areas to make them less extreme:
m_tempo <- update(m0, area/1e5 ~ tempo + (1|especie))
m_temposist <- update(m_tempo, area/1e5 ~ . + sistema)
aa <- allFit(m_temposist)
summary(aa)
## looks like Nelder_mead and nmkbw are slightly better (0.01 log-likelihood units)

m_temposist_NM <- update(m_temposist,
                         control=glmerControl(optimizer="Nelder_Mead"))
## not much practical difference

## warnings 
m_sist <- update(m_tempo, . ~ . + sistema)
##  Model failed to converge with max|grad| = 0.0937365 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
m_int <- update(m_tempo, . ~ tempo * sistema + (1|especie))