---
title: "Lecture 9: Interactions"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## Before we start...



# Categorical Variables, aka Factors
```{r GetCake, results='asis', echo=FALSE, message=FALSE}
library(lme4)
data("cake")
Sugar <- cake[cake$recipe!="B" & cake$temp=="225",] # remove recipe B & use temperature of 175C
Sugar$treatment <- c("Control", "Treatment")[1+grepl("A", Sugar$recipe)]
```

Start by looking again at 2 variables

# Two categorical Variables

We do not have to only consider one variable: just like multiple regression, we can consider several.

e.g. with the cakes, temperature was also a factor that was controlled

```{r CakesTwoWayData, echo=TRUE, comment=NA}
 # use temperatures of 175C & 225C
Cakes2 <- cake[cake$temp%in%c("175", "225") ,]
```


# Two categorical Variables

```{r PlotCakesTwoWayData, echo=TRUE, comment=NA, fig.height=5}
par(mar=c(4.1,6,1,1), cex=1.5)
boxplot(Cakes2$angle~Cakes2$recipe + Cakes2$temp, 
        horizontal=TRUE, las=1)
```

# Two categorical Variables

The contrasts for the temperature can be written in the same way

```{r TwoWayDataTempContrast, echo=TRUE, comment=NA}
Cakes2$tempF <- factor(Cakes2$temp)
contrasts(Cakes2$tempF)
```

# Two categorical Variables

The design matrix is now like this:

```{r DesignMatrixTwoWayData, echo=TRUE, comment=NA}
head(model.matrix(~recipe + tempF, data=Cakes2))

```

We now just have an extra column for temperature


# Another way of writing the model

In matrix form (and back soon to something simpler)

$$
\mathbf{Y} = X \mathbf{\beta} + \mathbf{\varepsilon}
$$



$$
\left( \begin{array}{c}
   y_1 \\
   y_2 \\
   y_3 \\
   y_4 \\
   y_5 \\
   y_6 \\
\vdots \\
  \end{array}  \right) = \left( \begin{array}{cccc}
   1 & 0 & 0 & 0 \\
   1 & 0 & 0 & 1 \\
   1 & 1 & 0 & 0 \\
   1 & 1 & 0 & 1 \\
   1 & 0 & 1 & 0 \\
   1 & 0 & 1 & 1 \\
\vdots & \vdots & \vdots \\
  \end{array}  \right)
  \left( \begin{array}{c}
   \beta_0 \\
   \beta_{1,2} \\
   \beta_{1,3} \\
   \beta_{2,2} \\
\vdots \\
  \end{array}  \right)
 +   \left( \begin{array}{c}
   \varepsilon_1 \\
   \varepsilon_2 \\
   \varepsilon_3 \\
   \varepsilon_4 \\
   \varepsilon_5 \\
   \varepsilon_6 \\
\vdots \\
  \end{array}  \right)
$$

# Parameters

Where 

- $\beta_0$ is the intercept
- $\beta_{1,2}$ and $\beta_{1,3}$ are the contrasts between the control and treatments 2 & 3 (hot chocolate & suger)
- $\beta_{2,2}$ is the cintrast between 225&deg; and 175&deg;

We could write this with $\beta_{1,1} = \beta_{2,1} = 0$ (the contrast of a treatment with itself)

#  Multiplying the matrices

Doing the matrix multiplication we get 

$$
y_1 = 1 \cdot \beta_0 + 0 \cdot \beta_{1,2} + 0 \cdot \beta_{1,3} + 0 \cdot \beta_{2,2} + \varepsilon_1 = \beta_0  + \varepsilon_1
$$
for the first row, and for the rest...

$$
\begin{aligned}
   y_2 &= 1 \cdot \beta_0 + 0 \cdot \beta_{1,2} + 0 \cdot \beta_{1,3} + 1 \cdot \beta_{2,2} + \varepsilon_2 
   &= \beta_0 + \beta_{2,2} + \varepsilon_2\\
   y_3 &= 1 \cdot \beta_0 + 1 \cdot \beta_{1,2} + 0 \cdot \beta_{1,3} + 0 \cdot \beta_{2,2} + \varepsilon_3
   &= \beta_0 + \beta_{1,2} + \varepsilon_3\\
   y_4 &= 1 \cdot \beta_0 + 1 \cdot \beta_{1,2} + 0 \cdot \beta_{1,3} + 1 \cdot \beta_{2,2} + \varepsilon_4
   &= \beta_0 + \beta_{1,2} + \beta_{2,2} + \varepsilon_4 \\
   y_5 &= 1 \cdot \beta_0 + 0 \cdot \beta_{1,2} + 1 \cdot \beta_{1,3} + 0 \cdot \beta_{2,2} + \varepsilon_5
   &= \beta_0 + \beta_{1,3} + \varepsilon_5\\
   y_6 &= 1 \cdot \beta_0 + 0 \cdot \beta_{1,2} + 1 \cdot \beta_{1,3} + 1 \cdot \beta_{2,2} + \varepsilon_6
   &= \beta_0 + \beta_{1,3} + \beta_{2,2} + \varepsilon_6
\end{aligned}
$$

- every data point has $\beta_0$ and 0 or 1 $\beta_{1,\cdot}$  and 0 or 1 $\beta_{2,\cdot}$


# Two categorical Variables: Fitting the model

We can fit the model, just like last time!

```{r FitTwoWayModel, echo=TRUE, comment=NA}
mod.2way <- lm(angle~recipe + tempF, data=Cakes2)
print(summary(mod.2way)$coefficients, digits=2)

```

Little effect of recipe, big effect of temperature

- higher temperature means cake bends more before breaking

# Why the intercept?

What happens if we remove the intercept?

```{r FitTwoWayModelNoInt, echo=TRUE, comment=NA}
mod.2way <- lm(angle~recipe + tempF-1, data=Cakes2)
print(summary(mod.2way)$coefficients, digits=2)

```

We get 3 levels of recipe, but still only one of tempF

# Why the intercept?

Without the intercept, recipe has 3 parameters that need to be estimated: $\beta^{*}_{1,1}$, $\beta^{*}_{1,2}$, $\beta^{*}_{1,3}$

- these are the means for the 3
- with the intercept we have $\beta_{0}$, $\beta_{1,2}$, and $\beta_{1,3}$

For the one-way model, this is fine. With the two-way model and no intercept, we would have 

- $\beta^{*}_{1,1}$, $\beta^{*}_{1,2}$, $\beta^{*}_{1,3}$ for recipe and 
- $\beta^{*}_{2,1}$, $\beta^{*}_{2,2}$, for temperature.

# Identifiability

The means for each combination are

$$
\begin{aligned}
\mu_{1,1} &= \beta^{*}_{1,1} + \beta^{*}_{2,1} \\
\mu_{1,2} &= \beta^{*}_{1,1} + \beta^{*}_{2,2}\\
\mu_{2,1} &= \beta^{*}_{1,2} + \beta^{*}_{2,1}\\
\mu_{2,2} &= \beta^{*}_{1,2} + \beta^{*}_{2,2}\\
\mu_{3,1} &= \beta^{*}_{1,3} + \beta^{*}_{2,1}\\
\mu_{2,2} &= \beta^{*}_{1,3} + \beta^{*}_{2,2}
\end{aligned}
$$

The problem now is that I can add $C$ to every $\beta^{*}_{1,\cdot}$ and subtract $C$ from every $\beta^{*}_{2,\cdot}$ and get the same $\mu$

Thus we need at least one constraint, hence the second factor is set up as a contrast

# Categorical and Continuous Data

The models for categorical data have been built by using the continuous data model, so thre is no reason why we cannot combine them.

For example


$$
\left( \begin{array}{c}
   y_1 \\
   y_2 \\
   y_3 \\
   y_4 \\
   y_5 \\
   y_6 \\
\vdots \\
  \end{array}  \right) = \left( \begin{array}{cccc}
   1 & 0 & 0 & 2.4 \\
   1 & 0 & 0 & 8.3 \\
   1 & 1 & 0 & 9.7 \\
   1 & 1 & 0 & 4.6 \\
   1 & 0 & 1 & 2.3 \\
   1 & 0 & 1 & 1.0 \\
\vdots & \vdots & \vdots \\
  \end{array}  \right)
  \left( \begin{array}{c}
   \beta_0 \\
   \beta_{1,2} \\
   \beta_{1,3} \\
   \beta_{2} \\
\vdots \\
  \end{array}  \right)
 +   \left( \begin{array}{c}
   \varepsilon_1 \\
   \varepsilon_2 \\
   \varepsilon_3 \\
   \varepsilon_4 \\
   \varepsilon_5 \\
   \varepsilon_6 \\
\vdots \\
  \end{array}  \right)
$$

So $\beta_2$ is the slope of the effect of the continuous covariate.

# Simulate Some Data

```{r ANCOVAData, echo=TRUE}
Cols=c(1,2,4)
x1 <- factor(rep(1:3, each=50))
x2 <- rnorm(length(x1), 1.5,1)
mu <- 1 - 2*(x1==2) -4*(x1==3) + (1.5)*x2
y <- rnorm(length(mu), mu, 1)
```

# Plot Some Data


```{r PlotANCOVAData, fig.width=5, fig.height=5, echo=FALSE}
plot(x2, y, col=Cols[x1], pch=20)
legend(0.9*min(x2),0.98*max(y), c("x1=0", "x1=1", "x1=2"), col=Cols, pch=20)
sapply(1:3, function(i) abline(lm(y[x1==i]~x2[x1==i]), col=Cols[i]))

```

# Fitting the model

```{r FitANCOVA, echo=TRUE, fig.width=5, fig.height=3.6}
mod.anc <- lm(y ~x1 + x2)
round(summary(mod.anc)$coefficients,3)
```

- x12 is the contrast between levels 1 & 2 of x1
- x2 is the slope of the relationship with x2

# Interactions

So far we have assumed that all effects act independently

e.g.

$$
y_i = \alpha + \beta_{1} x_{i1} + \beta_{2} x_{i2} + \varepsilon_i
$$

But this makes some strict assumptions about the model


# No Interactions

```{r NoInt, fig.width=5, fig.height=4}
alpha <- 0; beta <- c(0, 1.5, 0, -1, 2, 2); sigma <- 1
x <- expand.grid(x1=1:5, x2=0:1, rep=1:30)
mu <- alpha + beta[x$x1] + x$x2*beta[6]
y <- rnorm(nrow(x), mu, sigma)

Mean <- aggregate(y, list(x1=x$x1, x2=x$x2), mean)
Mean$at.x <- Mean$x1+0.1*Mean$x2 - 0.05
SD <- aggregate(y, list(x1=x$x1, x2=x$x2), sd)

plot(jitter(x$x1, amount=0.05), y, col=c("#ffcc33", "#9933ff")[1+x$x2], pch=1, cex=0.5, xlab="x1", las=1)
points(Mean$at.x, Mean$x, col=c("#b38f24", "#6b24b3")[1+x$x2], pch=16, cex=1)
segments(Mean$at.x, Mean$x+SD$x, Mean$at.x, Mean$x-SD$x, col=c("#b38f24", "#6b24b3")[1+x$x2], lwd=1.5)
lines(Mean$at.x[Mean$x2==0], Mean$x[Mean$x2==0], col="#b38f24", lwd=1.5)
lines(Mean$at.x[Mean$x2==1], Mean$x[Mean$x2==1], col="#6b24b3", lwd=1.5)

```

# But sometimes....

```{r Int, fig.width=5, fig.height=4}
# alpha <- 0; beta <- c(1,1.5); sigma <- 1
# x <- expand.grid(x1=0:1, x2=0:1, rep=1:30)
# mu <- alpha + as.matrix(x[,c("x1","x2")])%*%beta + gamma*apply(x[,c("x1","x2")], 1, prod)
# y <- rnorm(nrow(x), mu, sigma)

alpha <- 0; beta <- c(0, 1.5, 0, -1, 2, 2); gamma <- -1; sigma <- 1
x <- expand.grid(x1=factor(1:5), x2=factor(0:1), rep=1:30)
mu <- alpha + beta[x$x1] + (x$x2==1)*beta[6] + gamma*beta[x$x1]*(x$x2==1)
y <- rnorm(nrow(x), mu, sigma)


Mean <- aggregate(y, list(x1=x$x1, x2=x$x2), mean)
Mean$at.x <- as.numeric(Mean$x1)+0.1*(Mean$x2==1) - 0.05
SD <- aggregate(y, list(x1=x$x1, x2=x$x2), sd)

plot(jitter(as.numeric(x$x1), amount=0.05), y, col=c("#ffcc33", "#9933ff")[1+(x$x2==1)], 
     pch=1, cex=0.5, xlab="x1", las=1)
points(Mean$at.x, Mean$x, col=c("#b38f24", "#6b24b3")[1+(x$x2==1)], pch=16, cex=1)
segments(Mean$at.x, Mean$x+SD$x, Mean$at.x, Mean$x-SD$x, col=c("#b38f24", "#6b24b3")[1+(x$x2==1)], lwd=1.5)
lines(Mean$at.x[Mean$x2==0], Mean$x[Mean$x2==0], col="#b38f24", lwd=1.5)
lines(Mean$at.x[Mean$x2==1], Mean$x[Mean$x2==1], col="#6b24b3", lwd=1.5)

```

The effect of X1 depends on X2


# Interactions

These are called *interactions*

- this is a "first order interaction"
- can get more complicated


# Interactions in R

These are (too) easy to code. If we want the interaction between X1 and X2, we use

```{r TwoDay, eval=FALSE, echo=TRUE}
lm(Y ~ X1*X2)
```

we can also have interactions between 3 variables:

```{r ThreeDay, eval=FALSE, echo=TRUE}
lm(Y ~ X1*X2*X3)
```

...or even more


# Estimates

```{r TwoDayFit, echo=TRUE}
mod.2way <- lm(y ~ x$x1*x$x2)
summary(mod.2way)$coefficients

```

# Interpretation 

```{r TwoDayFitInterpret, echo=TRUE}
round(summary(mod.2way)$coefficients[c(1:3, 6:8),], 2)

```

- x1=0, x2=0: (Intercept)
- x1=1, x2=0: (Intercept) + x1
- x1=0, x2=1: (Intercept) \..... + x2
- x1=1, x2=1: (Intercept) + x1 + x2 + x1:x2

# Why the colon?

When we write A*B, R interprets it as the effect of A + the effect of B + the interaction

i.e. as A + B + A:B

So A:B is (just) the interaction

If we fit a model with just A:B, we assume that the average effect of A (and B) is exactly 0. This is almost never true (I have seen it once, but I had to prove it)

# Plot

```{r TwoDayPlot1, echo=TRUE, fig.height=3}
plot(jitter(as.numeric(x$x1), amount=0.05), y, 
     col=c("#ffcc33", "#9933ff")[1+(x$x2==1)], pch=1, cex=0.5, xlab="x1", las=1)
points(Mean$at.x, Mean$x, col=c("#b38f24", "#6b24b3")[1+(x$x2==1)], pch=16, cex=1)
segments(Mean$at.x, Mean$x+SD$x, Mean$at.x, Mean$x-SD$x, col=c("#b38f24", "#6b24b3")[1+(x$x2==1)], lwd=1.5)
lines(Mean$at.x[Mean$x2==0], Mean$x[Mean$x2==0], col="#b38f24", lwd=1.5)
lines(Mean$at.x[Mean$x2==1], Mean$x[Mean$x2==1], col="#6b24b3", lwd=1.5)

```

# Interactions with continuous variables

```{r ANCOVADataInt, echo=TRUE}
Cols=c(1,2,4)
x1 <- factor(rep(1:3, each=50))
x2 <- rnorm(length(x1), 1.5,1)
mu <- 1 - 2*(x1==2) -2*(x1==3) + (1.5 - 1.3*(x1==3))*x2
y <- rnorm(length(mu), mu, 1)
```

# Plot Some Data


```{r PlotANCOVADataInt, fig.width=5, fig.height=5, echo=FALSE, message=FALSE, results='hide'}
plot(x2, y, col=Cols[x1], pch=20)
legend(0.9*min(x2),0.98*max(y), c("x1=0", "x1=1", "x1=2"), col=Cols, pch=20)
sapply(1:3, function(i) abline(lm(y[x1==i]~x2[x1==i]), col=Cols[i]))

```

# Fitting the model

```{r FitANCOVAInt, echo=TRUE, fig.width=5, fig.height=3.6}
mod.anc <- lm(y ~x1*x2)
round(summary(mod.anc)$coefficients,3)
```

# A problem of Interpretation

```{r FitANCOVAInt2, echo=TRUE, fig.width=5, fig.height=3.6}
mod.anc <- lm(y ~x1*x2)
mod.ancS <- lm(y ~x1*scale(x2, scale=FALSE))
round(cbind(Orig=coef(mod.anc), Centred=coef(mod.ancS)),3)
```

# Plot Again


```{r PlotANCOVADataIntAgain, fig.width=5, fig.height=5, echo=FALSE, message=FALSE, results='hide'}
plot(x2, y, col=Cols[x1], pch=20)
legend(0.9*min(x2),0.98*max(y), c("x1=0", "x1=1", "x1=2"), col=Cols, pch=20)
sapply(1:3, function(i) abline(lm(y[x1==i]~x2[x1==i]), col=Cols[i]))

abline(v=0, lty=3, lwd=2, col="grey40")
abline(v=mean(x2), lty=3, lwd=2, col="grey40")

```



# Nested Effects

Another Complication!

An interaction means that any observation can take any value of all levels 

- all recipes can be cooked at all temperatures

But some data is nested

e.g. regions within countries

- Trøndelag is nested within Norway,  
- Hesse is nested within Germany
- Germany and Norway are nested within Europe

# Nested Effects leading to nests

Linnean classification is another example:

Animalia - Chordata - Aves - Psittaciformes - Psittacidae - Eclectus - E. roratus - E. r. solomonensis - Eric

There is a grey parrot called Eric, but he is not the same

# Nestedness

Take a model with A/B

Each level of B can only be in one level of A

# Fitting Nested Models

```{r BirdBrains}
data <- read.csv("../Data/BirdBrains.csv")
data.C <- data[grep('^[G]', data$Order),]
phylo <- lm(logBrainMass~Order/Family, data=data.C)
summ <- summary(phylo)$coefficients

round(summ[!is.na(summ[,"Estimate"]),],2)
```

No OrderGalliformes:FamilyGruidae effect!

# Next Week

Model Selection!
