options(digits = 4, scipen = 999)
knitr::opts_chunk$set(echo = TRUE, fig.width = 4, fig.height = 4)
library(tidyverse)
library(lmSupport)
library(psych)
library(car)
library(pwr)
library(effects)
library(lme4)
library(kableExtra)
options(knitr.kable.NA = '')

Updated to fix a typo

1

diabetes <- read.csv("https://whlevine.hosted.uark.edu/psyc5143/diabetes.csv")

# a: assessing nonlinearity
ggplot(data = diabetes, aes(y = BMI, x = Age)) +
    geom_point() +
    geom_smooth(method = "loess", se = FALSE)

# b: model
model1b <- lm(BMI ~ Age, diabetes)

# c: quadratic model
model1c <- lm(BMI ~ Age + I(Age^2), diabetes)

# d: mean-centered quadratic model
diabetes <- diabetes %>% 
    mutate(age.c = Age - mean(Age))

model1d <- lm(BMI ~ age.c + I(age.c^2), diabetes)

# e: linear slope of age when age = 60

# center age at 60 and fit the model!
diabetes <- diabetes %>% 
    mutate(age60 = Age - 60)

model1e <- lm(BMI ~ age60 + I(age60^2), diabetes)

See the graph below. There is a pretty clear non-linear component.

See the end of the answer key for some bonus information on how to quantify non-linearity. This is something I never covered in class.

The slope is (about) 0.026, which means that BMI is (according to this model) expected to increase by 0.026 per year of aging.
The linear slope is about 0.636, which means that at age 0, the quadratic function has a very strong positive slope at that point; because the sample includes no one younger than about 21, this is essentially meaningless. The quadratic slope is -0.008 or so, which means that the relationship between age and BMI is expected to get less positive as one ages (and eventually become negative).
The linear slope is about 0.123, which means that at age 33.4 or so (i.e., the mean age), the quadratic function still has a positive slope. The intercept is the predicted BMI for someone at this age. And the quadratic slope means the same thing here that it did in part c. The linear slope and the intercept are illustrated below.

# this may be useful for those who want to see these parameters

# first, store the parameter estimates for the tangent line
slope     <- coef(model1d)[2]
intercept <- coef(model1d)[1] + coef(model1d)[2]*(0 - mean(diabetes$Age))

ggplot(data = diabetes, aes(x = Age, y = BMI)) +
    geom_point(size = 0.5, color = "gray") +
    stat_smooth(method = "lm", formula = y ~ x + I(x^2), se = F) +
    geom_vline(xintercept = mean(diabetes$Age)) +
    geom_abline(slope = slope, intercept = intercept, color = "red") +
    theme_bw()

The linear slope at 60 years of age is about -0.285. BMI is expected to be going down as one ages at this point.
The linear slope in part b is the overall linear slope of the age-BMI relationship, in some sense the average of the slope across all values of age. By contrast, the linear slope in part d is restricted to only the mean of age in its interpretation; it’s a point slope because it’s the slope for only a particular value of age.

2

trash <- read.csv("https://whlevine.hosted.uark.edu/psyc5143/trash.csv", stringsAsFactors = TRUE)

# checking default coding
contrasts(trash$thoughts)
contrasts(trash$treatment) # control is baseline

# changing to contrast coding for the thoughts effect (will facilitate the
# second and third hypotheses)
contrasts(trash$thoughts) <- c(-1/2, 1/2)

# testing the interaction and the effect of thoughts in the control condition
model.2 <- lm(body.attitude ~ thoughts*treatment, trash)
summary(model.2)

# flipping the dummy coding around to test the effect of thoughts in the
# disposal condition and re-fitting the model
contrasts(trash$treatment) <- c(1, 0)
model.2.v2 <- lm(body.attitude ~ thoughts*treatment, trash)
summary(model.2.v2)

# some descriptive stats
trash %>% group_by(thoughts, treatment) %>% summarise(M = mean(body.attitude))

A linear model with thoughts (positive, negative), treatment (control, disposal), and their interaction as predictors of body attitude was fit. With respect to the first hypothesis, the interaction was significant, t(76) = 2.76, p = .007. In the control condition, there was a significant effect of thoughts, t(76) = 3.68, p = .0004, such that negative thoughts (M = 5.4) were associated with poorer body attitude than were positive thoughts (M = 6.4). By contrast, in the thought disposal condition, there was not a significant difference in body attitude between the negative thoughts condition (M = 5.6) and the positive thoughts condition (M = 5.6), t(76) = 0.21, p = .83.

3

rewards <- read.csv("https://whlevine.hosted.uark.edu/psyc5143/rewards.csv")

# I can't deal with these numbers for conditions, so I'm creating a factor version
rewards$condF <- factor(rewards$cond, labels = c("control1", "control2", "hungry", "thirsty", "sad"))

# b and e: creating contrasts,  verifying orthogonality, fitting the model
rewards <- rewards %>% 
  mutate(
    H1 = case_when(condF == "control1" ~  1/2,
                                 condF == "control2" ~ -1/2,
                                 TRUE ~ 0),
    H2 = ifelse(cond < 3, 3/5, -2/5),
    H3 = case_when(condF == "hungry"   ~  1/3,
                                 condF == "thirsty"  ~  1/3,
                                 condF == "sad"      ~ -2/3,
                                 TRUE ~ 0),
    H4 = case_when(condF == "hungry"   ~  1/2,
                                 condF == "thirsty"  ~ -1/2,
                                 TRUE ~ 0)
  )

cor(rewards[4:7])

model.3e <- lm(trials ~ H1 + H2 + H3 + H4, rewards)
summary(model.3e)

# some means will be super handy
aggregate(trials ~ condF, rewards, mean)

# g: power
modelEffectSizes(model.3e)                   # numerator for f-squared
model.3e.R2 <- summary(model.3e)$r.squared   # denominator for f-squared

# PRE for H3
PRE.H3 <- 0.0069


# f-squared
f2.H3 <- PRE.H3 / (1 - model.3e.R2)

# power analysis
pwr.f2.test(u = 1, v = NULL, f2 = f2.H3, power = .90)$v

# h: modeling residuals
rewards <- rewards %>% 
    mutate(e = resid(model.3e))

model3h <- lm(e ~ H1 + H2 + H3 + H4, rewards)
modelEffectSizes(model3h)

# i: dropping a contrast
model3i <- lm(trials ~ H1 + H3 + H4, rewards)
summary(model3i)
summary(model.3e) # for comparison

The best (but not especially convincing, frankly) argument in favor of carrying out the ANOVA here is probably that it’s conventional to do so. Also, there are some post-test procedures that assume a significant ANOVA to occur first before carrying out the post-tests (to properly keep familywise error rate at a desirable level). But the ANOVA does not tell us anything about the hypotheses in question.

The best argument I can think of against carrying out the ANOVA is that it’s superfluous; it does not provide an answer to any of the hypotheses that are deemed worth testing.
The contrasts are orthogonal because no group is compared to another group twice. Also, the correlations among the contrast lambdas in part e is 0.
I’m going to going with the Bonferroni correction because it’s easy. I’ll set \(\alpha\) = .05 / 4 = .0125.
Five groups means we need four predictors, no more, no less. My contrast variables (i.e,. H1, H2, H3, and H4) do the trick.
Group means and SDs appear below (for ease of interpretation):

## # A tibble: 5 × 3
##   condF        M    SD
##   <fct>    <dbl> <dbl>
## 1 control1    21  3   
## 2 control2    23  2.65
## 3 hungry       7  2.55
## 4 thirsty     12  2.55
## 5 sad         11  2.55

The two control groups do not differ significantly, t(20) = 1.19, p = .25.
The three experimental groups performed significantly better than the control groups, t(20) = 11.03, p = tiny.
The two singly-deprived groups do not differ significantly from the doubly-deprived group, t(20) = 1.03, p = .32.
The hungry group performed significantly better than the thirsty group , t(20) = 2.97, p = .008.

Tukey’s HSD or a full set of pairwise t-tests would allow for more comparisons to be made (10 in this case), and thus more fine-grained conclusions might be drawn (i.e., we would know more about which groups differ), but there would be a cost to power because alpha would be chopped up into smaller pieces.
Assuming f² = 0.053, we would need 201 observations to achieve .90 power in a replication effort.
All the SSR values are now zero. This is because all of the variability between groups has been removed from the data, leaving only error. And because the error is all that’s left in the data, SSE is the same.
The slopes of H1, H3, and H4 stayed the same, but their SEs were all quite a bit higher. The reason for this is that the SE for slopes is heavily dependent on model R², which was about .87 with H3 in the model and about .07 without it. Lots more error in the model means lots more imprecision in parameter estimates, even if their values stay the same.
If one tries to use the numbers that represent the conditions as a predictor as though they have quantitative meaning, a linear relationship between that predictor and the outcome is assumed. The graph below makes it clear that this is a hilariously bad idea.

rewards %>% ggplot(aes(y = trials, x = cond)) + geom_point() + geom_smooth(se = F) + theme_bw()

BONUS: Quantifying nonlinearity

To quantify nonlinearity in the age-BMI relationship for the data in part a, we need to fit a purely-linear model (which we already did in #1b). The SSR for the linear term (i.e., age) quantifies how much linearity there is. Then, to capture every possible relationship (linear + all the nonlinear possibilities), we need to fit a model that treats age like a grouping variable; that is, we need to do an ANOVA where age - which is a quantitative thing - is treated as a grouping variable. The SSR for this model will include both linearity and nonlinearity. Subtracting the SSRs for the two models leaves us with a quantification of nonlinearity. Here we go!

modelEffectSizes(model1b) # SSR = 65.621 on 1 df

# the weird age as group model
modelAgeAsGroup <- lm(BMI ~ factor(Age), diabetes)
modelEffectSizes(modelAgeAsGroup) # SSR = 2904.16 on 50 df
anova(modelAgeAsGroup)            # MS residuals = 40.1 or so on 663 df

Now we have the two SSRs and their associated dfs. SSR for the linear-only model is 65.621, estimated using 1 df. SSR for the age-as-group model is 2904.16, estimated (expensively) using 50 df. The difference between the two is the SSR due to nonlinearity, or 2838.539, which is estimated using 50 - 1 = 49 df.

Not only can we quantify nonlinearity, but we can do a significance test of it. The SSR for nonlinearity, 2838.539, divided by its associated df of 49, equals 57.9294, which is MS for nonlinearity. This can be divided by MS for residuals for the larger (i.e., age-as-groups) model, which is provided by the anova function above as about 40.1, estimated using 663 df.

Now we can divide MS for nonlinearity by MS for residuals to get F and then convert it to a p-value.

\(F = \frac{57.93}{40.1} = 1.445\)

Despite that this is a pretty small F-ratio, it’s got lots of df associated (49 and 663, to be precise). Its p-value is 0.0279. All of this is to say that there is significant nonlinearity in the data. We already knew this, but isn’t this a fun exercise? (We know from the model in #1c that SSR for the quadratic component is 1060 or so, which means that 1060 / 2838.5 or about .37 of the nonlinearity is quadratic in nature. Fun!)

Exam 1 Answer Key

Updated to fix a typo

1

2

3

BONUS: Quantifying nonlinearity