Appendix G — Solutions to section questions

Below are the solutions to exercises within sections, given here (which I appreciate is a weird place for them to be in the handout) so you don’t ‘accidentally’ see them. Not every answer is given here; I don’t need to give you answers to everything and so I don’t. Please come to the classes and ask me questions; you also have a dedicated question-and-answer section to go through questions. Come to those and engage with me: my experience of giving classes all the answers is they read the answers, convince themselves they could have answered the questions, and then fail the exam. Consider also that, in almost every case, the answers to the questions are in your handout: you are copy-pasting code from the handout and then reading the answers. This is by design and isn’t cheating! Thus reading me doing the same won’t help you as much as speaking to me will…

G.1 Philosophy of statistics

Fundamentals of probability definitions.

There is only one tomorrow but there are many volcanoes, which may or may not erupt on one day, so we can construct a frequentist probability. Equally, there are many days even if there were only one volcano.
Blaargian probability: $0.1 \times 0.9 = 0.09$. Floogian probability: $0.4 \times 0.6 = 0.24$. Thus the Floogian is most likely because it has the highest probability.
Two parameters are more than one, and we prefer simpler models (fewer parameters). Thus this makes you a little more likely to give credance to the Floogian distribution, but we don’t have a way to contrast likelihood/probabiltiy with model complexity (…yet…).
As you change the mean ($\mu$) to be more positive, the mean shifts to the right of the histogram.
hist(rnorm(n=10000, mean=0, sd=1)) and hist(rnorm(n=10000, mean=0, sd=10)) would be sufficient to show you that higher standard deviation (and so variance) are ‘flatter’ or ‘more spread out’ histograms.
hist(rpois(n=10000, lambda=1)) and hist(rpois(n=10000, lambda=30)) would be sufficient to show you that Poisson is quite skewed (long ‘tail’ hanging to the right of the distribution) and that this skew reduces as the value of $\lambda$ increases (indeed, it starts looking/behaving like a normal distribution somewhat).
Degrees of freedom and number of observations are related as outlined in the handout—estimating a parameter ‘soaks up’ a degree of freedom, and the total number ‘available’ is the number of observations. Playing with these parameters somewhat changes the degree of skew, as was seen with the Poisson.
I can’t really say any more here than is in the question—you’ve got the answer and the code for it there.
This answer is also left deliberately blank!

Maximum likelihood by hand.

prod(dnorm(c(120,130,110,105,121), mean=120, sd=5)).
95% CIs are 2 log-likelihood units away from each other, so we can try different values to see what’s within 2 units away from the maximum value. There’s no need to be too precise about it, so testing something in the following way would likely work:

data <- c(120,130,110,105,121)
sd <- sd(data)
log(prod(dnorm(data, mean=100, sd=sd)))
#> [1] -25.67212
log(prod(dnorm(data, mean=102, sd=sd))) # 2 log-units from MLE
#> [1] -23.99684
log(prod(dnorm(data, mean=105, sd=sd)))
#> [1] -21.87171
log(prod(dnorm(data, mean=110, sd=sd)))
#> [1] -19.36395
log(prod(dnorm(data, mean=115, sd=sd))) # MLE is around here
#> [1] -18.14886
log(prod(dnorm(data, mean=120, sd=sd)))
#> [1] -18.22641
log(prod(dnorm(data, mean=122, sd=sd)))
#> [1] -18.61938
log(prod(dnorm(data, mean=125, sd=sd))) # 2 log-units from MLE
#> [1] -19.59663
log(prod(dnorm(data, mean=130, sd=sd)))
#> [1] -22.25951

…of course this is all approximate. Note that you could also do the log-ing from within dnorm with something like sum(dnorm(data, mean=115, sd=sd, log=TRUE)). c. As above, but now with the SD varying.

data <- c(120,130,110,105,121)
mean <- mean(data)
log(prod(dnorm(data, mean=mean, sd=4)))
#> [1] -23.61366
log(prod(dnorm(data, mean=mean, sd=5))) # 2 log-units from MLE
#> [1] -20.37788
log(prod(dnorm(data, mean=mean, sd=6)))
#> [1] -18.92571
log(prod(dnorm(data, mean=mean, sd=8))) 
#> [1] -18.01378
log(prod(dnorm(data, mean=mean, sd=9))) # MLE is around here
#> [1] -17.96847
log(prod(dnorm(data, mean=mean, sd=10)))
#> [1] -18.04162
log(prod(dnorm(data, mean=mean, sd=12)))
#> [1] -18.36228
log(prod(dnorm(data, mean=mean, sd=15)))
#> [1] -18.9945
log(prod(dnorm(data, mean=mean, sd=20)))
#> [1] -20.05685
log(prod(dnorm(data, mean=mean, sd=25))) # 2 log-units from MLE
#> [1] -20.99851

…of course this is all approximate. d.

data <- c(120,130,110,105,121)
mean <- mean(data)
sd <- sd(data)
prnorm(100, mean=mean, sd=sd)

1-pnorm(150, mean=mean, sd=sd). You should absolutely spend less time with Michael, and inform the Senior Tutor about him.

You will need to draw a tree diagram with at least three ‘levels’ and so it should have eight terminal nodes (leaves). Check your answer by checking the probabilities you get with the answeres below.
Use the diagram from (a) above, making sure to sum up across all the terminal nodes (outcomes) that match what Sarah’s looking for (two reds).
For repeating (a), you could get this by recognising that it has to be a red each time and so calculating $\frac{2}{5}^3 = \frac{8}{125} = 0.064$. Note that giving fraction and not decimal answers would be fine in an exam. For repeating (b), use the Binomial formula: $\binom{3}{2} \times \frac{2}{5}^2 \times \frac{3}{5}^1 = 3 \times \frac{4}{25} \times \frac{3}{5} = \frac{3 \times 4 \times 3}{25 \times 5} = \frac{36}{125} = 0.288$.
dbinom(3, 3, .4) and dbinom(2, 3, .4) respectively.
You can get this lots of ways. sum(dbinom(30:50, 200, 30/170)) or pbinom(50, 200, 30/170) - pbinom(29, 200, 30/170) all work (note that it’s ‘29’ in the second answer here because otherwise you will remove all the 30 examples). The critical thing is, in the first, we’re summing up all the probabilities between 30 and 50, and in the second we’re taking the probabilities of less than 50 and subtracting from that 29 and below.
Doesn’t affect it: all that matters is that the ball isn’t red. We have defined ‘success’ in this trial as a red ball, so it doesn’t matter. There is something called the ‘multinomial distribution’, which generalised the binomial to more than one outcome, but you only need to use that if you’re interested in more than one kind of outcome.
The ‘with replacement’ is critical here because it means the balls are put back in the bag after each draw. Otherwise, you would have to re-calculate each time (…the number of red balls left in the bag would change, and so your probabilities would too) and you couldn’t use the Binomial distribution in this way. You could do this with tree diagrams but, if you were trying to do part (e), it would get very involved because you’ve need a tree 20 nodes deep. And now you know the origin story of Michael: once a promising student, driven to madness by his excitement with drawing from distributions without replacement.

G.2 Test statistics

Fundamentals of probability definitions.

The code is given to you. Look at the plots and eyeball them. I don’t see much of a sex effect (perhaps more variance in the males), but it does look like ‘they’ (tourists) caught bigger traps than ‘us’ (me and Helen¹).
with(data, t.test(diameter ~ sex)). We can see that the p-value is greater than 0.05 (5%), which suggests that, by chance, we would expect a test-statistic (in this case t value) this extreme more than 5% of the time. We’re not surprised by the difference we see among the variables. If you’re someone who likes talking about accepting/rejecting null hypotheses I’m fine with you describing this like that, and talking about ‘no significant difference’. Notice all we’ve done is replaced the word ‘boxplot’ with ‘t.test’ in the code.
with(data, t.test(diameter ~ source)). Now the p-value is less than 5%, and so we are surprised by this result: the difference is statistically significant and we reject the null hypothesis that there is no difference among the size of crabs caught by us and them.
Think about the biology here. Anything along the lines of ‘if they catch all the big crabs then there might not be any for other to catch’, ‘big crabs are older crabs and so they produce all the young so the population might die out’, ‘the crabs might evolve to be smaller and so less fun to catch’—anything that doesn’t sound absurd but shows you can vaguely think about biology—would be acceptable.

Use the code given in the question. It sure looks like they’re correlated, doesn’t it?
You could run with(data, summary(lm(mpg ~ acceleration))) or summary(lm(mpg ~ acceleration, data=data)) here and get the same result. You could also save out the model model <- lm(mpg ~ accelerate, data=data) and then summarise it summary(model) as we will be doing later in the course. At any rate, there is a positive slope value (1.20 or so) and the p-value is less than 5% (NOTE: the p-value we care about is not in the coefficients table on the same row as ‘acceleration’ but is at the bottom-right of the table—that is a contrast coefficient and is something different). See the answer above for the spiel about a p-value that is significant or surprising for how you could summarise this.
You can also run an ANOVA in different ways, as outlined in the answer above, but something like summary(lm(mpg ~ cylinders, data=data)) should give you output that shows that, as with the last example, there is an impact of having a cylinder. Remember that with ANOVA we care about the elephant in the room: the average MPG of something with less than or equal to 4 cylinder is 29.1, and the average of those with more than 4 is $29.1 + -11.8 = 29.1 - 11.8 = 17.3$.
Fit a t-test with with(data, t.test(mpg ~ cylinders)). I won’t go into the details of the extension footnote but, in the case that there are only two treatment levels, an ANOVA and a t-test are identical (note that there are many different kinds of t-test, so this isn’t always strictly true).
If you drove 10000 miles a year, then (10000/mean(data$mpg)) * 8887 would give you the number of grams of CO$_2$ you emit. (10000/(mean(data$mpg)+10)) * 8887 gives you the better car.
Hey, let’s do this using some intermediary variables (even though you can do this in a one-liner):

average <- (10000/mean(data$mpg)) * 8887
new <- (10000/(mean(data$mpg)+10)) * 8887
(average - new) * 2.54

…there we go. Apologies for using dollars in this exercise, but given the current political situation it’s the only way I can run this exercise without updating it every few days². 3.

my.t.test <- function(x, y){
  observation <- mean(x) - mean(y)
  std.err <- sqrt((sd(x)^2)/length(x) + (sd(y)^2)/length(y))
  t <- observation / expectation
  return(t)
}
calc.sst <- function(x){
  mean <- mean(x)
  squared.diffs <- (mean - x)^2
  return(sum(squared.diffs))
}
calc.ssm <- function(resp, explan){
  model <- lm(resp ~ explan)
  sst <- calc.sst(resp)
  residuals <- residuals(model)
  sse <- sum((residuals)^2)
  ssm <- sst - sse
  return(ssm)
}
# (If you have completed the extension section, try the following as well)
r2 <- function(resp, explan){
  sst <- calc.sst(resp)
  sse <- calc.sse(resp, explan)
  r2 <- ssm / sst
  return(r2)
}

summary(lm(diameter ~ source, data=crabs)) and summary(lm(mpg ~ acceleration, data=mileage)) show me that the $r^2$ of the mileage data is better. So, I suppose, that model is doing better—they have the same number of parameters, after all.
Neither of the above models have leverage plots (the last of the four) that suggest any particular outliers are having undue effects on model fit. The mileage plot is doing the best of the two on the Q-Q plot (i.e., has the most normally-distributed residuals), although neither Q-Q plot is something I would lose a lot of sleep over (the mileage plot is probably as bad as it could be without me being concerned—note that there’s no evidence of other issues, though, which the model has in its favour, along with its simplicity).
The crab data are plotted against a categorical explanatory variable in the above model, and so the predictions are essentially just two values (‘them’ or ‘us’) and so we have the strange straight lines in the plots accordingly.

G.3 Multiple models

Fundamentals of probability definitions.

summary(lm(rating ~ year, data=data)): yes, there is a positive slope (and a significant p-value of the overall modell).
summary(lm(rating ~ year + kickstarter, data=data)): yes, the contrast of kickstarter is positive.
You need to do this in steps, comparing your alternative model (the interaction) with the null (the additive model). So something like:

interaction <- lm(rating ~ year + kickstarter, data=data)
additive <- lm(rating ~ year + kickstarter, data=data)
anova(interaction, additive)

…which tells us there isn’t an interaction here. d. Again, remember to set up a sensible alternative model and a sensitive null as well. Something like:

alternative <- lm(rating ~ year + kickstarter + co.op + complexity, data=data)
null <- lm(rating ~ year + kickstarter, data=data)
anova(alternative, null) #1
summary(alternative)     #2

…which tells us that the additional variance explained by these two terms is more than we would expect by chance (which we can see in line-highlight 1), and from looking at the model output we can see that Michael is only partially correct: board gamers like more complex games, but they also prefer co-operative games and so don’t only care about winning. 2. a. Heavens if I know what’s going on here: there are too many variables floating around. Note that the I(year^2) nonsense fits a squared (quadratic; curved) term to the model. No, really, that’s the answer here: summary(top.model) is interpretable but heaven help you figuring it all out. We need to summarise this information somehow… b. The following code will run everything through for you, and show you a few things like games that more people own are rated higher (…makes sense…), and weightier (more complex) games are very popular.

models <- dredge(top.model)
summary(model.avg(models, subset=delta<4, fit=TRUE))
sw(model.avg(models, subset=delta<4, fit=TRUE))

The following code will make things even easier by giving us standard effect sizes. This makes it clear that the major drivers are whether or not it’s a recent game that is complex and owned by a lot of people (possible circularity in that last one).

# You could scale each column manually by hand, or do it faster like this
s.data <- data
for(i in seq_along(names(s.data)))
  s.data[,i] <- as.numeric(scale(s.data[,i],scale=sd(s.data[,i])*2))
s.data$rating <- data$rating # No need to re-scale the rating; let's put it back
top.model <- lm(rating ~ year + I(year^2) + max_players +
log.play.time + two.players + three.players + weight +
mech.pca.1 + mech.pca.2 + mech.pca.3 +  kickstarter +
log.owned + log.n.ratings + co.op,
data=s.data, na.action="na.pass")
# ... look above - do you see why using data= statements is useful now?
models <- dredge(top.model)
summary(model.avg(models, subset=delta<4, fit=TRUE))

back.transform <- function(coefficient, orig.x){
  sd.orig <- sd(x)
  mean.orig <- mean(x)
  rescaled <- coefficient * sd.orig + mean.orig
  return(rescaled)
}

The residuals aren’t unimodal (show a single mean), let alone normally distributed, so we have a problem here. Also the $r^2$ isn’t looking very clever.

model <- lm(damaged ~ diameter + sex, data=data)
hist(residuals(model))

model <- glm(damaged ~ diameter + sex, data=data, family=binomial). Bigger crabs are more likely to be damaged, as are male crabs.

model <- glm(damaged ~ diameter + sex, data=data, family=binomial)
pred <- with(data,
  data.frame(diameter=seq(min(diameter), max(diameter)), sex="female")
) # change the above line is you want to see the effect for males
pred$damaged <- predict(model, pred, type="response")
with(data, boxplot(diameter ~ damaged, horizontal=TRUE))
with(pred, lines(diameter, damaged+1, col="red", lwd=3))

model <- lm(species ~ latitude + elevation + habitat, data=data). Fewer ants in lower latitudes and lower elevations; more ants in forests.
model <- glm(species ~ latitude + elevation + habitat, data=data, family=poisson). Much the same as above!

model <- glm(species ~ latitude + elevation + habitat,
  data=data, family=poisson)
pred <- with(data,
  data.frame(
    latitude=seq(min(latitude), max(latitude), by=.1),
    elevation=mean(elevation), habitat="forest"
  )
)
pred$species <- predict(model, pred, type="response")
with(data, plot(species ~ latitude))
with(pred, lines(species ~ latitude, col="red", lwd=3))

G.4 Hierarchical (Bayesian) models

You read everything? Great.

model <- lm(human ~ temp.change + trump.percent, data=data)
summary(model)
null <- lm(human ~ trump.percent, data=data)
anova(model, null)

More likely to believe if temperature has changed more (quick model comparison below tells us the difference isn’t significant), seems like voting for Trump means you’re less likely. Likely not a reliable model, though, because we’re pseudo-replicating across states (the voting data are the same for every county within a state). b.

library(lme4)
model <- lmer(human ~ temp.change + trump.percent + (1|State_abb),
  data=data)
summary(model)
# No p-values, so let's load the package... sigh...
library(lmerTest)
summary(model)

…same results as before, if anything with larger effect sizes, and with the same significances. c.

state.cols <- setNames(rainbow(length(unique(data$State_abb))),
  unique(data$State_abb))
with(data, plot(human ~ trump.percent, col=state.cols[State_abb], pch=20))
for(i in seq)
abline(coef(model))

Who is an epidemiologist now, I believe, and I imagine is not still talking about these f**king crabs quite some years later…↩︎
At the time of writing our Prime Minister is Liz Truss, and everything is apparently Putin’s fault… Update: The previous statement has been left as a time-capsule…↩︎

# Solutions to section questions Below are the solutions to exercises within sections, given here (which I appreciate is a weird place for them to be in the handout) so you don't 'accidentally' see them. Not every answer is given here; I don't need to give you answers to everything and so I don't. Please come to the classes and ask me questions; you also have a dedicated question-and-answer section to go through questions. Come to those and engage with me: my experience of giving classes all the answers is they read the answers, convince themselves they could have answered the questions, and then fail the exam. Consider also that, in almost every case, the answers to the questions are in your handout: you are copy-pasting code from the handout and then reading the answers. This is by design and isn't cheating! Thus reading me doing the same won't help you as much as speaking to me will... ## Philosophy of statistics 1. **Fundamentals of probability definitions.** a. There is only one tomorrow but there are many volcanoes, which may or may not erupt on one day, so we can construct a frequentist probability. Equally, there are many days even if there were only one volcano. b. Blaargian probability: $0.1 \times 0.9 = 0.09$. Floogian probability: $0.4 \times 0.6 = 0.24$. Thus the Floogian is most likely because it has the highest probability. c. Two parameters are more than one, and we prefer simpler models (fewer parameters). Thus this makes you a little more likely to give credance to the Floogian distribution, but we don't have a way to contrast likelihood/probabiltiy with model complexity (...yet...). d. As you change the mean ($\mu$) to be more positive, the mean shifts to the right of the histogram. e. `hist(rnorm(n=10000, mean=0, sd=1))` and `hist(rnorm(n=10000, mean=0, sd=10))` would be sufficient to show you that higher standard deviation (and so variance) are 'flatter' or 'more spread out' histograms. f. `hist(rpois(n=10000, lambda=1))` and `hist(rpois(n=10000, lambda=30))` would be sufficient to show you that Poisson is quite skewed (long 'tail' hanging to the right of the distribution) and that this skew reduces as the value of $\lambda$ increases (indeed, it starts looking/behaving like a normal distribution somewhat). g. Degrees of freedom and number of observations are related as outlined in the handout—estimating a parameter 'soaks up' a degree of freedom, and the total number 'available' is the number of observations. Playing with these parameters somewhat changes the degree of skew, as was seen with the Poisson. h. I can't really say any more here than is in the question—you've got the answer and the code for it there. i. This answer is also left deliberately blank! 2. **Maximum likelihood by hand.** a. `prod(dnorm(c(120,130,110,105,121), mean=120, sd=5))`. b. 95% CIs are 2 log-likelihood units away from each other, so we can try different values to see what's within 2 units away from the maximum value. There's no need to be too precise about it, so testing something in the following way would likely work: ```{r} data <- c(120,130,110,105,121) sd <- sd(data) log(prod(dnorm(data, mean=100, sd=sd))) log(prod(dnorm(data, mean=102, sd=sd))) # 2 log-units from MLE log(prod(dnorm(data, mean=105, sd=sd))) log(prod(dnorm(data, mean=110, sd=sd))) log(prod(dnorm(data, mean=115, sd=sd))) # MLE is around here log(prod(dnorm(data, mean=120, sd=sd))) log(prod(dnorm(data, mean=122, sd=sd))) log(prod(dnorm(data, mean=125, sd=sd))) # 2 log-units from MLE log(prod(dnorm(data, mean=130, sd=sd))) ``` ...of course this is all approximate. Note that you could also do the log-ing from within `dnorm` with something like `sum(dnorm(data, mean=115, sd=sd, log=TRUE))`. c. As above, but now with the SD varying. ```{r} data <- c(120,130,110,105,121) mean <- mean(data) log(prod(dnorm(data, mean=mean, sd=4))) log(prod(dnorm(data, mean=mean, sd=5))) # 2 log-units from MLE log(prod(dnorm(data, mean=mean, sd=6))) log(prod(dnorm(data, mean=mean, sd=8))) log(prod(dnorm(data, mean=mean, sd=9))) # MLE is around here log(prod(dnorm(data, mean=mean, sd=10))) log(prod(dnorm(data, mean=mean, sd=12))) log(prod(dnorm(data, mean=mean, sd=15))) log(prod(dnorm(data, mean=mean, sd=20))) log(prod(dnorm(data, mean=mean, sd=25))) # 2 log-units from MLE ``` ...of course this is all approximate. d. ```{r} #| eval: false data <- c(120,130,110,105,121) mean <- mean(data) sd <- sd(data) prnorm(100, mean=mean, sd=sd) ``` e. `1-pnorm(150, mean=mean, sd=sd)`. You should absolutely spend less time with Michael, and inform the Senior Tutor about him. 3. a. You will need to draw a tree diagram with at least three 'levels' and so it should have eight terminal nodes (leaves). Check your answer by checking the probabilities you get with the answeres below. b. Use the diagram from (a) above, making sure to sum up across all the terminal nodes (outcomes) that match what Sarah's looking for (two reds). c. For repeating (a), you could get this by recognising that it has to be a red each time and so calculating $\frac{2}{5}^3 = \frac{8}{125} = 0.064$. Note that giving fraction and not decimal answers would be fine in an exam. For repeating (b), use the Binomial formula: $\binom{3}{2} \times \frac{2}{5}^2 \times \frac{3}{5}^1 = 3 \times \frac{4}{25} \times \frac{3}{5} = \frac{3 \times 4 \times 3}{25 \times 5} = \frac{36}{125} = 0.288$. d. `dbinom(3, 3, .4)` and `dbinom(2, 3, .4)` respectively. e. You can get this lots of ways. `sum(dbinom(30:50, 200, 30/170))` or `pbinom(50, 200, 30/170) - pbinom(29, 200, 30/170)` all work (note that it's '29' in the second answer here because otherwise you will remove all the 30 examples). The critical thing is, in the first, we're summing up all the probabilities between 30 and 50, and in the second we're taking the probabilities of less than 50 and subtracting from that 29 and below. f. Doesn't affect it: all that matters is that the ball isn't red. We have defined 'success' in this trial as a red ball, so it doesn't matter. There is something called the 'multinomial distribution', which generalised the binomial to more than one outcome, but you only need to use that if you're interested in more than one kind of outcome. g. The 'with replacement' is critical here because it means the balls are put back in the bag after each draw. Otherwise, you would have to re-calculate each time (...the number of red balls left in the bag would change, and so your probabilities would too) and you couldn't use the Binomial distribution in this way. You could do this with tree diagrams but, if you were trying to do part (e), it would get very involved because you've need a tree 20 nodes deep. And now you know the origin story of Michael: once a promising student, driven to madness by his excitement with drawing from distributions without replacement. ## Test statistics 1. **Fundamentals of probability definitions.** a. The code is given to you. Look at the plots and eyeball them. I don't see much of a sex effect (perhaps more variance in the males), but it does look like 'they' (tourists) caught bigger traps than 'us' (me and Helen^[Who is an epidemiologist now, I believe, and I imagine is not still talking about these f**king crabs quite some years later...]). b. `with(data, t.test(diameter ~ sex))`. We can see that the p-value is greater than 0.05 (5%), which suggests that, by chance, we would expect a test-statistic (in this case t value) this extreme more than 5% of the time. We're not surprised by the difference we see among the variables. If you're someone who likes talking about accepting/rejecting null hypotheses I'm fine with you describing this like that, and talking about 'no significant difference'. Notice all we've done is replaced the word 'boxplot' with 't.test' in the code. c. `with(data, t.test(diameter ~ source))`. Now the p-value is less than 5%, and so we are surprised by this result: the difference is statistically significant and we reject the null hypothesis that there is no difference among the size of crabs caught by us and them. d. Think about the biology here. Anything along the lines of 'if they catch all the big crabs then there might not be any for other to catch', 'big crabs are older crabs and so they produce all the young so the population might die out', 'the crabs might evolve to be smaller and so less fun to catch'—anything that doesn't sound absurd but shows you can vaguely think about biology—would be acceptable. 2. a. Use the code given in the question. It sure looks like they're correlated, doesn't it? b. You could run `with(data, summary(lm(mpg ~ acceleration)))` or `summary(lm(mpg ~ acceleration, data=data))` here and get the same result. You could also save out the model `model <- lm(mpg ~ accelerate, data=data)` and then summarise it `summary(model)` as we will be doing later in the course. At any rate, there is a positive slope value (1.20 or so) and the p-value is less than 5% (NOTE: the p-value we care about is **not** in the coefficients table on the same row as 'acceleration' but is at the bottom-right of the table—that is a contrast coefficient and is something different). See the answer above for the spiel about a p-value that is significant or surprising for how you could summarise this. c. You can also run an ANOVA in different ways, as outlined in the answer above, but something like `summary(lm(mpg ~ cylinders, data=data))` should give you output that shows that, as with the last example, there is an impact of having a cylinder. Remember that with ANOVA we care about **the elephant in the room**: the average MPG of something with less than or equal to 4 cylinder is 29.1, and the average of those with more than 4 is $29.1 + -11.8 = 29.1 - 11.8 = 17.3$. d. Fit a t-test with `with(data, t.test(mpg ~ cylinders))`. I won't go into the details of the extension footnote but, in the case that there are only two treatment levels, an ANOVA and a t-test are identical (note that there are many different kinds of t-test, so this isn't always strictly true). e. If you drove 10000 miles a year, then `(10000/mean(data$mpg)) * 8887` would give you the number of grams of CO$_2$ you emit. `(10000/(mean(data$mpg)+10)) * 8887` gives you the better car. f. Hey, let's do this using some intermediary variables (even though you can do this in a one-liner): ```{r} #| eval: false average <- (10000/mean(data$mpg)) * 8887 new <- (10000/(mean(data$mpg)+10)) * 8887 (average - new) * 2.54 ``` ...there we go. Apologies for using dollars in this exercise, but given the current political situation it's the only way I can run this exercise without updating it every few days^[At the time of writing our Prime Minister is Liz Truss, and everything is apparently Putin's fault... Update: The previous statement has been left as a time-capsule...]. 3. ```{r} #| eval: false my.t.test <- function(x, y){ observation <- mean(x) - mean(y) std.err <- sqrt((sd(x)^2)/length(x) + (sd(y)^2)/length(y)) t <- observation / expectation return(t) } calc.sst <- function(x){ mean <- mean(x) squared.diffs <- (mean - x)^2 return(sum(squared.diffs)) } calc.ssm <- function(resp, explan){ model <- lm(resp ~ explan) sst <- calc.sst(resp) residuals <- residuals(model) sse <- sum((residuals)^2) ssm <- sst - sse return(ssm) } # (If you have completed the extension section, try the following as well) r2 <- function(resp, explan){ sst <- calc.sst(resp) sse <- calc.sse(resp, explan) r2 <- ssm / sst return(r2) } ``` 4. a. `summary(lm(diameter ~ source, data=crabs))` and `summary(lm(mpg ~ acceleration, data=mileage))` show me that the $r^2$ of the mileage data is better. So, I suppose, that model is doing better—they have the same number of parameters, after all. b. Neither of the above models have leverage plots (the last of the four) that suggest any particular outliers are having undue effects on model fit. The mileage plot is doing the best of the two on the Q-Q plot (*i.e.*, has the most normally-distributed residuals), although neither Q-Q plot is something I would lose a lot of sleep over (the mileage plot is probably as bad as it could be without me being concerned—note that there's no evidence of other issues, though, which the model has in its favour, along with its simplicity). c. The crab data are plotted against a categorical explanatory variable in the above model, and so the predictions are essentially just two values ('them' or 'us') and so we have the strange straight lines in the plots accordingly. ## Multiple models 1. **Fundamentals of probability definitions.** a. `summary(lm(rating ~ year, data=data))`: yes, there is a positive slope (and a significant p-value of the overall modell). b. `summary(lm(rating ~ year + kickstarter, data=data))`: yes, the contrast of kickstarter is positive. c. You need to do this in steps, comparing your alternative model (the interaction) with the null (the additive model). So something like: ```{r} #| eval: false interaction <- lm(rating ~ year + kickstarter, data=data) additive <- lm(rating ~ year + kickstarter, data=data) anova(interaction, additive) ``` ...which tells us there isn't an interaction here. d. Again, remember to set up a sensible alternative model and a sensitive null as well. Something like: ```{r} #| eval: false alternative <- lm(rating ~ year + kickstarter + co.op + complexity, data=data) null <- lm(rating ~ year + kickstarter, data=data) anova(alternative, null) #1 summary(alternative) #2 ``` ...which tells us that the additional variance explained by these two terms is more than we would expect by chance (which we can see in line-highlight 1), and from looking at the model output we can see that Michael is only partially correct: board gamers like more complex games, but they also prefer co-operative games and so don't only care about winning. 2. a. Heavens if I know what's going on here: there are too many variables floating around. Note that the `I(year^2)` nonsense fits a squared (quadratic; curved) term to the model. No, really, that's the answer here: `summary(top.model)` is interpretable but heaven help you figuring it all out. We need to summarise this information somehow... b. The following code will run everything through for you, and show you a few things like games that more people own are rated higher (...makes sense...), and weightier (more complex) games are very popular. ```{r} #| eval: false models <- dredge(top.model) summary(model.avg(models, subset=delta<4, fit=TRUE)) sw(model.avg(models, subset=delta<4, fit=TRUE)) ``` c. The following code will make things even easier by giving us standard effect sizes. This makes it clear that the major drivers are whether or not it's a recent game that is complex and owned by a lot of people (possible circularity in that last one). ```{r} #| eval: false # You could scale each column manually by hand, or do it faster like this s.data <- data for(i in seq_along(names(s.data))) s.data[,i] <- as.numeric(scale(s.data[,i],scale=sd(s.data[,i])*2)) s.data$rating <- data$rating # No need to re-scale the rating; let's put it back top.model <- lm(rating ~ year + I(year^2) + max_players + log.play.time + two.players + three.players + weight + mech.pca.1 + mech.pca.2 + mech.pca.3 + kickstarter + log.owned + log.n.ratings + co.op, data=s.data, na.action="na.pass") # ... look above - do you see why using data= statements is useful now? models <- dredge(top.model) summary(model.avg(models, subset=delta<4, fit=TRUE)) ``` 3. ```{r} #| eval: false back.transform <- function(coefficient, orig.x){ sd.orig <- sd(x) mean.orig <- mean(x) rescaled <- coefficient * sd.orig + mean.orig return(rescaled) } ``` 4. a. The residuals aren't unimodal (show a single mean), let alone normally distributed, so we have a problem here. Also the $r^2$ isn't looking very clever. ```{r} #| eval: false model <- lm(damaged ~ diameter + sex, data=data) hist(residuals(model)) ``` b. `model <- glm(damaged ~ diameter + sex, data=data, family=binomial)`. Bigger crabs are more likely to be damaged, as are male crabs. ```{r} #| eval: false model <- glm(damaged ~ diameter + sex, data=data, family=binomial) pred <- with(data, data.frame(diameter=seq(min(diameter), max(diameter)), sex="female") ) # change the above line is you want to see the effect for males pred$damaged <- predict(model, pred, type="response") with(data, boxplot(diameter ~ damaged, horizontal=TRUE)) with(pred, lines(diameter, damaged+1, col="red", lwd=3)) ``` 5. a. `model <- lm(species ~ latitude + elevation + habitat, data=data)`. Fewer ants in lower latitudes and lower elevations; more ants in forests. b. `model <- glm(species ~ latitude + elevation + habitat, data=data, family=poisson)`. Much the same as above! c. ```{r} #| eval: false model <- glm(species ~ latitude + elevation + habitat, data=data, family=poisson) pred <- with(data, data.frame( latitude=seq(min(latitude), max(latitude), by=.1), elevation=mean(elevation), habitat="forest" ) ) pred$species <- predict(model, pred, type="response") with(data, plot(species ~ latitude)) with(pred, lines(species ~ latitude, col="red", lwd=3)) ``` ## Hierarchical (Bayesian) models 1. You read everything? Great. 2. a. ```{r} #| eval: false model <- lm(human ~ temp.change + trump.percent, data=data) summary(model) null <- lm(human ~ trump.percent, data=data) anova(model, null) ``` More likely to believe if temperature has changed more (quick model comparison below tells us the difference isn't significant), seems like voting for Trump means you're less likely. Likely not a reliable model, though, because we're pseudo-replicating across states (the voting data are the same for every county within a state). b. ```{r} #| eval: false library(lme4) model <- lmer(human ~ temp.change + trump.percent + (1|State_abb), data=data) summary(model) # No p-values, so let's load the package... sigh... library(lmerTest) summary(model) ``` ...same results as before, if anything with larger effect sizes, and with the same significances. c. ```{r} #| eval: false state.cols <- setNames(rainbow(length(unique(data$State_abb))), unique(data$State_abb)) with(data, plot(human ~ trump.percent, col=state.cols[State_abb], pch=20)) for(i in seq) abline(coef(model)) ```