03: Significance tests

This week, you’ll revisit the same model you fit last week, focusing now on manually recreating the significance tests reported in summary(). You’ll also practice writing up the model results using APA-style reporting.

NoteGet set up

1. Open RStudio.
2. Create a new .Rmd file for this week’s exercises.
3. Save it somewhere you can find it again.
4. Give it a clear name (for example, dapr2_lab03.Rmd).
5. In the first code chunk, load the packages you’ll need this week (and install them if you don’t have them already):
   • tidyverse

---

Research question (RQ): Is there an association between wellbeing and time spent outdoors, when controlling for the effect that social interaction has on wellbeing?

Data dictionary:

variable	description
age	Age in years of respondent
outdoor_time	Self report estimated number of hours per week spent outdoors
social_int	Self report estimated number of social interactions per week (both online and in-person)
routine	Binary 1=Yes/0=No response to the question 'Do you follow a daily routine throughout the week?'
wellbeing	Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS), a self-report measure of mental health and wellbeing. The scale is scored by summing responses to each item, with items answered on a 1 to 5 Likert scale. The minimum scale score is 14 and the maximum is 70
location	Location of primary residence (City, Suburb, Rural)
steps_k	Average weekly number of steps in thousands (as given by activity tracker if available)

NoteMore detail about this dataset

From the Edinburgh & Lothians, 100 city/suburb residences and 100 rural residences were chosen at random and contacted to participate in the study. The Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS) was used to measure mental health and wellbeing.

Participants filled out a questionnaire including items concerning: estimated average number of hours spent outdoors each week, estimated average number of social interactions each week (whether on-line or in-person), whether a daily routine is followed (yes/no). For those respondents who had an activity tracker app or smart watch, they were asked to provide their average weekly number of steps.

Recap from last week

Question 1

Read in the data from https://uoepsy.github.io/data/wellbeing_rural.csv and store it in a variable named mwdata (“mw” stands for “mental wellbeing”).

Use the function lm() to fit the linear model represented by the mathematical model formulation below, and name the result mdl.

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{outdoor\_time}) + (\beta_2 \cdot \text{social\_int}) + \epsilon \]

mwdata <- read_csv("https://uoepsy.github.io/data/wellbeing_rural.csv")
mdl    <- lm(wellbeing ~ outdoor_time + social_int, data = mwdata)

Test significance of \(\beta_1\)

Question 2

Our RQ is interested in the association between hours spent outdoors and wellbeing, holding social interaction constant. Thus the coefficient we’re interested in testing for this RQ is \(\beta_1\).

We’ll start by defining our null hypothesis H0 and our alternative hypothesis H1.

Write the null and alternative hypotheses for the significance test for \(\beta_1\), using mathematical notation. Also write the hypotheses out in one sentence each, using plain English.

🗂️ See Specifying hypotheses flash card.

\[ \begin{align} H_0 &: \beta_1 = 0 \\ H_1 &: \beta_1 \neq 0\\ \end{align} \]

Write:

$$
\begin{align}
H_0 &: \beta_1 = 0 \\
H_1 &: \beta_1 \neq 0\\
\end{align}
$$

In plain English:

• The null hypothesis is that \(\beta_1\) is equal to zero (in other words, that the slope is zero, or that the line is flat).
• The alternative hypothesis is that \(\beta_1\) is not equal to zero.

Question 3

Run the code summary(mdl) to produce the model summary.

The summary shows that for the outdoor_time coefficient, which corresponds to \(\beta_1\), the estimated coefficient is 0.19909 and the estimated standard error is 0.05060.

Use these two values to calculate the observed t-statistic for \(\beta_1\).

Does it match the t-statistic presented in the model summary?

🗂️ See Reconstructing model estimates > t value flash card.

summary(mdl)

Call:
lm(formula = wellbeing ~ outdoor_time + social_int, data = mwdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7611  -3.1308  -0.4213   3.3126  18.8406 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.62018    1.48786  19.236  < 2e-16 ***
outdoor_time  0.19909    0.05060   3.935 0.000115 ***
social_int    0.33488    0.08929   3.751 0.000232 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.065 on 197 degrees of freedom
Multiple R-squared:  0.1265,    Adjusted R-squared:  0.1176 
F-statistic: 14.26 on 2 and 197 DF,  p-value: 1.644e-06

The formula for the t-statistic for \(\beta_1\) is

\[t = \frac{\beta_1}{SE(\beta_1)} = \frac{0.19909}{0.05060}\]

Let’s use R to divide those numbers:

0.19909 / 0.05060

[1] 3.934585

The t-value for \(\beta_1\) is 3.934585, which indeed rounds to the 3.935 shown in the model summary.

Question 4

This t-statistic must be compared to a t-distribution with degrees of freedom equal to \(n - k - 1\).

• \(n\) = sample size, i.e., number of observations
• \(k\) = number of predictors, i.e., number of \(\beta\) coefficients not including the intercept

Calculate how many degrees of freedom this coefficient’s t-distribution should have.

To work out \(n\), we can look at the number of observations in the RStudio environment, or alternatively:

nrow(mwdata)

[1] 200

The model has two predictors, so \(k = 2\).

Thus:

\[ \begin{align} df &= n - k - 1 \\ &= 200 - 2 - 1 \\ &= 197 \end{align} \]

Question 5

Assume \(\alpha = .05\) (the standard in Psychology) and a two-tailed test.

Find the critical t-values (that is, the t-values which establish the boundaries for statistical significance) on the t-distribution with the degrees of freedom determined in Q4.

Does the observed t-value for \(\beta_1\) fall beyond these critical values? Can we reject the H0?

🗂️ See Reconstructing model estimates > t value flash card.

c(
  lower = round(qt(0.025, 197), 3),
  upper = round(qt(0.975, 197), 3)
)

 lower  upper 
-1.972  1.972 

If our observed t-value for \(\beta_1\) is smaller than –1.972 or greater than 1.972, then we conclude that \(\beta_1\) is significantly different from zero.

Our observed t-value is 3.935, which is greater than 1.972. Thus \(\beta_1\) is significantly different from zero, and we can reject the H0.

Variance explained

Question 6

All model summaries present two \(R^2\) values. For this model, which one should we pay attention to? Why?

🗂️ See Assessing model fit > R-squared and adjusted R-squared flash card.

The model summary presents both “Multiple R-squared” and “Adjusted R-squared”.

Multiple \(R^2\) is appropriate for linear models with one single predictor. Adjusted \(R^2\) is appropriate for linear models with more than one predictor. We should therefore pay attention to the adjusted \(R^2\).

Why is \(R^2\) appropriate for LMs with more than one predictor? Because we will always be able to explain the data better if we add in more predictors. This means that multiple \(R^2\) will always increase, the more predictors we add. To take this into account, adjusted \(R^2\) modifies this value depending on how many predictors the model contains.

Question 7

How much variance in wellbeing scores does our model explain?

🗂️ See Assessing model fit > R-squared and adjusted R-squared flash card.

summary(mdl)

Call:
lm(formula = wellbeing ~ outdoor_time + social_int, data = mwdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7611  -3.1308  -0.4213   3.3126  18.8406 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.62018    1.48786  19.236  < 2e-16 ***
outdoor_time  0.19909    0.05060   3.935 0.000115 ***
social_int    0.33488    0.08929   3.751 0.000232 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.065 on 197 degrees of freedom
Multiple R-squared:  0.1265,    Adjusted R-squared:  0.1176 
F-statistic: 14.26 on 2 and 197 DF,  p-value: 1.644e-06

Our model expalins 11.76% of the variance in wellbeing scores (based on the adjusted \(R^2\)).

Report regression models

Question 8

Imagine you’re writing a report to describe this analysis. Write one sentence that describes the structure of your model: what’s the outcome? what are the predictors? Try to also include the units/measurement scale of each variable.

This model predicts people’s WEMWBS wellbeing score as a function of how many hours they spend outdoors as well as how many social interactions they have per week.

Question 9

There are official APA guidelines on how to report regression coefficient estimates.

• For reporting an unstandardised coefficient, use this template structure (filled in with nonsense values for the sake of illustration): (\(b\) = 0.54, 95% CI [0.32, 0.76], \(p\) = .041).
  • For a standardised coefficient, replace \(b\) with \(\beta\).
• For p-values, you should report exact p-values to three decimal places (like \(p\) = .041) unless \(p\) is smaller than .001, in which case you can just write \(p\) < .001.

Based on these guidelines, write a couple sentences which

• interpret the outdoor_time coefficient,
• include APA-formatted estimates of the relevant statistical quantities, and
• address whether the coefficient is significantly different from zero and what this means for the H0.

(Tip: use confint(mdl) to get each coefficient’s 95% CI.)

🗂️ See the APA Numbers and Statistics Guide.

confint(mdl) |> round(2)

             2.5 % 97.5 %
(Intercept)  25.69  31.55
outdoor_time  0.10   0.30
social_int    0.16   0.51

According to the model, increasing one hour of outdoor time is associated with an increase of 0.20 points on the WEMWBS scale, while holding the number of social interactions constant. This estimate is significantly different from zero (\(b\) = 0.20, 95% CI [0.10, 0.30], \(p\) < .001), and we can therefore reject the H0.

(Note: the coefficient interpretation is only complete if it mentions holding the other predictor(s) constant.)

Bonus conceptual question

Question 10

You’ll notice that in the model summary, all \(\beta\) coefficients (intercept AND slopes) are associated with a t-value and a p-value. The H0 being tested for all of these parameters is that their estimates are equal to zero.

This H0 makes a lot of sense for the slope coefficients: if a slope is equal to zero, then the line is flat, meaning there’s no association between \(x\) and \(y\).

Does this H0 also make sense for the intercept? If our intercept estimate is significantly different from zero, does that tell us anything interesting?

The H0 for the intercept is that its estimate is equal to zero. In other words, that the average estimated outcome value, when all predictors are equal to zero, is zero.

This is not really a plausible or interesting H0 … for our wellbeing example, we would never really expect wellbeing scores to be equal to zero. (Especially since all the scores we actually observe are between 22 and 59.) So when the significance test for the intercept tells us that our intercept estimate is significantly different from zero, it’s like … thanks I guess? We didn’t really think it was going to be equal to zero anyway.

In short: Comparing the intercept coefficient to zero is not really that informative or interesting when we have a continuous outcome variable that we’d expect to be different from zero anyway. The intercept will almost always be significantly different from zero, and we don’t usually care.