#change your working directory to where you have the seedling survival file

#Generalized linear mixed effects models (GLMMs)
# EXAMPLE using SEEDLING SURVIVAL DATA

#read in data file
sdata=read.table("seedling.survival.data.txt",header=T,as.is=T,sep="\t")

#look at first few lines to see column names and some values
head(sdata)


#explore your data a bit more:

table(sdata$STATUS) #how many lived & died

range(sdata$CONS) #look at ranges and distributions of your predictors
hist(sdata$CONS,breaks=50)

range(sdata$LIGHT)
hist(sdata$LIGHT,breaks=10)

range(sdata$HEIGHT)
hist(sdata$HEIGHT,breaks=50)

#NOTE THAT YOUR PREDICTORS HAVE VERY DIFFERENT RANGES
# SO, IF YOU WANT TO DIRECTLY COMPARE THE EFFECTS OR
# TEST FOR INTERACTIONS, YOU SHOULD RESCALE THEM BY
# SUBTRACTING THE MEAND AND DIVIDING BY THE STANDARD DEVIATION:

sdata$LIGHT.ADJ=(sdata$LIGHT-mean(sdata$LIGHT))/sd(sdata$LIGHT)
sdata$HEIGHT.ADJ=(sdata$HEIGHT-mean(sdata$HEIGHT))/sd(sdata$HEIGHT)
sdata$CONS.ADJ=(sdata$CONS-mean(sdata$CONS))/sd(sdata$CONS)


# START ANALYZING:

#first load library to do mixed effects models
library(lme4)

# a simple model with only height as a predictor
surv.glm=glm(sdata$STATUS~sdata$HEIGHT.ADJ,family=binomial)
surv.glmer=lmer(sdata$STATUS~sdata$HEIGHT.ADJ+(1|sdata$PLOT),family=binomial)

summary(surv.glm)
summary(surv.glmer)

# Is it necessary to include PLOT as a random effect?
# Test using likelihood ratio test, as recommended by Bolker et al. (2009), Fig.1
# Look back at your notes to remember how to do a likelihood ratio test
# Hint, to extract the log likelihood value for your models in R use logLik()


# now compare adding species as a fixed vs. a random effect
surv.glm.spp=glm(sdata$STATUS~sdata$SPECIES+sdata$HEIGHT.ADJ,family=binomial)
surv.glmer.spp=lmer(sdata$STATUS~sdata$HEIGHT.ADJ+(1|sdata$PLOT)+(1|sdata$SPECIES),family=binomial)

summary(surv.glm.spp)
summary(surv.glmer.spp)

#now try playing around with adding other fixed effects (eg. light, conspecific neighbors)
#how do your results change when you add/remove PLOT as a random effect?

########################################################################

#Linear mixed effects models (=normal errors)
# EXAMPLE using FRUIT PRODUCTION DATA WITH REPEATED MEASURES

#lmer() function does not return p values because of uncertainties in 
# calculating the appropriate degrees of freedom

#read in data
fruitdata=read.table("fruit.data.txt",header=T,as.is=T,sep="\t")

#fruit production is log normally distributed, so we can make it normally
#distributed by log transforming it

hist(fruitdata$FRUIT)

fruitdata$LOGFRUIT=log(fruitdata$FRUIT+1)
hist(fruitdata$LOGFRUIT)

plot(fruitdata$DBH,fruitdata$LOGFRUIT)

#We have 3 years of fruit production data, with the same individuals measured in each year
# (ie. repeated measures) which we can account for by including individual as a random effect

fruit.lmer=lmer(fruitdata$LOGFRUIT~as.factor(fruitdata$YEAR)+fruitdata$DBH+(1|fruitdata$TAG))

summary(fruit.lmer)

#NOTE that for factors, the Estimates are relative to the Intercept, which is the value for the 
#first factor (in this case, YEAR 1)

#In this case, we do not test whether we should include individual (TAG) in the model
# because we know we need to control for repeated measures, so we should keep it in

#Since lmer does not give p values, we can evaluate significance using confidence intervals
# based on simulations, using two functions: mcmcsamp & HPDinterval

#mcmcsamp generates a sample of size n from the posterior distribution of the parameters  
# of our fitted model using Markov Chain Monte Carlo methods. 

fruit.mcmc=mcmcsamp(fruit.lmer,n=5000) 

#to then generate 95% "confidence" intervals, we use HPDinterval with prob=0.95  
# for each parameter, it returns the shortest interval with a 95% probability content in the empirical distribution. 
# HPDinterval is in two packages (code & lme4), and only the lme4 version seems to be working, so we have to specific lme4

lme4::HPDinterval(fruit.mcmc, prob=0.95)  

#Also, there is this function, in the languageR package, which does all of this for you in  one step:

library('languageR')
pvals.fnc(fruit.lmer,nsim = 10000)

# for more details that you really want on this issue, see: http://tolstoy.newcastle.edu.au/R/help/06/08/32409.html and related posts