Simple inference: t-test, ANOVA, lm

t-test, one-way ANOVA, regression with lm; introducing ggplot2

Author

Affiliation

George Michaelides

Norwich Business School, University of East Anglia

Published

July 7, 2026

Learning objectives

By the end of this tutorial you will be able to:

run one-sample, independent-samples, and paired t-tests in R and interpret the output
run a one-way ANOVA with aov() and follow up with pairwise comparisons using TukeyHSD()
fit simple and multiple regression models with lm() and read the summary() output
produce your first ggplot2 scatter plot with a fitted regression line
check regression assumptions with the base plot(fit) diagnostic panels

Setup

This tutorial uses the PS dataset we imported in Tutorial 03, plus ggplot2 for the first time. We also load psych for descriptive statistics. The setup chunk sets the knit root directory to the project folder so the data/ paths work.

library(haven)
library(ggplot2)
library(psych)

We import the PS dataset from its SPSS file, strip the labels, convert -99 to NA, and turn gender and age into factors — the same workflow you met in Tutorial 03. For the rest of this tutorial we use two derived scores as our running examples:

# Mean positive affect score per participant
ps$pa <- (ps$pa1 + ps$pa2 + ps$pa3) / 3
# Mean creativity score per participant
ps$creat <- (ps$creat1 + ps$creat2 + ps$creat3) / 3
# Mean job control score per participant
ps$con <- (ps$con1 + ps$con3 + ps$con6 + ps$con8 + ps$con9 + ps$con10) / 6

Each is the average of the items that make up the scale, which gives us a single number per participant for that construct. We will use these throughout the tutorial.

t-tests

A t-test compares the means of two groups (or one group against a reference value) and asks whether the difference is statistically significant. R’s t.test() function handles all three flavours.

One-sample t-test

The one-sample t-test compares the mean of a single sample against a hypothesised population mean. Is the mean positive affect score in our sample different from the scale midpoint of 3.5?

t.test(ps$pa, mu = 3.5)

## 
##  One Sample t-test
## 
## data:  ps$pa
## t = -0.16857, df = 331, p-value = 0.8662
## alternative hypothesis: true mean is not equal to 3.5
## 95 percent confidence interval:
##  3.372792 3.607127
## sample estimates:
## mean of x 
##   3.48996

The output gives us the observed mean (mean of x), the t-statistic, degrees of freedom, p-value, the 95% confidence interval for the mean, and a reminder of the hypothesised mu. Read it as: t(r) = …, p = …, with the 95% CI for the mean.

Independent-samples t-test

The independent-samples t-test compares the means of two unrelated groups. Do male and female participants differ in their mean positive affect? We pass a formula — pa ~ gender — which R reads as “model pa as a function of gender”.

t.test(pa ~ gender, data = ps, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  pa by gender
## t = 2.6522, df = 310.12, p-value = 0.008408
## alternative hypothesis: true difference in means between group Male and group Female is not equal to 0
## 95 percent confidence interval:
##  0.08165575 0.55105985
## sample estimates:
##   mean in group Male mean in group Female 
##             3.636704             3.320346

var.equal = FALSE requests Welch’s t-test, which does not assume equal variances in the two groups. This is the safer default and is what R uses unless you say otherwise. Setting var.equal = TRUE gives the classical Student’s t-test, which assumes equal variances — only use this if you have checked the assumption with var.test() or a Levene’s test.

Note

The formula interface y ~ group is the same syntax you will use throughout R for regression, ANOVA, and (from Tutorial 05 onward) lavaan models. Get comfortable with the tilde ~ now — it means “is modelled by” or “depends on”.

Paired t-test

The paired t-test compares two measurements from the same participants — for example, pre- and post-intervention scores. The PS dataset is cross-sectional, so we do not have a natural paired example, but the syntax is:

t.test(pre_scores, post_scores, paired = TRUE)

The paired = TRUE argument tells R to compute the difference score for each participant and test whether the mean difference is zero.

One-way ANOVA

Analysis of variance (ANOVA) generalises the t-test to three or more groups. We test whether mean positive affect differs across the five age bands in the PS dataset.

Fitting the model

fit_aov <- aov(pa ~ age, data = ps)
summary(fit_aov)

##              Df Sum Sq Mean Sq F value Pr(>F)
## age           4    8.1   2.017   1.735  0.142
## Residuals   325  377.8   1.162               
## 13 observations deleted due to missingness

The summary() output shows two rows: age (the between-groups effect) and Residuals (the within-groups variation). The F-statistic tests whether the group means are all equal. The p-value tells you whether to reject that null hypothesis.

Pairwise comparisons

A significant ANOVA tells you that some group differences exist, but not which groups differ. Follow up with Tukey’s Honest Significant Difference test, which performs all pairwise comparisons while controlling the family-wise error rate:

TukeyHSD(fit_aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pa ~ age, data = ps)
## 
## $age
##            diff        lwr        upr     p adj
## 2-1 -0.45873016 -1.1945650 0.27710467 0.4291283
## 3-1 -0.50836551 -1.2121410 0.19540999 0.2773531
## 4-1 -0.62160778 -1.3361364 0.09292082 0.1216173
## 5-1 -0.70158730 -1.5179228 0.11474815 0.1298360
## 3-2 -0.04963535 -0.5010203 0.40174962 0.9981919
## 4-2 -0.16287762 -0.6308515 0.30509626 0.8749829
## 5-2 -0.24285714 -0.8551087 0.36939445 0.8126295
## 4-3 -0.11324227 -0.5289905 0.30250598 0.9451098
## 5-3 -0.19322179 -0.7665446 0.38010103 0.8873041
## 5-4 -0.07997952 -0.6664522 0.50649318 0.9958196

The output lists every pairwise comparison, the difference in means, the lower and upper bounds of the 95% confidence interval, and the adjusted p-value (p adj). A comparison is significant at the .05 level if p adj is below .05.

ANOVA requires a factor

If your grouping variable is numeric (e.g. age coded 1, 2, 3, 4, 5), aov() will treat it as a continuous predictor and run a regression instead of an ANOVA. Always confirm that the grouping variable is a factor with class(ps$age) before fitting the model. We converted age to a factor in the setup chunk above for exactly this reason.

Two-way ANOVA

To test the effects of two grouping variables and their interaction, list both on the right-hand side of the formula separated by *:

fit_2way <- aov(pa ~ gender * age, data = ps)
summary(fit_2way)

The * expands to gender + age + gender:age, where the colon denotes the interaction term. We will not pursue two-way ANOVA further here; the principle is the same as one-way, with one extra row in the summary table per additional effect.

Simple linear regression with `lm()`

Regression models a continuous outcome as a linear function of one or more predictors. The workhorse function is lm() — “linear model”.

Fitting the model

Does positive affect predict creativity? We fit the regression of creat on pa:

fit_lm <- lm(creat ~ pa, data = ps)
summary(fit_lm)

## 
## Call:
## lm(formula = creat ~ pa, data = ps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01806 -0.60540  0.00231  0.57600  2.00231 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.46712    0.15462   15.96  < 2e-16 ***
## pa           0.17685    0.04231    4.18 3.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8355 on 330 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.05028,    Adjusted R-squared:  0.0474 
## F-statistic: 17.47 on 1 and 330 DF,  p-value: 3.743e-05

The summary() output for an lm object has four sections worth understanding:

Call — the model formula, repeated for reference.
Residuals — a five-number summary of the differences between observed and predicted values. Symmetric residuals (median near zero) are a good sign.
Coefficients — a table with one row per parameter. The (Intercept) row is the predicted creat when pa is zero; the pa row is the slope — the change in creat for a one-unit increase in pa. The columns are the estimate, its standard error, the t-statistic, and the two-tailed p-value.
Model fit — Multiple R-squared is the proportion of variance in creat explained by pa. Adjusted R-squared corrects for the number of predictors. The F-statistic tests the overall model — whether the slope is significantly different from zero.

Confidence intervals for coefficients

confint(fit_lm)

##                  2.5 %    97.5 %
## (Intercept) 2.16294819 2.7712944
## pa          0.09361827 0.2600913

This gives the 95% confidence interval for each coefficient. If the interval for the slope does not include zero, the effect is statistically significant at the .05 level — equivalent to the p-value in the summary() output being below .05.

Plotting the regression with `ggplot2`

This is our first ggplot2 plot. The grammar takes getting used to, but the payoff is publication-quality figures with consistent styling.

ggplot(ps, aes(x = pa, y = creat)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, colour = "blue") +
  labs(x = "Positive affect (mean of 3 items)",
       y = "Creativity (mean of 3 items)",
       title = "Creativity as a function of positive affect") +
  theme_minimal()

Reading the code:

ggplot(ps, aes(x = pa, y = creat)) — initialise the plot with the data frame and the aesthetic mapping (which variable goes on which axis).
geom_point(alpha = 0.5) — add a scatter of points; alpha = 0.5 makes them semi-transparent so overlapping points are visible.
geom_smooth(method = "lm", se = TRUE) — add a fitted regression line; se = TRUE shades the 95% confidence band around the line.
labs(...) — label the axes and title.
theme_minimal() — a clean white background with subtle gridlines.

The ggplot2 mental model

A ggplot2 figure is built up in layers connected by +. You always start with ggplot(data, aes(...)), then add geoms (point, line, smooth, bar, boxplot), then scales, labels, and a theme. Every ggplot2 plot in this series follows this same layer-by-layer pattern. We will use ggplot2 for all plots from this tutorial onward.

Residual diagnostics

The base plot(fit_lm) call produces four diagnostic plots. We display them in a 2 × 2 grid with par(mfrow = c(2, 2)) — base R graphics, which we use here because these diagnostic plots come straight from lm itself.

par(mfrow = c(2, 2))
plot(fit_lm)

par(mfrow = c(1, 1))

The four panels are:

Residuals vs Fitted — checks linearity. A roughly horizontal red line with no clear pattern is what you want. A curved line suggests a non-linear relationship.
Normal Q–Q — checks normality of residuals. Points should lie close to the diagonal dashed line. Systematic departure from the line at the tails indicates non-normal residuals.
Scale–Location — checks homoscedasticity (equal variance). A roughly horizontal line with evenly spread points is what you want. A funnel shape indicates heteroscedasticity.
Residuals vs Leverage — flags influential observations. Points with high leverage and large residuals (top-right corner) are worth inspecting; they may disproportionately influence the fitted line.

Try it

Fit a simple regression of creat on con (job control as predictor of creativity). Is the slope significant? What is the R²? Produce a ggplot2 scatter with the fitted line.

Solution

fit_con <- lm(creat ~ con, data = ps)
summary(fit_con)

## 
## Call:
## lm(formula = creat ~ con, data = ps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.06211 -0.70988 -0.07051  0.61085  1.97988 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.26370    0.25919  12.592   <2e-16 ***
## con         -0.05040    0.06778  -0.744    0.458    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8562 on 340 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.001623,   Adjusted R-squared:  -0.001313 
## F-statistic: 0.5529 on 1 and 340 DF,  p-value: 0.4577

ggplot(ps, aes(x = con, y = creat)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", colour = "blue") +
  labs(x = "Job control (mean of 6 items)",
       y = "Creativity (mean of 3 items)") +
  theme_minimal()

The slope of con and its p-value appear in the coefficients table of the summary() output. The multiple R-squared is reported further down in the same summary.

Multiple regression

Multiple regression extends lm() to two or more predictors. The formula uses + to add predictors:

fit_multi <- lm(creat ~ pa + con, data = ps)
summary(fit_multi)

## 
## Call:
## lm(formula = creat ~ pa + con, data = ps)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9628 -0.6282  0.0228  0.5525  1.8611 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.78140    0.27995   9.935  < 2e-16 ***
## pa           0.18514    0.04271   4.335 1.94e-05 ***
## con         -0.09132    0.06785  -1.346    0.179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8344 on 329 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.05548,    Adjusted R-squared:  0.04974 
## F-statistic: 9.662 on 2 and 329 DF,  p-value: 8.363e-05

Each coefficient is now interpreted as the effect of that predictor holding the other predictors constant. The slope for pa is the change in creat for a one-unit increase in pa, with con held constant — and vice versa.

Compare the slope of pa in this multiple regression to its slope in the simple regression above. It will usually be different (often smaller) because some of the variance pa shared with con has been statistically controlled.

Confidence intervals

confint(fit_multi)

##                  2.5 %    97.5 %
## (Intercept)  2.2306774 3.3321313
## pa           0.1011256 0.2691498
## con         -0.2247964 0.0421528

Comparing nested models

To test whether adding a predictor significantly improves the model fit, use anova() to compare the nested models:

anova(fit_lm, fit_multi)

## Analysis of Variance Table
## 
## Model 1: creat ~ pa
## Model 2: creat ~ pa + con
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    330 230.33                           
## 2    329 229.07  1    1.2613 1.8115 0.1792

The output tests the null hypothesis that the additional predictor’s coefficient is zero. A significant p-value means adding con to a model that already contains pa significantly improves the fit.

Note

This anova() call is different from the aov() we used earlier. Here, anova() takes two (or more) fitted lm objects and performs a likelihood-ratio test between them. With aov(), you fit a single ANOVA model and summary() partitions the variance.

From regression to path analysis

The multiple regression model creat ~ pa + con is, in fact, a tiny path model — a system of regression equations with multiple predictors of a single outcome. In Tutorial 09 we will extend this idea to a model with multiple outcomes and multiple equations, fitted with lavaan::sem() instead of lm(). The jump from lm() to sem() is smaller than it looks: the model syntax changes from a formula to a string, but the underlying idea — modelling an outcome as a linear function of predictors — is identical.

Similarly, the regression of creat on pa is the simplest possible building block of the structural models we will fit in Tutorials 10 and 12. Keep the structure of summary(fit_lm) in mind: the coefficient table, the R², and the F-test all reappear in lavaan output, just in slightly different formatting.

Common errors

not all arguments have the same length

You passed vectors of different lengths to t.test() — usually because one vector has NAs that the other does not, or because you subsetted one but not the other. Use the formula interface (y ~ group, data = df) instead of passing separate vectors; R then handles missing data consistently.

`contrasts can be applied only to factors with 2 or more levels`

You tried to run an ANOVA with a grouping variable that has only one level — perhaps because you accidentally filtered the data down to a single group first. Check table(ps$age) before fitting the model.

`NA` coefficients in `lm()` output

A coefficient shows as NA in the summary() table when the predictor is a linear combination of other predictors (perfect multicollinearity). This happens most often when you include a categorical predictor as both a numeric variable and a factor in the same model.

Exercises

Run an independent-samples t-test comparing male and female participants on mean creativity (creat). Report the t-statistic, degrees of freedom, p-value, and 95% CI in APA format. Use var.equal = FALSE (Welch’s test).
Run a one-way ANOVA testing whether mean job control (con) differs across the five age bands. Is the overall F-test significant? If so, run TukeyHSD() and report which pairwise comparisons are significant at the .05 level.
Fit a multiple regression model predicting creat from pa, con, and age (as a factor). Interpret each coefficient. Use confint() to get the 95% CIs. Which predictors are statistically significant?
Use ggplot2 to produce a scatter plot of con (x-axis) against creat (y-axis), with points coloured by gender. Add a separate fitted regression line for each gender using geom_smooth(method = "lm", aes(colour = gender)).

Solutions

Exercise 1

t.test(creat ~ gender, data = ps, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  creat by gender
## t = 0.29736, df = 320.07, p-value = 0.7664
## alternative hypothesis: true difference in means between group Male and group Female is not equal to 0
## 95 percent confidence interval:
##  -0.1578933  0.2141217
## sample estimates:
##   mean in group Male mean in group Female 
##             3.097378             3.069264

Read the t-statistic, df, p-value, and 95% CI from the output. In APA format: t(df) = value, p = value, 95% CI [lower, upper].

Exercise 2

fit_age <- aov(con ~ age, data = ps)
summary(fit_age)

##              Df Sum Sq Mean Sq F value Pr(>F)
## age           4   2.67  0.6679   1.443  0.219
## Residuals   325 150.39  0.4628               
## 13 observations deleted due to missingness

TukeyHSD(fit_age)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = con ~ age, data = ps)
## 
## $age
##           diff        lwr       upr     p adj
## 2-1 0.06190476 -0.4023827 0.5261922 0.9961691
## 3-1 0.12998713 -0.3140720 0.5740462 0.9295799
## 4-1 0.24526370 -0.2055802 0.6961076 0.5682286
## 5-1 0.30000000 -0.2150807 0.8150807 0.4999212
## 3-2 0.06808237 -0.2167266 0.3528914 0.9654838
## 4-2 0.18335893 -0.1119171 0.4786350 0.4332485
## 5-2 0.23809524 -0.1482153 0.6244058 0.4411274
## 4-3 0.11527657 -0.1470468 0.3776000 0.7482399
## 5-3 0.17001287 -0.1917349 0.5317606 0.6980693
## 5-4 0.05473630 -0.3153086 0.4247812 0.9942803

The summary() output gives the F-test for the overall effect of age on job control. If Pr(>F) is below .05, the omnibus test is significant. TukeyHSD() then lists all pairwise comparisons; significant differences are those with p adj below .05.

Exercise 3

fit_full <- lm(creat ~ pa + con + age, data = ps)
summary(fit_full)

## 
## Call:
## lm(formula = creat ~ pa + con + age, data = ps)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0302 -0.5468  0.0386  0.5961  1.7806 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.00638    0.32060   9.377  < 2e-16 ***
## pa           0.17090    0.04207   4.062 6.12e-05 ***
## con         -0.06463    0.06668  -0.969 0.333175    
## age2        -0.20649    0.20139  -1.025 0.305974    
## age3        -0.10678    0.19318  -0.553 0.580806    
## age4        -0.42067    0.19736  -2.132 0.033802 *  
## age5        -0.79722    0.22557  -3.534 0.000469 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8053 on 323 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1296, Adjusted R-squared:  0.1134 
## F-statistic: 8.012 on 6 and 323 DF,  p-value: 4.49e-08

confint(fit_full)

##                   2.5 %      97.5 %
## (Intercept)  2.37565348  3.63710502
## pa           0.08811995  0.25367033
## con         -0.19581637  0.06655931
## age2        -0.60269825  0.18971160
## age3        -0.48682286  0.27326059
## age4        -0.80894708 -0.03240084
## age5        -1.24099429 -0.35344875

Each coefficient’s significance is judged by its p-value in the summary() output (or, equivalently, whether its 95% CI in confint() excludes zero). The pa and con slopes are interpreted as the change in creat for a one-unit increase in that predictor, holding the others constant. The age coefficients are differences relative to the reference level (the first level of the factor, here age band 1).

Exercise 4

ggplot(ps, aes(x = con, y = creat, colour = gender)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Job control (mean of 6 items)",
       y = "Creativity (mean of 3 items)",
       colour = "Gender") +
  theme_minimal()

The colour = gender aesthetic inside ggplot() maps gender to point colour; passing the same aesthetic to geom_smooth() makes it fit a separate line for each gender.

Summary & next tutorial

You can now run t-tests, one-way ANOVA, and regression with lm(), read the summary() output, and produce your first ggplot2 figures. The regression framework you have met here — a formula y ~ x, a fitted object, a summary() call — is the direct ancestor of the lavaan models we will fit from Tutorial 05 onward. In the next tutorial we make the jump: we introduce lavaan’s model syntax (=~, ~, ~~, :=) and the three functions (cfa(), sem(), lavaan()) that you will use for the rest of this series.

References

Dalgaard, P. (2008). Introductory statistics with R (2nd ed.). Springer.

R Core Team. (2024). R: A language and environment for statistical computing. https://www.R-project.org/

Wickham, H. (2019). Advanced R (2nd ed.). CRC Press.

Learning objectives

Setup

t-tests

One-sample t-test

Independent-samples t-test

Paired t-test

One-way ANOVA

Fitting the model

Pairwise comparisons

Two-way ANOVA

Simple linear regression with lm()

Fitting the model

Confidence intervals for coefficients

Plotting the regression with ggplot2

Residual diagnostics

Multiple regression

Confidence intervals

Comparing nested models

From regression to path analysis

Common errors

contrasts can be applied only to factors with 2 or more levels

NA coefficients in lm() output

Exercises

Summary & next tutorial

Further reading

References

Simple linear regression with `lm()`

Plotting the regression with `ggplot2`

`contrasts can be applied only to factors with 2 or more levels`

`NA` coefficients in `lm()` output