## **************************************************
## 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=1)
hist(pred, col='red')

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

response <- rnorm(num, mean=effect.size.1*pred +
                         effect.size.2*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)

## now, lets create our response variable.
effect.size <- 2
effect.sd <- 0.25
response <- rnorm(num.total, mean=effect.size*pred, sd=effect.sd)

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)

## order rows within block (by pred)
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)
## plot ordered values within blocks
plot(dd$response, pch=16)
## --------------------------------------------------

## --------------------------------------------------
## 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

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




## **************************************************
## 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)

## now, lets create our response variable.
effect.size <- 0
effect.sd <- 0.25
response <- rnorm(num.total, mean=effect.size*pred, sd=effect.sd)

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
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)
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)
## --------------------------------------------------

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

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