rm(list=ls())
library(lme4)

## 
## EXAMPLE 1: fitting and plotting a quadratic model
## 

set.seed(1)

## --------------------------------------------------
## simulate data - we're simulating with very strong trends here so
## that everything is easy to see!
##
## Lets create a simple numeric predictor
num <- 100
pred <- rnorm(100, mean=0, sd=2)
hist(pred, col='red')

## Now, lets create our response variable.  Let's assume a strong
## positive effect of pred.
alpha          <- 2
beta.linear    <- 2
beta.quadratic <- 1
## sd around this response
effect.sd <- 0.5

response <- rnorm(num,
                  mean=alpha + beta.linear*pred + beta.quadratic*pred^2,
                  sd=effect.sd)

## package into a data.frame
dd <- data.frame(pred=pred,
                 pred2=pred^2,
                 response=response)
head(dd)
## --------------------------------------------------

## --------------------------------------------------
## plot the raw data
plot(x=dd$pred, y=dd$response, xlab='pred', ylab='response',
     col='red', pch=16, las=1)
## --------------------------------------------------

## --------------------------------------------------
## fit a model
out <- lm(response ~ pred + pred2, data=dd)
summary(out)

## extract model coefficients
coeffs <- coef(out)
b0 <- coeffs['(Intercept)']
b1 <- coeffs['pred']
b2 <- coeffs['pred2']

## add the best fit line
curve(b0 + b1*x + b2*x^2,
      from=min(dd$pred), to=max(dd$pred), col='blue', add=TRUE)
## --------------------------------------------------



## 
## EXAMPLE 2: random effect (neglecting it obscures trend)
## 

set.seed(1)

## --------------------------------------------------
## simulate data
##
num.blocks <- 20
num.per.block <- 5
num.total <- num.blocks*num.per.block

block <- rep(1:num.blocks, each=num.per.block)

## create predictor
pred <- rnorm(num.total, mean=0, sd=1)

## simulate response variable (no block effect yet)
effect.size <- 2
effect.sd <- 0.25
response <- rnorm(num.total, mean=effect.size*pred, sd=effect.sd)
## simulate a block effect
block.effect <- rnorm(num.blocks, mean=0, sd=10)
block.effect.by.ind <- rep(block.effect, each=num.per.block)
## add block effect to the response
response <- response+block.effect.by.ind

## package all this up into a data.frame
dd <- data.frame(block=block, pred=pred, response=response)
head(dd, 20)

## order rows within block (by pred) to facilitate plotting
dd <- dd[order(dd$pred),]
dd <- dd[order(dd$block),]
## --------------------------------------------------

## --------------------------------------------------
## plot the raw data
##
## plot with 'pred' as x-axis
plot(x=dd$pred, y=dd$response, pch=16, col=dd$block)
## plot ordered values within blocks
plot(dd$response, pch=16, col=dd$block)
## --------------------------------------------------

## --------------------------------------------------
## fit a model without any random effects
out <- lm(response ~ pred, data=dd)
summary(out)
drop1(out, test='F') ## no significant effect of pred detected and effect size is incorrect

## fit a model with random effect
out <- lmer(response ~ pred + (1|block), data=dd)
summary(out)
drop1(out, test='Chisq') ## significant (and accurate) effect of pred correctly detected!
## --------------------------------------------------

## --------------------------------------------------
## could also now plot best fit lines for each block (using out from
## the second model above)
intercept.block <- coef(out)$block[['(Intercept)']]
intercept.total <- fixef(out)['(Intercept)']
slope <- fixef(out)['pred']
plot(x=dd$pred, y=dd$response, pch=16, col=dd$block)
## add farm-specific lines
for(ii in 1:length(intercept.block))
  curve(intercept.block[ii] + slope * x,
        from=min(dd$pred), to=max(dd$pred), lwd=0.3, add=TRUE, col=ii)
## add mean line across all farms
curve(intercept.total + slope * x,
      from=min(dd$pred), to=max(dd$pred), lwd=2, add=TRUE)
## --------------------------------------------------


## 
## EXAMPLE 3: random effect (neglecting it creates trend)
## 

set.seed(1)

## --------------------------------------------------
## simulate data
##
num.blocks <- 3
num.per.block <- 20
num.total <- num.blocks*num.per.block

block <- rep(1:num.blocks, each=num.per.block)

## create predictor
pred <- rnorm(num.total, mean=0, sd=1)

## simulate response variable (no block effect yet)
effect.size <- 0
effect.sd <- 0.25
response <- rnorm(num.total, mean=effect.size*pred, sd=effect.sd)
## simulate a block effect
block.effect <- rnorm(num.blocks, mean=0, sd=10)
block.effect.by.ind <- rep(block.effect, each=num.per.block)
## add block effect to both the predictor and the response (so farms
## differ in their predictor values and their response values)
pred <- pred+block.effect.by.ind
response <- response+block.effect.by.ind

## package all this up into a data.frame
dd <- data.frame(block=block, pred=pred, response=response)
head(dd)

## order rows within block (by pred) to facilitate plotting
dd <- dd[order(dd$pred),]
dd <- dd[order(dd$block),]
## --------------------------------------------------

## --------------------------------------------------
## plot the raw data
##
## plot ordered values within blocks
plot(dd$response, pch=16)
## plot with 'pred' as x-axis
plot(x=dd$pred, y=dd$response, pch=16, col=dd$block)
## --------------------------------------------------

## --------------------------------------------------
## fit a model without any random effects
out <- lm(response ~ pred, data=dd)
summary(out)
drop1(out, test='F') ## significant effect of pred (incorrect) and effect size is incorrect

## fit a model with random effect
out <- lmer(response ~ pred + (1|block), data=dd)
summary(out)
drop1(out, test='Chisq') ## significant effect of pred (correct)
## --------------------------------------------------


## 
## EXAMPLE 4: exploring sphericity
## 

set.seed(1)

## --------------------------------------------------
## simulate data
##
num.individuals <- 5
num.treatments <- 3
num.total <- num.individuals*num.treatments

## create individual
individual <- rep(1:num.individuals, num.treatments)
## create treatment
treatment <- rep(1:num.treatments, each=num.individuals)

## simulate response variable (no block effect yet)
effect.size <- c(0,1,-3)
effect.sd <- 0.25
response <- rnorm(num.total, mean=effect.size[treatment], sd=effect.sd)
## simulate an individual effect
individual.effect <- rnorm(num.total, mean=0, sd=3)
## add individual effect to the response
response <- response + individual.effect

## package all this up into a data.frame
dd <- data.frame(individual=individual,
                 treatment=treatment,
                 response=response)
dd

## --------------------------------------------------
## plot the raw data
##
## plot with 'pred' as x-axis
plot(x=dd$treatment, y=dd$response,
     pch=as.character(dd$individual),
     col=dd$individual,
     xlab='Treatment', ylab='Response', las=1)
## --------------------------------------------------

## quick way to figure out all unique combinations
mm <- matrix(NA, nrow=num.treatments, ncol=num.treatments)
combns <- which(upper.tri(mm), arr.ind=TRUE)
diffs <- apply(combns, 1, function(x)
  response[treatment==x[1]]-response[treatment==x[2]])
diffs
apply(diffs, 2, var)
## sphericity is the assumption that the above variances are all
## approximately equal
## --------------------------------------------------