Rust inhibitors

Data from Kutner et al. (2005).

Ten experimental units assigned to each of four brands of rust inhibitors.

link <- url("https://people.stat.sc.edu/gregorkb/data/KNNLrust.txt")
rust <- read.csv(link,col.names=c("score","brand","rep"),sep = "", header = FALSE)
head(rust)

  score brand rep
1  43.9     1   1
2  39.0     1   2
3  46.7     1   3
4  43.8     1   4
5  44.2     1   5
6  47.7     1   6

Do the brands differ in effectiveness? Is there a best brand?

Rust inhibitors

Rust inhibitors

Visually compare the treatment group means with boxplots or dotplots.

boxplot(score ~ brand, data = rust)

Rust inhibitors

stripchart(score ~ brand, data = rust, vertical = T)

Rust inhibitors

stripchart(score ~ brand, data = rust, vertical = T, method = "jitter")

Randomized experiments comparing treatments

Start with \(N\) experimental units (EUs), e.g. subjects, mice, etc.

Randomly assign each EU to one of \(a\) treatment groups.

Measure on each EU after treatment a response \(Y\).

Compute the average of the responses in each treatment group

Questions we’d like to answer:

• Is the response mean the same in all treatment groups?
• If not, then which pairs of means are different?

One-way treatment effects model

Consider the model

\[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where

• \(Y_{ij}\) is the response for EU \(j\) in treatment group \(i\).
• \(\mu\) represents an overall or baseline mean.
• \(\tau_i\) is the treatment effect for treatment \(i\).
• The \(\varepsilon_{ij}\) are independent \(\mathcal{N}(0,\sigma^2)\) error terms.
• The \(n_i\) are the numbers of replicates in the treatment groups.

Of central interest are the hypotheses \[ \text{$H_0$: $\tau_1 = \dots = \tau_a$ \quad versus \quad $H_1$: $\tau_i \neq \tau_{i'}$ for some $i \neq i'$.} \] If we reject \(H_0\), we may wish to sort/compare the treatments.

Alternative “cell means” setup

An alternate version of the model is

\[ Y_{ij} = \mu_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where

• \(Y_{ij}\) is the response for EU \(j\) in treatment group \(i\).
• \(\mu_i\) represents the mean of treatment group \(i\).
• The \(\varepsilon_{ij}\) are error terms distributed as \(\mathcal{N}(0,\sigma^2)\).

In this version of the model the central hypotheses become \[ \text{$H_0$: $\mu_1 = \dots = \mu_a$ \quad versus \quad $H_1$: $\mu_i \neq \mu_i'$ for some $i \neq i'$}. \]

Goals in one-way ANOVA

In the one-way treatment effects model \[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where \(\varepsilon_{ij}\overset{\operatorname{ind}}{\sim}\mathcal{N}(0,\sigma^2)\), we wish to

1. Visualize the data.
2. Estimate the parameters \(\mu,\tau_1,\dots,\tau_a\).
3. Estimate the error term variance \(\sigma^2\).
4. Decompose the variation in the \(Y_{ij}\) as signal plus noise.
5. Test whether there is any difference in treatment group means.
6. Sort/compare the treatment means if there is any difference.
7. Check whether the model assumptions are satisfied.

Important quantities

• Treatment group means: \(\displaystyle \bar Y_{i.} = \frac{1}{n_i}\sum_{j=1}^{n_i} Y_{ij}\), \(i=1,\dots,a\).

• Treatment group standard deviations: \(\displaystyle S_i = \frac{1}{n_i - 1}\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{i.})^2\), \(i=1,\dots,a\).

• Overall mean: \(\displaystyle \bar Y_{..} = \frac{1}{N}\sum_{i=1}^a\sum_{j=1}^{n_i} Y_{ij}\), where \(N = n_1 + \dots + n_a\).

Here we use the convention of placing a “.” in the position of an index over which we have taken the sum.

Rust inhibitors

aggregate(score ~ brand, mean, data = rust) # mean of each treatment group

  brand score
1     1  43.1
2     2  89.4
3     3  68.0
4     4  40.5

aggregate(score ~ brand, sd, data = rust) # standard deviation

  brand score
1     1  3.00
2     2  2.22
3     3  2.17
4     4  2.44

Rust inhibitors

The mean of all response values was \(\bar Y_{..} = 60.25\).

Estimating the mean and treatment effects parameters

We would like estimators \(\hat \mu\) and \(\hat \tau_1,\dots,\hat \tau_a\) of \(\mu\) and \(\tau_1,\dots,\tau_a\) which satisfy \[\begin{align*} \hat \mu + \hat \tau_1& = \bar Y_{1.} \\ \hat \mu + \hat \tau_2& = \bar Y_{2.} \\ &\vdots \\ \hat \mu + \hat \tau_a& = \bar Y_{a.} \end{align*}\]

However, we cannot uniquely determine the values of \(\hat \mu\) and \(\hat \tau_1,\dots,\hat \tau_a\) from the above because there are only \(a\) equations (and \(a+1\) unknowns).

A way out is to simply set \(\hat \tau_1 = 0\) (the default strategy in R).

Estimating the mean and treatment effects parameters

Setting \(\hat \tau_1 = 0\) gives \[\begin{align*} \hat \mu & = \bar Y_{1.} \\ \hat \mu + \hat \tau_2& = \bar Y_{2.} \\ &\vdots \\ \hat \mu + \hat \tau_a& = \bar Y_{a.} \end{align*}\]

After rearranging we obtain \[\begin{align*} \hat \mu & = \bar Y_{1.} \\ \hat \tau_2& = \bar Y_{2.} - \bar Y_{1.} \\ &\vdots \\ \hat \tau_a& = \bar Y_{a.} - \bar Y_{1.} \end{align*}\]

So \(\hat \tau_2,\dots,\hat \tau_a\) represent differences in the mean from the first or “baseline” treatment group.

Rust inhibitors

For the rust inhibitor data we have \(\bar Y_{1.} = 43.14\), \(\bar Y_{2.} = 89.44\), \(\bar Y_{3.} = 67.95\), and \(\bar Y_{4.} = 40.47\), so we set \[\begin{align*} \hat \mu & = 43.14 \\ \hat \mu + \hat \tau_2& =89.44 \\ \hat \mu + \hat \tau_3& =67.95 \\ \hat \mu + \hat \tau_4& =40.47. \end{align*}\]

Rearranging gives \[\begin{align*} \hat \mu & = 43.14 \\ \hat \tau_2& = 89.44 - 43.14 = 46.3\\ \hat \tau_3& = 67.95 - 43.14 = 24.81\\ \hat \tau_4& = 40.47 - 43.14 = -2.67. \end{align*}\]

Rust inhibitors

Use lm() with as.factor() to fit the one-way ANOVA model.

# use as.factor() to designate brand as a "factor"
lm_out <- lm(score ~ as.factor(brand), data = rust)
lm_out

Call:
lm(formula = score ~ as.factor(brand), data = rust)

Coefficients:
      (Intercept)  as.factor(brand)2  as.factor(brand)3  as.factor(brand)4  
            43.14              46.30              24.81              -2.67  

Estimation of the error term variance \(\sigma^2\)

As in linear regression, define the

• fitted values \(\hat Y_{ij}\) as \(\hat Y_{ij} = \bar Y_{i.}\) for \(j = 1,\dots,n_i\), and the
• residuals \(\hat \varepsilon_{ij}\) as \(\hat \varepsilon_{ij} = Y_{ij} - \bar Y_{i.}\)

for \(j=1,\dots,n_i\), \(i = 1,\dots,a\).

Then an unbiased estimator of \(\sigma^2\) is given by \[ \hat \sigma^2 = \frac{1}{N-a} \sum_{i=1}^a (n_i - 1)S_i^2 =\frac{1}{N - a} \sum_{i=1}^a\sum_{j=1}^{n_i} \hat \varepsilon_{ij}^2. \] So \(\hat \sigma^2\) pools the group variances into a single estimate.

Divide by \(N-a\) since the \(N\) residuals depend on \(a\) estimated quantities

Rust inhibitors

Consider how to compute \(\hat \sigma^2\) from the data below:

One obtains \(\hat \sigma^2 = 6.14\).

Rust inhibitors

The value of \(\hat \sigma\) is printed in the summary() output:

summary(lm_out)

Call:
lm(formula = score ~ as.factor(brand), data = rust)

Residuals:
   Min     1Q Median     3Q    Max 
-4.270 -1.597  0.395  1.275  4.730 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         43.140      0.784   55.06   <2e-16 ***
as.factor(brand)2   46.300      1.108   41.78   <2e-16 ***
as.factor(brand)3   24.810      1.108   22.39   <2e-16 ***
as.factor(brand)4   -2.670      1.108   -2.41    0.021 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.48 on 36 degrees of freedom
Multiple R-squared:  0.986, Adjusted R-squared:  0.985 
F-statistic:  866 on 3 and 36 DF,  p-value: <2e-16

Sums of squares in the one-way ANOVA model

As in linear regression we decompose the variation in the \(Y_{ij}\) by defining:

• Total sum of squares: \(\operatorname{SS}_{\operatorname{Tot}} = \sum_{i=1}^a\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{..})^2\)
• Treatment sum of squares: \(\operatorname{SS}_{\operatorname{Trt}} = \sum_{i=1}^a n_i(\bar Y_{i.} - \bar Y_{..})^2\)
• Error sum of squares: \(\operatorname{SS}_{\operatorname{Error}} = \sum_{i=1}^a\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{i.})^2\)

We have \(\operatorname{SS}_{\operatorname{Tot}} = \operatorname{SS}_{\operatorname{Trt}} + \operatorname{SS}_{\operatorname{Error}}\).

Note that \(\operatorname{SS}_{\operatorname{Trt}}\) is computed just like \(\operatorname{SS}_{\operatorname{Reg}}\) in linear regression.

We again define \(\displaystyle R^2 = \frac{\operatorname{SS}_{\operatorname{Trt}}}{\operatorname{SS}_{\operatorname{Tot}}}\).

Sampling distributions of our sums of squares

The SS, appropriately scaled, follow chi-square distributions:

• \(\displaystyle \operatorname{SS}_{\operatorname{Tot}}/\sigma^2 \sim \chi^2_{N-1}(\phi_{\operatorname{Tot}})\)
• \(\displaystyle \operatorname{SS}_{\operatorname{Trt}}/\sigma^2 \sim \chi^2_{a-1}(\phi_{\operatorname{Trt}})\)
• \(\displaystyle \operatorname{SS}_{\operatorname{Error}}/\sigma^2 \sim \chi^2_{N-a}\),

where \(\phi_{\operatorname{Tot}}\) and \(\phi_{\operatorname{Trt}}\) are noncentrality parameters.

The mean squares in the one-way ANOVA model

Dividing \(\operatorname{SS}_{\operatorname{Trt}}\) and \(\operatorname{SS}_{\operatorname{Error}}\) by their dfs, we define:

• Treatment mean square: \(\displaystyle \operatorname{MS}_{\operatorname{Trt}}= \frac{\operatorname{SS}_{\operatorname{Trt}}}{a-1}\)
• Error mean square: \(\displaystyle \operatorname{MS}_{\operatorname{Error}} = \frac{\operatorname{SS}_{\operatorname{Error}}}{N-a}\)

The ratio \(\displaystyle \operatorname{F}_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}\) has an F distribution.

The Analysis of Variance (ANOVA) table

We often present the SS, df, and MS values in a table like this:

In the table \(\displaystyle \operatorname{F}_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}\).

The p-value is based on \(F \sim F_{a-1,N-a}\).

Rust inhibitors

Obtain the ANOVA table “manually”:

Y <- rust$score
Y.. <- mean(Y)
Yi. <- aggregate(score ~ brand, mean, data = rust)[,2]
nn <- aggregate(score ~ brand, length, data = rust)[,2]
a <- length(Yi.)
N <- sum(nn)

SStot <- sum((Y - Y..)**2)
SStrt <- sum(nn*(Yi. - Y..)**2)
SSerror <- SStot - SStrt

MStrt <- SStrt / (a - 1)
MSerror <- SSerror / (N - a)

Ftest <- MStrt / MSerror
pval <- 1 - pf(Ftest, a-1, N - a)

Rust inhibitors

Consider how to compute \(\operatorname{SS}_{\operatorname{Trt}}\) from the data. Use the fact that \(\bar Y_{..} = 60.25\).

Rust inhibitors

Obtain the ANOVA table with the anova() function on the lm() output.

anova(lm_out)

Analysis of Variance Table

Response: score
                 Df Sum Sq Mean Sq F value Pr(>F)    
as.factor(brand)  3  15953    5318     866 <2e-16 ***
Residuals        36    221       6                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing whether there is any difference in treatment means

In the one-way treatment effects model we wish to test \[ \text{$H_0$: $\tau_1 = \dots = \tau_a$ \quad versus \quad $H_1$: $\tau_i \neq \tau_{i'}$ for some $i \neq i'$.} \]

We use the overall F test of significance:

1. Compute \(\displaystyle \operatorname{F}_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}\)
2. Reject \(H_0\) at \(\alpha\) if \(\operatorname{F}_{\operatorname{test}}> F_{a-1,N-a,\alpha}\).
3. Obtain p-value as \(P(F > \operatorname{F}_{\operatorname{test}})\), where \(F \sim F_{a-1,N-a}\).

The value of \(\displaystyle \operatorname{F}_{\operatorname{test}}\) and the p-value are printed in the summary() output.

Interpretation of F statistic

Note that \(\operatorname{F}_{\operatorname{test}}\) is a ratio of the form \(\displaystyle \frac{\text{Between treatment variation}}{\text{Within treatment variation}}\).

Exercise

For which data set will the F-statistic be largest/smallest?

Some CI formulas (without familywise adjustment)

In the cell-means formulation of the model \[ Y_{ij} = \mu_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where \(\mu_i = \mu + \tau_i\), we have the following CI formulas:

Rust inhibitors

Compute 95% CIs for \(\mu_1\) and \(\mu_2 - \mu_1\).

alpha <- 0.05
me <- qt(1-alpha/2,N-a) * sqrt(MSerror) / sqrt(nn[1])
lo1 <- Yi.[1] - me
up1 <- Yi.[1] + me

A 95% CI for \(\mu_1\) is \([41.551,44.729]\).

me <- qt(1-alpha/2,N-a) * sqrt(MSerror) * sqrt(1/nn[1] + 1/nn[2])
lo21 <- Yi.[2] - Yi.[1] - me
up21 <- Yi.[2] - Yi.[1] + me

A 95% CI for \(\mu_2 - \mu_1\) is \([44.053,48.547]\).

Post-hoc comparisons of means

• If we reject \(H_0\): \(\mu_1 = \dots = \mu_a\), then we may wish to compare means.

• Call such comparisons post-hoc as we do them after the F-test.

• We may wish to compare several pairs of means, which is like testing several hypotheses at once.

• When several hypotheses are tested at once, the familywise Type I error rate is the probability that any Type I error is committed.

• We discuss two methods for post-hoc comparisons of means which control the familywise Type I error rate.

Comparing all pairs of means

• We want to build a CI for \(\mu_i - \mu_{i'}\) for all pairs \(i \neq i'\).

• Suppose the design is balanced, i.e. \(n_i = n\) for all \(i = 1,\dots,a\).

• If we build for all \(i \neq i'\) the ordinary \((1-\alpha)\times 100\%\) CIs \[ \bar Y_{i.} - \bar Y_{i'.} \pm t_{a(n-1),\alpha/2} \hat \sigma \sqrt{2/n}, \] each one will cover its target with probability \(1-\alpha\).

• But now we want simultaneous coverage with probability \(1-\alpha\), i.e. \[ P(\cap_{i \neq i'} \{\text{CI for $\mu_i - \mu_{i'}$ captures target}\}) = 1-\alpha. \]

• Above probability is called the familywise coverage.

The venerable John Tukey

John Tukey, 1915 – 2000

Multiple comparisons of means with Tukey’s HSD

• Suppose the design is balanced, i.e. \(n_i = n\) for all \(i = 1,\dots,a\).

• Suppose we could find the value \(q_{a,a(n-1),\alpha}\) such that \[ P\left(\max_{i \neq i'} \left\{ \frac{|(\bar Y_{i.} - \bar Y_{i'.}) - (\mu_i - \mu_{i'}) |}{\hat \sigma /\sqrt{n}}\right\} \leq q_{a,a(n-1),\alpha} \right) = 1-\alpha. \]

• Then with probability \(1-\alpha\) the CIs \[ \bar Y_{i.} - \bar Y_{i'.} \pm q_{a,a(n-1),\alpha} \hat \sigma / \sqrt{n} \] will simultaneously cover the targets \(\mu_i - \mu_{i'}\) for all \(i \neq i'\). Show!

• Tukey made tables of the values \(q_{a,a(n-1),\alpha}\).

• Can use the simultaneous intervals to sort/compare the means.

Tukey’s table

Table A.6 from Mohr, Wilson, and Freund (2021)

Rust inhibitors

For the rust data we have \(n = 10\) and \(a = 4\).

At \(\alpha = 0.05\) we have \(q_{a,a(n-1),\alpha} = q_{4,36,0.05} \approx 3.85\) from table.

Obtain exact value with qtukey(.95,4,36) = 3.809.

Build the Tukey HSD CI for \(\mu_2 - \mu_1\).

n <- nn[1] # n_i = n for all i
me <- qtukey(.95,a,a*(n-1)) * sqrt(MSerror) / sqrt(n)
lo21 <- Yi.[2] - Yi.[1] - me
up21 <- Yi.[2] - Yi.[1] + me

We obtain the interval \([43.316,49.284]\).

Rust inhibitors

Use TukeyHSD() on aov() output to obtain the simultaneous CIs.

# must use the aov() function instead of the lm() function
aov_out <- aov(score ~ as.factor(brand), data = rust)
TukeyHSD(aov_out)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = score ~ as.factor(brand), data = rust)

$`as.factor(brand)`
      diff    lwr     upr p adj
2-1  46.30  43.32  49.284 0.000
3-1  24.81  21.83  27.794 0.000
4-1  -2.67  -5.65   0.314 0.093
3-2 -21.49 -24.47 -18.506 0.000
4-2 -48.97 -51.95 -45.986 0.000
4-3 -27.48 -30.46 -24.496 0.000

Rust inhibitors

plot(TukeyHSD(aov_out))

Comparison of treatments with a baseline treatment

• It may be that not all pairwise comparisons are of interest.

• Then Tukey’s method is too conservative (CIs wider than necessary).

• Say we want to compare all treatments to a “baseline” treatment.

• Build CIs for \(\mu_i - \mu_1\), \(i = 2,\dots,a\), \(1\) the baseline treatment.

• This makes \(a-1\) CIs instead of \({a\choose 2}\) CIs.

• Can use Dunnett’s method, Dunnett (1964).

The equally venerable Charles Dunnett

Charles Dunnett, 1921 – 2007 (Canadian, served in WWII, photo taken in Belgium)

Dunnett’s method for comparisons with a baseline

• Assume \(n_i = n\) for all \(i\) (balanced case).

• Given a value \(d_{a,a(n-1),\alpha}\) such that \[ P\left( \max_{2 \leq i \leq a} \Big| \frac{(\bar Y_{i.} - \bar Y_{1.}) - (\mu_i - \mu_1)}{\hat \sigma\sqrt{2/n}} \Big| \leq d_{a,a(n-1),\alpha}\right) = 1 - \alpha, \] with probability \(1-\alpha\) the CIs \[ \bar Y_{i.} - \bar Y_{1.} \pm d_{a,a(n-1),\alpha} \hat \sigma\sqrt{2/n} \] will simultaneously cover the targets \(\mu_i - \mu_1\) for all \(i=2,\dots,a\).

• Dunnett made tables of the values \(d_{a,a(n-1),\alpha}\).

• Cannot sort all the means after Dunnett’s.

Dunnett’s table

Table A.5 from Mohr, Wilson, and Freund (2021)

Rust inhibitors

For the rust data we have \(n = 10\) and \(a = 4\).

At \(\alpha = 0.05\) we have \(d_{a,a(n-1),\alpha} = d_{4,36,0.05}\).

Use value \(2.44\) in the table (should be close).

Treat Brand 1 as the baseline and make comparisons with Dunnett’s.

me <- 2.44 * sqrt(MSerror) * sqrt(2/n) # margin of error for Dunnett's
lo21 <- Yi.[2] - Yi.[1] - me
up21 <- Yi.[2] - Yi.[1] + me

The \(95\%\) Dunnett’s CI for \(\mu_2 - \mu_1\) is \([43.596,49.004]\).

Follow the same procedure for the other comparisons.

Rust inhibitors

Use DunnettTest() from R package DescTools.

library(DescTools) # first time run install.packages("DescTools")
Dunnett_out <- DunnettTest(score ~ as.factor(brand), data = rust, control = "1")
Dunnett_out

  Dunnett's test for comparing several treatments with a control :  
    95% family-wise confidence level

$`1`
     diff lwr.ci  upr.ci   pval    
2-1 46.30  43.58 49.0175 <2e-16 ***
3-1 24.81  22.09 27.5275 <2e-16 ***
4-1 -2.67  -5.39  0.0475 0.0549 .  

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rust inhibitors

plot(Dunnett_out)

Dunnett’s vs Tukey’s

• Tukey’s is for comparisons between all pairs of means.
• Dunnett’s is for comparison of means with a baseline.
• So Tukey’s must make greater adjustments to control the familywise Type I error.
• Therefore Tukey intervals will be wider than Dunnett intervals.
• Tukey’s allows you to sort the means, while Dunnett’s does not.
• Both methods assume a balanced design, i.e. \(n_i = n\) for all \(i\). Modifications for unbalanced designs exist, but are not straightforward to implement in R.

Bonferroni correction

If building \(m\) CIs you can ALWAYS use the Bonferroni correction:

• Build each CI ordinarily, but use \(\alpha/m\) instead of \(\alpha\).
• Ensures simultaneous coverage of all CIs with probability \(\geq 1-\alpha\).
• True prob of simultaneous coverage may be greater than \(1-\alpha\)
• Bonferroni-corrected CIs will be wider than Dunnett’s and wider than Tukey’s if used for making those same comparisons.
• Use when we do not know how to adjust for multiple comparisons.

Rust inhibitors

Compare Brand 3 to 4 and Brand 1 to 3, using the Bonferroni correction to control the familywise error rate.

alpha <- 0.05
m <- 2
me <- qt(1 - (alpha/m)/2,a*(n-1)) * sqrt(MSerror) * sqrt(2/n)
tab <- rbind(c(Yi.[3] - Yi.[4] - me,Yi.[3] - Yi.[4] + me),
             c(Yi.[1] - Yi.[3] - me,Yi.[1] - Yi.[3] + me))
rownames(tab) <- c("3-4","1-3")
colnames(tab) <- c("lower","upper")
tab

    lower upper
3-4  24.9  30.1
1-3 -27.4 -22.2

Checking model assumptions

Validity of the foregoing analyses depends on these assumptions:

1. The responses are normally distributed around the treatment means (Check QQ plot of residuals).

2. The response has the same variance in all treatment groups (Check residuals vs fitted values plot).

3. The response values are independent of each other (No way to check; must trust experimental design).

Rust inhibitors

plot(lm_out,which = 2)

Rust inhibitors

plot(lm_out,which = 1, add.smooth = F) # we don't want the red line

Perception of slope

Do axis re-scalings affect how we perceive an x-y relationship?

For a single data set with data pairs \((X_i,Y_i)\), with \(X_i \sim \mathcal{N}(0,1)\) and \(Y_i = \mathcal{N}(X_i,1)\) for \(i = 1,\dots,50\), three scatterplot treatments were constructed:

1. “Control” used x and y plotting limits given by the range of the data.
2. “X” extended the x-limits by 1.5 in each direction.
3. “Y” extended the y-limits by 1.5 in each direction.

Each student in a class was randomly assigned a scatterplot and told to draw with a ruler the best-fitting line through the data. The slope of each student-drawn line was measured and recorded as the response.

Is the response mean the same in the three treatment groups?

Perception of slope

An artifact from each treatment group:

slope <- c(1.23,1.80,1.81,1.29,2.89,1.58,0.99,1.24,
           1.26,1.57,1.27,1.19,1.82,1.76,1.91,1.25,
           1.09,1.29,1.12,1.51,2.13,1.16,0.62,1.04)
trt <- c("X","Y","X","X","Y","X","Y","C",
         "Y","C","C","C","Y","C","X","Y",
         "X","X","Y","C","Y","X","Y","C")

Perception of slope

boxplot(slope ~ trt)

Perception of slope

lm_slope <- lm(slope ~ as.factor(trt))
summary(lm_slope)

Call:
lm(formula = slope ~ as.factor(trt))

Residuals:
   Min     1Q Median     3Q    Max 
-0.922 -0.285 -0.129  0.263  1.348 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       1.3686     0.1816    7.54  2.1e-07 ***
as.factor(trt)X   0.0514     0.2487    0.21     0.84    
as.factor(trt)Y   0.1737     0.2421    0.72     0.48    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.48 on 21 degrees of freedom
Multiple R-squared:  0.0261,    Adjusted R-squared:  -0.0666 
F-statistic: 0.282 on 2 and 21 DF,  p-value: 0.757

Perception of slope

plot(lm_slope,which = 2)

Perception of slope

plot(lm_slope,which = 1)

Levene’s test for equality of variances

Checks if the mean magnitude of the residuals is equal across groups:

1. Obtain the residuals \(\hat \varepsilon_{ij}\) from the one-way ANOVA model.
2. Treat the absolute values \(|\hat \varepsilon_{ij}|\) of the residuals as new responses.
3. Test for equal means of the new responses with the F test.

So, do the ordinary F-test with the \(|\hat \varepsilon_{ij}|\) as the responses.

Perception of slope

Perform Levene’s test:

ehat <- lm_slope$residuals
lm_levene <- lm(abs(ehat) ~ as.factor(trt))
summary(lm_levene)

Call:
lm(formula = abs(ehat) ~ as.factor(trt))

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2914 -0.1277 -0.0498  0.0822  0.7986 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)       0.2098     0.0935    2.24    0.036 *
as.factor(trt)X   0.0502     0.1281    0.39    0.699  
as.factor(trt)Y   0.3393     0.1247    2.72    0.013 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.247 on 21 degrees of freedom
Multiple R-squared:  0.303, Adjusted R-squared:  0.237 
F-statistic: 4.57 on 2 and 21 DF,  p-value: 0.0226

Perception of slope

Can also use the leveneTest() function in the R package car.

library(car)
leveneTest(slope~as.factor(trt),center = mean)

Levene's Test for Homogeneity of Variance (center = mean)
      Df F value Pr(>F)  
group  2    4.57  0.023 *
      21                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We conclude that the variances are not equal across treatment groups.