rm(list=ls())

library(lme4)
library(visreg)

## **************************************************
## Repeatability of a sexual signal trait
## **************************************************

## --------------------------------------------------
## Read and examine the data

## 1.
ff <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/flycatcher.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## 2.
head(ff)

## 3.
interaction.plot(response=ff$patch,
                 x.factor=ff$year,
                 trace.factor=ff$bird,
                 legend=FALSE,
                 lty=1,
                 xlab='Year', ylab='Patch size', 
                 type='b', pch=16, las=1)

## --------------------------------------------------
## Fit a linear mixed-effects model

## 1.
out <- lmer(patch ~ 1 + (1|bird), data=ff)
## also create a summary object (useful for quick extraction of
## parameter estimates)
summ <- summary(out)
summ

## 4.
##
## repeatability is a measure of how variable birds are relative to
## our measurements.  If measurements are highly variable and birds
## are not that different, we would have low repeatability (e.g., most
## of the variance in our samples would simply be due measurement
## error).  If measurements error is small and birds differ a lot, we
## would have high repeatability (e.g., most of the variance in our
## samples would be due to actual differences in birds).  Thus, to
## calculate repeatability, we can calculate the ratio:
##
## variance due to actual differences / (variance due to actual
## differences + variance due to measurement error)
##
## for our model output, this is given by
1.243 / (1.243 + 0.358)
## =0.776

## 6.
plot(x=fitted(out), y=resid(out), pch=16,
     xlab='Fitted value', ylab='Residual', las=1)

## **************************************************
## Goldie's vision
## **************************************************

## --------------------------------------------------
## Read and examine the data

## 1.
gg <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/goldfish.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## 2.
interaction.plot(response=gg$sensitivity,
                 x.factor=gg$wavelength,
                 trace.factor=gg$fish,
                 legend=FALSE,
                 lty=1,
                 xlab='Wavelength', ylab='Sensitivity', 
                 type='b', pch=16, las=1)

## or make the plot manually
plot(NA, xlim=c(1,9), ylim=range(gg$sensitivity),
     xlab='Wavelength', ylab='Sensitivity', las=1)
lines(gg$sensitivity[gg$fish=='fish1'], type='l')
lines(gg$sensitivity[gg$fish=='fish2'], col='red')
lines(gg$sensitivity[gg$fish=='fish3'], col='green4')
lines(gg$sensitivity[gg$fish=='fish4'], col='blue')
lines(gg$sensitivity[gg$fish=='fish5'], col='purple')

## 3.
##
## It is a "subjects-by-treatment" repeated measures design, since
## each fish is measured once under each treatment. It is essentially
## the same as a randomized complete block design (think of the
## individual fish as "blocks").

## --------------------------------------------------
## Fit a linear mixed-effects model

## 1.
out <- lmer(sensitivity ~ wavelength + (1|fish), data=gg)
summ <- summary(out)
summ

gg <- gg[which(gg$wavelength=='nm585'),]
out <- lmer(sensitivity ~ wavelength + (1|fish), data=gg)
summ <- summary(out)

## 2.
visreg(out, xvar='wavelength')

## 3.
##
## raw data
gg$sensitivity
## fitted values + residuals
fitted(out)+resid(out)
## plot, just to check
plot(x=gg$sensitivity, y=fitted(out)+resid(out),
     xlab='data', ylab='fitted value + residual', las=1, pch=16)
abline(a=0,b=1)

## 4.
plot(x=fitted(out), y=resid(out),
     xlab='fitted value', ylab='residual', las=1, pch=16)
abline(a=0, b=0, lty=2, col='red')

## 5.
fixef(out)

## 6.
VarCorr(out) ## gives us our variances (these are also presented as
             ## part of the model summary)

## note, we can calculate the residual standard deviation manually,
## using our residuals - we divide by the number of degrees of freedom
## which is 45 (num obs) - 8 (num params for estimated fixed effect -
## 1 (num params for estimated random effect).
sqrt(sum(resid(out)^2)/(45-8-1))

## typically the quantity
sd(ranef(out)$fish[,'(Intercept)'])
## would give the standard deviation of the random effect, but here
## our variance is so small that rounding error, etc is too large

## **************************************************
## Yukon yarrow
## **************************************************

## --------------------------------------------------
## Visualize the data

## 1.
kk <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/kluane.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## 2.
head(kk)

## 3.
## The experiment used a split-plot design, in which whole plots were
## randomly assigned different treatments, and then different levels
## of a second treatment (duration) were assigned to plot halves.

## 4.
## here's a ggplot solution to mix things up
library(tidyverse)
ggplot(data=kk, aes(x=duration, y=phen.ach, colour=duration)) + 
  facet_grid(. ~ treatment, switch="x") +
  geom_line(aes(group=plot), colour="black") +
  geom_point(size=3) +
  labs(x="treatment", y="phen.ach") +
  theme_classic() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

## --------------------------------------------------
## Fit a linear mixed-effects model

## 1.
out <- lmer(log(phen.ach) ~ treatment + duration + (1|plot), data=kk)
summ <- summary(out)
summ

## check a couple model diagnostics
##
## plot residuals
plot(x=fitted(out), y=resid(out),
     xlab='Fitted Values', ylab='Residuals', las=1, pch=16)
abline(h=0, col='red')
## not a great fit - residuals appear to be positive for intermediate
## values, and lower for extreme values

## qq-plot
qqnorm(as.vector(resid(out)), pch=16, main='', xlab='', ylab='')
qqline(as.vector(resid(out)))
## this one looks ok

## plot fitted values with visreg
visreg(out, xvar='duration', by='treatment')

## 2.
out <- lmer(log(phen.ach) ~ treatment * duration + (1|plot), data=kk)
summ <- summary(out)
summ
visreg(out, xvar='duration', by='treatment')

## In contrast to (1) above, which had no interaction, here the effect
## size of duration can differ between treatments
##
## Based on a visual inspection, this model with interaction fits the
## data much better, as the fitted values appear to take on more
## appropriate values.

## 3.
plot(x=fitted(out), y=resid(out),
     xlab='Fitted Values', ylab='Residuals', las=1, pch=16)
abline(h=0, col='red')
## this definitely looks better than before

## qq-plot
qqnorm(as.vector(resid(out)), pch=16, main='', xlab='', ylab='')
qqline(as.vector(resid(out)))
## alos looks fine

## 4.
fixef(out)
ranef(out)
VarCorr(out)

## 5.
drop1(out, test='Chisq')

## could use 'nlme' package which will return p-values
library(nlme)
out <- lme(fixed=log(phen.ach) ~ treatment + duration,
           random= ~1|plot,
           data=kk)
summary(out)