Chapter 4 R support

Base R for group comparisons and a factorial model

Use these examples beside Chapter 4 and its lecture, or adapt their short commands for an exercise, Computer Lab or homework. Computer Lab 4 selects one connected analysis for its final classroom activity round; it does not ask you to work through every section here. This page is reusable support, available before class, and completing it beforehand is not required. Choose one software route: R here, desktop R, or the lab’s SPSS instructions. The lab contains its activity and worked answers.

Every example on this site uses base R with the standard stats and graphics commands; no package is loaded or installed. Choose Run All to run a complete cell and Start Over to restore its starting code. Every cell repeats the short setup it needs, so any one of them can run on its own in a fresh session. Calculate at full precision and report statistical results rounded once to four decimal places, with trailing zeros dropped (3.1, not 3.1000); rounded inputs or different software rounding can change the last digits. For the buttons, the keyboard shortcuts and what to do when R reports an error, see how to run code.

Your task	Useful starting point
New to R, or returning after a break	Before you start
Enter or read the 30 comprehension scores of the background-sound study, or the 24 quiz scores of the feedback-by-practice study	The background-sound data or the feedback-by-practice data
Core Exercises 1-2 or homework 9	One-way ANOVA
What separation, spread and group size do to F	What moves the F statistic
Core Exercises 3-4: a declared comparison, or all pairwise differences	Contrasts and Tukey
One comparison under several families	One estimate, five families
The chance of a false rejection among several comparisons	The familywise error rate
Computer Lab 4, core Exercises 5-6 or homework 10	The full factorial model, then the interaction plot and declared comparison
Checking the fitted model against its residuals	Residuals from the factorial fit
What leaving practice out does to the error term	Practice left out
How much variation a source represents	Eta squared for each source
A common variance in doubt	Welch’s one-way ANOVA
What unequal spreads do to the two F tests	Unequal spreads with unequal group sizes
Private course data or another complete supplied script	Programming help and your workspace

Use the task’s own data, model and requested precision. The workspace and the demonstration cells have separate R objects, so copy a whole cell, including its setup, when you move code. Fixed code and output remain readable if the browser R runtime cannot start.

JavaScript is off, so the code boxes on this page cannot run. Read the fixed code and output for Chapter 4, or run the code in desktop R.

Before you start

A command is one instruction for R, such as mean(sound_data$score). An object is a named value R keeps for later, such as n. The assignment arrow <- creates one: n <- 30 stores 30 under the name n, and typing n on its own line prints it. c() combines several values into one vector, so c("Silence", "Music with lyrics") holds two labels. A table of records is an object too: sound_data <- read.csv("data/background_sound.csv") reads the 30 records of the background-sound study into a table called sound_data, and factor() with levels = then fixes the order the conditions are listed in. $ takes one named column out of a table: sound_data$score is the score column. A function is a command that does a job, written with its arguments inside the brackets. In aov(score ~ sound, data = sound_data) the formula is the first argument and data = sound_data is a named one; commas separate them, and text in quotes is a label rather than an object. Names must match exactly, so Score is not score. For the buttons, the keyboard shortcuts and what to do when R reports an error, see how to run code.

The background-sound data and a first plot

Chapter 4 works with two invented studies. The first supplies the one-way analysis, the contrasts and Welch’s one-way ANOVA. Thirty first-year university students each read the same passage for 20 minutes, in an individual session, and then answered a 20-item reading-comprehension quiz scored from 0 to 20. Each student was randomly assigned to one of three background-sound conditions, ten to each: Silence, Instrumental music, or Music with lyrics playing during the reading. The second study, 24 students in a feedback-by-practice layout, opens the full factorial model below. The data register records the identity of the file and its invented setting.

The 30 records are stored in a CSV file with the columns student, sound and score. read.csv() reads that file into a data frame. A CSV holds plain text with no order of its own, so the factor() line sets the display and model order of the three conditions explicitly, with Silence first because it is the reference for the declared comparisons. Without that line R would sort the labels alphabetically and list Instrumental music first. Tukey’s pairs involving Silence would then print as Silence minus each music condition, and the weights c(-1, 1, 0) and c(-1, 0, 1) in the contrast cell would compare different conditions from the ones their labels name. head() shows the first six rows, str() reports each column’s type and a factor’s level order, and table() counts the students in each condition.

The table has 30 rows and three columns. str() reports student as text (chr), sound as a factor with three levels in the declared order and score as whole numbers (int), and the counted table shows ten students in each condition. Every one-way cell on this page repeats these two setup lines. On desktop R, use read.csv("background_sound.csv") with the file beside your script.

The same 30 scores typed in

You do not need the file. Typing the scores keeps the whole record visible, which helps when you want to see exactly what is being analysed, or to change one score and rerun. conditions holds the three labels in the declared order. data.frame() builds the table: sprintf("S%02d", 1:30) writes the identifiers S01 to S30, rep(conditions, each = 10) repeats each condition label ten times, and the scores follow condition by condition, ten to a line. The factor() line after it, the same setup line as in the CSV cells with the labels taken from conditions, fixes the level order.

This produces the same 30 records, in the same order, as the CSV, so either starting point can feed the analyses below. The cells that follow all start from the CSV, because that is the shorter setup to repeat.

Counts, means and spreads

table() counts the students in each condition, and aggregate(score ~ sound, ...) applies mean within each condition. tapply() applies a function within each group and returns a named vector, here the three SDs.

Each condition has ten students. The means are 12, 10 and 8 comprehension points, two points apart, and the SDs, 2.7080, 2.4037 and 2.5820, are similar. A difference between means summarizes the groups; it does not say that every student in one condition outscored every student in another.

Boxplots with every observation

boxplot(score ~ sound, ...) draws one box per condition: the thick line is the median and the box spans the middle half of the scores. horizontal = TRUE lays the boxes on their side so the condition labels fit, las = 1 keeps those labels horizontal and par(mar = ...) widens the left margin for them. With the boxes on their side the score axis runs across the page, and ylim = c(0, 20) still sets it, here to the full quiz range; col and border draw pale boxes so the points stand out, and title() labels the score axis. stripchart() then adds every student: method = "jitter" moves each point up or down by a small random amount, at most jitter = .15, so equal scores stay visible, and add = TRUE draws the points over the boxes instead of starting a new plot. set.seed(12345) fixes that random movement, so every run places the points alike. This cell draws the boxes and points of the left panel of the idea-of-ANOVA figure in Chapter 4, Section 2.1, laid on their side; the group means come from the summary cell above.

The thick lines mark the medians, 12.5, 10.0 and 8.0 points, the boxes the middle half of each condition’s scores, and the points keep every student visible. The ranges, 8 to 16 for Silence, 6 to 14 for Instrumental music and 4 to 12 for Music with lyrics, overlap widely: some students who read with lyrics playing, scoring 11 or 12, outscored some who read in silence, scoring 8 or 9. A chart of the three means alone would hide this within-group variation, and the F test below weighs the separation of the means against exactly that variation.

One-way ANOVA

aov() fits the model, and its two arguments state what is compared. In the formula score ~ sound, ~ separates the outcome from its grouping factor, so the model has one mean for each sound condition (Chapter 4, Section 2.2), and data = sound_data names the 30 records of the background-sound study. summary() prints the ANOVA table: a row for the factor and a row for the within-group residuals. Its sums of squares are sequential: each source is credited with the variation it accounts for after the sources listed above it. With a single factor these are the classical between-group and within-group sums of squares of Section 2.2, whatever the group sizes; in a model with two factors the sequential and classical sums agree when every cell holds the same number of students, as in the full factorial below.

The table gives F(2, 27) = 6.0674 with p = .0067; R’s table rounds them to 6.067 and 0.00667. The null is that the three sound-condition means are equal; the alternative is that at least two differ. Under independent observations, equal group error variances and Normal errors, the null ratio has an F(2, 27) reference. This result concerns the 30 invented background-sound scores only.

Where the df and mean squares come from

The F statistic has two kinds of degrees of freedom, as Chapter 4, Section 2.3, explains: 2 between-groups df, the numerator df, and 27 within-groups df, the denominator df. deviance() returns the within-group (residual) sum of squares and df.residual() its 27 df, so their ratio is the within-group mean square MS_W. The between-group sum of squares 80 divided by its 2 df gives MS_B; the F statistic is the ratio of the two mean squares. In pf(F, df1 = 2, df2 = df_W, lower.tail = FALSE), df1 and df2 are the numerator and denominator df of the F(2, 27) reference, always in that order, and lower.tail = FALSE asks for the area to the right of the statistic: the right tail of Section 2.3, because only greater separation of the means enlarges F. qf(.95, df1 = 2, df2 = df_W) answers the reverse question: it returns the value below which .95 of F(2, 27) lies, the .05 critical value that the F-distribution figure of Section 2.3 marks. The cell first prints each condition’s variance with tapply(), so the pooling can be followed by hand.

SS_W is 178 with 27 df, so MS_W is 6.5926. This is Chapter 3’s pooled variance extended to three groups: the group variances averaged with weights equal to their df. The three variances, 7.3333, 5.7778 and 6.6667 (SDs 2.7080, 2.4037 and 2.5820), each rest on 9 df. Multiplying each by its 9 df gives 66, 52 and 60, which add to the within-group sum of squares 178, and dividing by the 27 df gives 6.5926; with equal group sizes this is the plain average of the three variances. The three group means carry 2 df because the third is fixed once the grand mean and the other two are known, and each group of ten residuals sums to zero, leaving 30 - 3 = 27. MS_B is 40, so F = 40 / 6.5926 = 6.0674 with p = .0067, as in the table. The .05 critical value is 3.3541, and the observed 6.0674 lies beyond it, so comparing F with the critical value and comparing p with .05 give the same decision. Both mean squares estimate the same error variance under the null; their ratio is not expected to be exactly 1 in an observed sample.

What moves the F statistic

Three quantities move F: the separation of the group means, the within-group spread in MS_W and the group size. Chapter 4, Section 2.3, changes one of them at a time in a table that starts from these thirty scores, and this cell rebuilds its three changed rows. They are constructed variations of the invented scores, not further studies.

ave(sound_data$score, sound_data$sound) returns each student’s own group mean, one value per student. transform() returns a copy of the table with the score column replaced. In doubled, group_mean + 2 * (score - group_mean) doubles every student’s distance from the group mean, so the three means stay at 12, 10 and 8 while the spread doubles. In halved, score - (group_mean - 10) / 2 moves each group mean halfway to the grand mean 10, to 11, 10 and 9, while every student keeps the same distance from the group mean. Each version is refitted with the formula of the one-way fit above, and anova() prints its table: the same sequential sums of squares that summary() prints, in a table whose columns can also be selected by name, as the eta-squared cell below does.

The third row doubles the group size to twenty with the same means and group SDs. No second set of sixty scores exists, so the cell computes that row from the group means and variances: with twenty students per group, SS_B_twenty is 20 times the sum of the squared distances of the group means from their average, SS_W_twenty is 19 times the sum of the group variances, because each variance now rests on 19 df, and the F reference has 2 and 60 - 3 = 57 df. eta_twenty is the share of the total sum of squares that the between-group part represents.

Doubling the spread leaves MS_B at 40 and multiplies MS_W by four, to 26.3704. Halving the separation leaves MS_W at 6.5926 and divides MS_B by four, to 10. Each change alone quarters F, from 6.0674 to F(2, 27) = 1.5169, and the p-value rises from .0067 to .2375: F reads the separation of the means relative to the spread, never either alone. Doubling the group size leaves every boxplot looking the same but doubles MS_B, so F doubles to F(2, 57) = 12.1348 with p < .0001; more students supply more information about the same pattern. R prints the last line in scientific notation because the p-value is so small: e+01 multiplies by 10 and e-05 divides by 100,000, so 1.213483e+01 is 12.1348 and 4.068856e-05 is below .0001. Eta squared barely moves, .2986 against .3101 for the thirty scores (see eta squared below), because it describes the share of variation rather than the amount of evidence.

Contrasts and multiple comparisons

A contrast gives each group mean a weight, with the weights summing to zero. The study plan declares two comparisons, each music condition against Silence: c(-1, 1, 0) is Instrumental music minus Silence and c(-1, 0, 1) is Music with lyrics minus Silence, in the declared order Silence, Instrumental music, Music with lyrics. A negative estimate means fewer items answered correctly than in silence. rbind() stacks the two weight rows; contrast_weights %*% group_means multiplies them out and drop() returns a plain named vector.

The SE of a contrast is sqrt(MS_W * sum(c^2 / n)). It extends Chapter 3’s two-group SE, sp * sqrt(1 / n1 + 1 / n0): MS_W takes the place of the two-group pooled variance, and the weights are squared because variances of independent means add. sweep(contrast_weights^2, 2, group_n, "/") divides each squared weight by its own group count, and rowSums() adds them into weight_sum, one sum for each comparison. qt(1 - .05 / (2 * 2), df = df.residual(one_way)) is the Bonferroni critical value: in its first argument the first 2 allows two tails and the second allows the two declared comparisons, giving at least 95% joint coverage under the model (Chapter 4, Section 3.3), and df names the 27 error df of MS_W, the variance estimate that every contrast SE uses. State the family before choosing the multiplier.

The estimates are -2 and -4 points. Both comparisons have SE 1.1483 with 27 error df, and the Bonferroni critical value is 2.3734, so each margin is 2.7253 points. The intervals are [-4.7253, 0.7253] for Instrumental music minus Silence and [-6.7253, -1.2747] for Music with lyrics minus Silence. Under the stated conditions the second excludes zero and the first does not; an interval containing zero does not establish equal means. The two intervals share the Silence mean, so they are positively related.

All three pairs with Tukey

TukeyHSD() answers a different question: all three pairwise differences at once. Its family is larger than the two declared comparisons, so its multiplier and intervals differ. Run it in its own cell so the two families stay visibly separate. In the call, which = "sound" names the factor whose pairs of levels form the family, conf.level = .95 sets the simultaneous level at which the three intervals cover their differences together, and ordered = FALSE, the default, keeps the declared level order, so each difference is a later level minus an earlier one.

Instrumental music minus Silence, Music with lyrics minus Silence and Music with lyrics minus Instrumental music are -2, -4 and -2, with intervals [-4.847, 0.847], [-6.847, -1.153] and [-4.847, 0.847] and adjusted p-values .2085, .0047 and .2085. Covering a third pair widens the two shared intervals from a half-width of 2.7253 to 2.8470 points. Choose a family because it matches the question you asked, never because its intervals exclude zero.

One estimate, five families

Chapter 4, Section 3.3, sets out the families of comparisons in a table and gives two critical values for one difference between two sound conditions, on 27 error df: 2.0518 for a single planned comparison and 2.5900 for Scheffe’s family of every contrast, with the declared pair, Tukey’s all pairs and Bonferroni’s three pairs in between. This cell computes all five critical values and the five intervals for one comparison, Music with lyrics minus Silence. The estimate, 8 - 12 = -4 points, and its SE, from MS_W = 178 / 27 and two groups of ten, are the same in every row; only the critical value changes with the family the interval must cover.

c() with names builds one named vector of critical values, one per family. qt(.975, df = 27) is the pointwise value for a single comparison fixed in advance. qt(1 - .05 / (2 * 2), df = 27) is the Bonferroni value for the two declared comparisons of the contrast cell, and qt(1 - .05 / (2 * 3), df = 27) the Bonferroni value for all three pairs. In qtukey(.95, nmeans = 3, df = 27), .95 is the simultaneous level, nmeans = 3 names the three means whose pairs form Tukey’s family and df = 27 the error df; qtukey() returns a quantile of the range of the means in units of one mean’s SE, and dividing by sqrt(2) converts it to units of the SE of a difference. In Scheffe’s multiplier, qf(.95, df1 = 2, df2 = 27) names the simultaneous level, the two-dimensional space of all contrasts among three means and the error df. cbind() places the critical values and the two endpoints side by side, and round(..., 4) rounds them for display.

The SE is 1.1483 points in every row. The critical values grow from 2.0518 for a single comparison to 2.3734 for the two declared comparisons, 2.4794 for Tukey’s three pairs, 2.5525 for Bonferroni’s three pairs and 2.5900 for Scheffe’s family of every contrast among the three means, and the intervals widen from [-6.3561, -1.6439] to [-6.9740, -1.0260]. The second row is the Bonferroni interval of the contrast cell and the third is Tukey’s interval for the same pair. Here every interval excludes zero, so the family changes the width of the uncertainty statement but not its direction. The family comes from the comparisons for which coverage is owed, never from which interval excludes zero. The chapter’s table also names Dunnett’s method, for comparing each group with one reference condition; its critical value needs a multivariate t quantile that base R does not supply, so it is left out here. For the declared pair it would give slightly narrower intervals than Bonferroni’s.

The familywise error rate

Section 3.3 explains why a family needs its own calibration. Each comparison tested at .05 rejects a true null hypothesis with probability .05, but the chance that at least one of several comparisons does so grows with their number. This cell computes the two lines of the chapter’s familywise figure for one to ten hypothetical comparisons whose null hypotheses are all true.

m <- 1:10 holds the numbers of comparisons. For m independent comparisons, each tested at .05, (1 - .05)^m is the chance that none of them rejects wrongly, because the chances of independent events multiply, so 1 - (1 - .05)^m is the chance of at least one false rejection: the familywise error rate. .05 / m is Bonferroni’s allocation, which tests each comparison at .05 divided by the number in the family, and 1 - (1 - .05 / m)^m is the familywise rate under that allocation. cbind() places the three columns side by side, and round(..., 4) rounds them for display.

Tested at .05 each, the familywise rate rises with the number of comparisons: .0500, .1426, .2262 and .4013 for one, three, five and ten. Under Bonferroni’s allocation it stays at or just below .05 for every m: .0492 and .0490 for three and five comparisons. Independence is what lets the chances multiply in this cell; Bonferroni’s guarantee needs no independence, because the chance that at least one of several events occurs is at most the sum of their chances, m times .05 / m, which is .05. Its price is the wider intervals of the previous cell.

The full factorial model

The second study supplies the factorial analysis. It gives 24 different Grade 7 students a fractions quiz scored from 0 to 20. Each student appears once, in one of six cells formed by three feedback conditions and two practice schedules, with four students per cell. No one supplies an observation in both schedules, so these are six separate groups of four different students. The data register records the identity of the file and its invented setting.

The feedback-by-practice data

The 24 records are stored in a CSV file with the columns student, feedback, practice and score. As with the background-sound file, the CSV has no order of its own, so the two factor() lines set each factor’s display and model order explicitly; without them R would sort alphabetically and put “Answer feedback” first. table() with two factors counts the six cells.

The table has 24 rows and four columns. str() confirms that feedback is a factor with three levels in the declared order and practice a factor with two, and the counted table shows six cells of four students. Every factorial cell on this page repeats these same three setup lines. On desktop R, use read.csv("feedback_practice.csv") with the file beside your script.

The same 24 scores typed in

Typing the scores is again a starting point that needs no file. feedback_levels and practice_levels hold the labels in the declared order. rep(feedback_levels, each = 8) repeats each feedback label eight times; schedule holds Massed four times and then Spaced four times, and rep(schedule, times = 3) repeats that block once for each feedback condition. The two factor() lines then fix the level order, as in the CSV cell.

This produces the same 24 records, in the same order, as the CSV. The factorial cells below start from the CSV.

Cell means and the full model

In the model formula, feedback * practice includes both main effects and their interaction, so the model does not assume that the schedule difference is the same in every feedback condition (Chapter 4, Section 4.4), and data = feedback_data names all 24 records. summary() again prints sequential sums of squares; because every cell holds four students they equal the classical sums of squares, and writing the two factors in the other order would not change them. tapply() with list() of the two factors returns the table of six cell means; rowMeans() and colMeans() give the equally weighted marginal means.

Six cells of four students leave 24 - 6 = 18 error df. Feedback gives F(2, 18) = 38.4, practice gives F(1, 18) = 28.8, and the interaction gives F(2, 18) = 9.6 with p = .0015. The feedback marginal means are 6, 10 and 14 and the practice marginal means are 8 and 12; these averages use equal weights across the other factor’s levels. With an interaction present, the cell means are needed to see the pattern that those averages summarize. The one-way analysis above came from the other study: 30 background-sound students in three groups, leaving 30 - 3 = 27 error df rather than 18.

Residuals from the factorial fit

A residual is one student’s score minus the fitted value for that student’s cell, so the residuals are what the six cell means leave unexplained. fitted() returns the fitted cell mean for each of the 24 students and resid() the matching residual. Plotting one against the other puts the cells side by side: plot() draws the points, abline(h = 0, lty = 2) adds a dashed horizontal line at zero, round(resid(full_fit), 4) removes tiny computing remainders from the residuals, and table(residual) counts how often each rounded value occurs.

The points fall in five vertical columns, at fitted values 6, 8, 10, 12 and 18. Two cells share the fitted mean 6, so that column holds eight students; their residuals coincide exactly, so only four dots are visible there. Every cell contributes the residuals -2, -1, 1 and 2, which is why the counted table shows each of those four values six times while the plot shows twenty dots.

Look for three things. First, the points in each column should be centred on the dashed line; they are, because fitting cell means forces the residuals within a cell to sum to zero. Second, the vertical spread should be about the same at every fitted value: a fan that widens as the fitted mean grows would question the equal error variance the F tests assume. Third, a point far from all the others would mark an unusual observation worth checking in the data. With four scores per cell this display can show only a coarse pattern, and these invented scores were built to be regular, so the plot shows what to inspect rather than evidence that the conditions hold.

Practice left out

The error row of the full model holds only the variation of students around their own cell means. Chapter 4, Section 4.4, shows what happens to that row if practice is left out of the model, so that each feedback condition becomes one group of eight students, and this cell repeats that display and prints the two tables behind its numbers. It is not an analysis the design calls for: the full model remains the analysis of this study. It shows what the second factor removes from the error term.

boxplot() draws one box per feedback condition. points() adds symbols to an existing plot at the positions you give: as.integer() turns each student’s feedback condition into its box position 1, 2 or 3, and ifelse() shifts each Massed student slightly left and each Spaced student slightly right and gives them different symbols, filled circles for Massed and triangles for Spaced. aov(score ~ feedback, data = feedback_data) states the one decision that differs from the full fit: its formula names feedback as the only factor, so every schedule difference, and every change in the feedback comparison across schedules, stays in its Residuals row. anova() prints the two tables, and tapply() the three group SDs.

Leaving practice out changes neither the feedback means, 6, 10 and 14, nor the feedback sum of squares, 256; it changes the error term. The feedback-only Residuals row holds 220 = 60 + 96 + 64 on 21 df, the within-cell error together with the practice and interaction sums of squares, so its mean square is 10.4762 against the full model’s 3.3333 on 18 df, and the feedback F falls from 38.4000 to 12.2182, now on 2 and 21 df (p = .0003). The plot shows where the extra variation comes from: each group mixes two schedules, so the Explanatory feedback group splits into two clusters, Massed scores from 8 to 12 and Spaced scores from 16 to 20, and its SD, 4.5981, is much larger than Rereading’s 1.6903, whose two schedules do not differ; Answer feedback’s SD is 2.7255. One pooled variance would misdescribe all three groups. A factor that influences the outcome belongs in the model.

Eta squared for each source

Eta squared is a descriptive summary of magnitude: one source’s sum of squares divided by the total sum of squares. It reports the share of the observed squared variation that goes with that source in this sample. anova() on a fitted model returns the table whose "Sum Sq" column holds those sums, the same sequential sums that summary() prints. They equal the classical sums in the one-way fit, which has a single factor, and in the factorial fit, because every cell holds four students. rownames() supplies the source names, so full_SS / sum(full_SS) divides each one by the total. The cell covers both studies, and each fit follows its own data setup.

In the one-way fit of the background-sound study, the between-group sum of squares is 80 and the total is 258, so eta squared is .3101. In the full factorial fit the sums of squares are 256 for feedback, 96 for practice, 64 for the interaction and 60 for the residuals, a total of 476; the four shares are .5378, .2017, .1345 and .1261, and the exact shares add to one, although the four rounded values shown add to 1.0001 because rounding raises three of them.

Read these as the fraction of variation among these observed scores, not as the fraction of comprehension or learning caused by background sound or feedback. A value depends on the study, its design and the sources the model contains, so the one-way .3101 and the factorial .5378 are two different quantities rather than one quantity measured twice. If software reports partial eta squared instead, its denominator is that term’s own sum of squares plus the residual sum of squares, not the total; say which one you are reporting.

The interaction plot and one declared comparison

interaction.plot(x.factor, trace.factor, response) calculates the cell means and joins them with one line per trace level. type = "b" draws points and lines, lty and col distinguish the three lines, trace.label names the legend and ylim fixes the score axis.

Each line joins means of different students; the lines are not repeated measurements. Nonparallel lines are descriptive evidence that the schedule comparison changes across feedback conditions. An inferential claim also needs uncertainty and a stated reference distribution.

Simple differences and a declared difference of differences

cell_means[, "Spaced"] - cell_means[, "Massed"] subtracts two columns of the cell-mean table, giving the schedule difference within each feedback condition. Subtracting two of those differences asks how the schedule comparison itself changes between two conditions. All four cell means come from different students, so under the full model their variances add and SE = sqrt(MS_E * (1 / 4 + 1 / 4 + 1 / 4 + 1 / 4)).

Chapter 4, Section 4.5, declares this comparison in the study plan before any scores are seen, as the study’s one interaction question. It is then a family of one, and qt(.975, df = error_df) supplies the ordinary t multiplier: .975 leaves 2.5% in each tail for a two-sided 95% interval, and df = error_df names the 18 error df of the full model. unname() drops an inherited label when printing; the value is unchanged. The last two lines store the interval and display it: round(interval, 4) rounds each endpoint once to four decimal places, and format(..., digits = 15, scientific = FALSE, trim = TRUE) shows the rounded endpoints without trailing zeros, both with the same number of decimals. It returns text; interval keeps its full precision.

The schedule differences are 0, 4 and 8 points for Rereading, Answer feedback and Explanatory feedback. The declared difference of differences is 4 with SE 1.8257, the multiplier is 2.1009, and the 95% interval is [0.1643, 7.8357], which excludes zero: under the model, spacing helps more under Answer feedback than under Rereading. The interaction test, p = .0015, says that some departure from additivity exists but not which one; this interval answers the one question declared in advance. Had the comparison been picked after looking at the plot, its interval would need a family covering every comparison the plot could have suggested, with a larger multiplier (Chapter 4, Section 3.3). The simple effect within Answer feedback, 12 - 8 = 4 points, is a different target: it compares two cells, so its SE is 1.291 rather than 1.8257.

Welch’s one-way ANOVA

The classical F test pools the three group variances into MS_W. Welch’s one-way ANOVA, the counterpart of Chapter 3’s Welch t test, keeps each group’s own variance instead: it weights each group mean by the precision that its own variance and size imply, and it estimates the denominator df of an approximate F reference from the group variances (Chapter 4, Section 5.2). oneway.test() runs either version. In the two calls below, the formula score ~ sound and data = sound_data name the outcome, the grouping factor and the 30 records of the background-sound study, as in aov(). var.equal = TRUE pools the three group variances into MS_W and reproduces the classical F of the one-way table; var.equal = FALSE, the function’s default, keeps each group’s own variance and gives Welch’s test. The cell first prints the three group SDs with tapply(), so the evidence about spread sits beside the two results.

The SDs are 2.7080, 2.4037 and 2.5820, the largest 1.1266 times the smallest. The classical call gives F(2, 27) = 6.0674 and p = .0067, the one-way table’s result; Welch’s call gives F(2, 17.9561) = 5.5144 and p = .0136, which R prints rounded as 17.956 and 0.01358. With similar spreads and ten students in every group the two versions lead to the same conclusion at .05, although Welch’s p is about twice as large: its F is slightly smaller, and it estimates its denominator df from the group variances, 17.9561 in place of 27. The choice between the two belongs before the analysis, made from the design and what is known about the spreads, not from whichever p-value falls below .05.

Unequal spreads with unequal group sizes

Chapter 4, Section 5.2, shows by simulation what happens to the two F tests when the group sizes differ and the spreads may differ too. This cell reruns its first design: groups of 10, 20 and 40 students with SDs 8, 4 and 4, so the smallest group varies most. The three population means are equal, so every rejection is a rejection of a true null.

set.seed(12345) fixes the random draws, and B is the number of simulated studies. The loop draws one group at a time: each column of x is that group in one simulated study, colMeans() gives its mean and apply(x, 2, var) its variance, stored one study per row of means and vars. MS_B / MS_W repeats the classical F of the one-way table. The weights w, each group’s n_j / s_j^2 as in Section 5.2, their total W, the correction term S, F_welch and nu compute Welch’s weighted statistic and its estimated denominator df. qf() is written with df1 and df2, the two df of each reference, and the last line gives the share of studies in which each test rejects at level .05. A run takes a few seconds.

R prints 0.15724 and 0.05106: the classical F rejects the true null in .1572 of the 50,000 studies, more than three times the nominal .05, and Welch’s test in .0511, the first row of the chapter’s table. To reproduce the second row, change only the line s <- c(8, 4, 4) to s <- c(4, 4, 8) and rerun: 0.01256 and 0.04990. Restore it, then change only n <- c(10, 20, 40) to n <- c(20, 20, 20) for the third row, 0.06096 and 0.05020. The classical F errs in the same directions as Chapter 3’s equal-variance t test: far too often when the smallest group varies most, too rarely when the largest group does, and only a little with equal group sizes. Welch’s test stays near .05 in all three designs. Another seed moves each rate by about .001 to .002, the simulation’s own uncertainty, while the contrast between .16 and .05 remains. This cell is a demonstration to rerun and vary, not a Computer Lab task.

Run a supplied script

Paste the whole script for an exercise, Computer Lab question or homework into your workspace, including its data-loading and factor() lines. The workspace has its own R session, so objects created in the cells above are not available there. The two public CSVs on this page are available in the workspace as data/background_sound.csv and data/feedback_practice.csv; on desktop R, keep that subfolder or change the path to your downloaded file. A script that types the scores needs no file at all.

The complete starting script below runs in order in a fresh session. For the background-sound study it reproduces the one-way ANOVA table; for the feedback-by-practice study, the factorial table, the cell means, the schedule differences and the declared interval. Copy all of it, then change only what your own task requires.

sound_data <- read.csv("data/background_sound.csv")
sound_data$sound <- factor(sound_data$sound,
  levels = c("Silence", "Instrumental music", "Music with lyrics"))
one_way <- aov(score ~ sound, data = sound_data)
summary(one_way)
feedback_data <- read.csv("data/feedback_practice.csv")
feedback_data$feedback <- factor(feedback_data$feedback,
  levels = c("Rereading", "Answer feedback", "Explanatory feedback"))
feedback_data$practice <- factor(feedback_data$practice,
  levels = c("Massed", "Spaced"))
full_fit <- aov(score ~ feedback * practice, data = feedback_data)
summary(full_fit)
cell_means <- tapply(feedback_data$score,
  list(feedback_data$feedback, feedback_data$practice), mean)
cell_means
schedule_difference <- cell_means[, "Spaced"] - cell_means[, "Massed"]
schedule_difference
declared <- schedule_difference["Answer feedback"] -
            schedule_difference["Rereading"]
error_df <- df.residual(full_fit)
MS_E <- deviance(full_fit) / error_df
declared_SE <- sqrt(MS_E * (1 / 4 + 1 / 4 + 1 / 4 + 1 / 4))
critical <- qt(.975, df = error_df)
c(estimate = unname(declared), SE = declared_SE, multiplier = critical)
interval <- unname(declared) + c(-1, 1) * critical * declared_SE
format(round(interval, 4), digits = 15, scientific = FALSE, trim = TRUE)

Commands used in the supplied scripts include data.frame() and rep() for entering scores, read.csv() and factor() for a data file and its level order, aov() and summary() for a model, anova() for a fitted model’s table, deviance() and df.residual() for the error term, TukeyHSD() for all pairs, oneway.test() for Welch’s one-way ANOVA, qt(), qtukey() and qf() for critical values, ave() and transform() for changed versions of the scores, boxplot(), stripchart(), points() and interaction.plot() for base graphs, ifelse() for choosing a value by condition, set.seed() for a reproducible jitter, and round() and format() for reporting: format(round(x, 4), digits = 15, scientific = FALSE, trim = TRUE) rounds once to four decimal places and drops trailing zeros. format() returns text; keep the numeric objects for further arithmetic. Use the family, weights and factor order declared by your task.

Open a private course CSV

Choose Open a CSV in the workspace and select the course-supplied hsb.csv. It remains in this browser session and is not uploaded. Read it with the first line of the task’s complete script:

hsb <- read.csv("hsb.csv")
names(hsb)

names() shows the column names. Keep the course’s numeric codes and run the supplied group and coding setup, weights and model. Opening a file does not by itself create the R object. If you reload the page, reopen both the CSV and your saved script; saving a script does not save the CSV records.

Function and software references

The R manuals explain data frames, factors, ANOVA models, model summaries, Tukey comparisons, studentized range quantiles, group averages, changed copies of a table, conditional values, one-way tests, t critical values, F critical values, interaction plots, boxplots, stripcharts, added points, rounding and formatting results. Every command on this page belongs to base R, stats or graphics, so desktop R needs no installation step for them.

For the Computer Lab’s SPSS route: full-factorial GLM, profile plots, one-way contrasts, cell format and recoding. Use the task’s group order, model and error term with either route. The data register describes the public records.

Your own code and the fixed output

Run a complete script, open a course CSV and save your work in the R workspace; the same page carries how to run code and the fixed code and output for this chapter, which stay readable when the browser runtime is unavailable.

runtime_watch = window.qmsbrRuntime ? window.qmsbrRuntime.watch(webROjs.webRPromise) : ""