## Problem 4.1

The data have to be loaded in the workspace first. This can be done as follows:

```> WCGS <- read.table("wcgs.txt",header=TRUE,sep="\t")
> WCGS\\$arcus_f <- as.factor(WCGS\\$arcus)
> WCGS.lm1 <- lm(chol~arcus_f,data=WCGS)
> summary(WCGS.lm1)
Call:
lm(formula = chol ~ arcus_f, data = WCGS)

Residuals:
Min       1Q   Median       3Q      Max
-119.906  -28.906   -2.906   27.094  422.094

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 222.9060     0.9186 242.652  < 2e-16 ***
arcus_f1     11.5140     1.6807   6.851 8.82e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 43.11 on 3138 degrees of freedom
(14 observations deleted due to missingness)
Multiple R-Squared: 0.01474,    Adjusted R-squared: 0.01442
F-statistic: 46.93 on 1 and 3138 DF,  p-value: 8.816e-12```

So, we can see that mean effect is estimated to be 11.5: This is the difference between response measured for level 1 and level 0 of the arcus factor. In this case, the intercept is simply the baseline value, that is the mean response for level 0.

Now, if we relevel the factor such that reference level becomes level 1, we will see a different pattern of results:

```> levels(WCGS\\$arcus_f)
[1] "0" "1"
> WCGS\\$arcus_f2 <- relevel(WCGS\\$arcus_f,"1")
> levels(WCGS\\$arcus_f2)
[1] "1" "0"
> WCGS.lm2 <- lm(chol~arcus_f2,data=WCGS)
> summary(WCGS.lm2)
Call:
lm(formula = chol ~ arcus_f2, data = WCGS)

Residuals:
Min       1Q   Median       3Q      Max
-119.906  -28.906   -2.906   27.094  422.094

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  234.420      1.407 166.552  < 2e-16 ***
arcus_f20    -11.514      1.681  -6.851 8.82e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 43.11 on 3138 degrees of freedom
(14 observations deleted due to missingness)
Multiple R-squared: 0.01474,    Adjusted R-squared: 0.01442
F-statistic: 46.93 on 1 and 3138 DF,  p-value: 8.816e-12
> with(WCGS, tapply(chol,arcus_f,mean,na.rm=T))
0        1
222.9060 234.4200
> with(WCGS, tapply(chol,arcus_f2,mean,na.rm=T))
1        0
234.4200 222.9060
> table(WCGS\\$arcus_f2)

1    0
941 2211```

At present, Intercept takes the mean response for the reference value (234.4) while the (negative) slope represents the mean difference between the two conditions.

Obviously, regarding the slope paramater, this recoding doesn't affect other aspect of estimation procedure, such as standard error (indeed, they are the same in the two models), confidence interval, or associated p-value. However, we will notice some difference for the intercept, since `beta`0 takes a different value and is estimated on a different sample size (941 instead of 2211): As a consequence of reduced n, SE will be increased, and t-value will be decreased.

Finally, note that we will get strictly the same results when treating arcus as a quantitative variable (results not shown), when using 0/1 coding. However, this will not hold if we recode the 0/1 value as 1/2. Indeed, if our variable now takes only two integer values, ranging form 1 to 2, both the intercept and slope parameters will differ from what was observed in the preceding. It results from simple geometrical consideration about the estimation of slope and intercept in linear regression with continuous outcome and predictor [Fig. 1].

```> WCGS\\$arcus2 <- WCGS\\$arcus
> WCGS\\$arcus2[WCGS\\$arcus2==0] <- 2
> WCGS.lm3 <- lm(chol~arcus2,data=WCGS)
> summary(WCGS.lm3)
Call:
lm(formula = chol ~ arcus2, data = WCGS)

Residuals:
Min       1Q   Median       3Q      Max
-119.906  -28.906   -2.906   27.094  422.094

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  245.934      2.961  83.056  < 2e-16 ***
arcus2       -11.514      1.681  -6.851 8.82e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 43.11 on 3138 degrees of freedom
(14 observations deleted due to missingness)
Multiple R-Squared: 0.01474,    Adjusted R-squared: 0.01442
F-statistic: 46.93 on 1 and 3138 DF,  p-value: 8.816e-12```

First, we notice that the sign of the slope coefficient has been reversed, as before (releveling the categorical variable). Meanwhile, the value for intercept is now 245.9 (vs. 222.9 for the 0/1 coding) and its standard error is three times that obtained using 0/1 coding (result not shown). Wherefore such a difference? As can be seen from the following figure, in the case of a linear regression where slope is estimated to be unity (and we are considering centered variable for the purpose of the illustration), reversing the values taken by X “move” the relation in the opposite direction: Slope become -1 while intercept is twice the preceding (precisely, `beta`0 times (x2-x1)).

Figure 1: A simple illustration of regression slope estimation in the continuous case.

Other coding scheme could be considered, for example (-1,+1) for the two level. In this case, the difference between the two level would be largely decreased, as can be seen in the following:

```> WCGS\\$arcus3 <- WCGS\\$arcus
> WCGS\\$arcus3[WCGS\\$arcus==0] <- -1
> WCGS.lm5 <- lm(chol~arcus3,data=WCGS)
> summary(WCGS.lm5)
Call:
lm(formula = chol ~ arcus3, data = WCGS)

Residuals:
Min       1Q   Median       3Q      Max
-119.906  -28.906   -2.906   27.094  422.094

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 228.6630     0.8404 272.098  < 2e-16 ***
arcus3        5.7570     0.8404   6.851 8.82e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 43.11 on 3138 degrees of freedom
(14 observations deleted due to missingness)
Multiple R-Squared: 0.01474,    Adjusted R-squared: 0.01442
F-statistic: 46.93 on 1 and 3138 DF,  p-value: 8.816e-12```

## Problem 4.3

If we consider again the output of lm(chol~arcus,data=WCGS), we can get the fitted value with the predict or fitted functions.

```> yhat <- fitted(WCGS.lm3)
> length(yhat)
[1] 3140
> length(WCGS\\$chol[!is.na(WCGS\\$arcus)])
[1] 3152```

Note, however, that the length of vector of fitted values differ from what is expected. This comes from the fact that some observations were removed before fitting the model.

```> print(idx <- attr(WCGS.lm3\\$model,"na.action"))
97  494  619  746  910 1273 1661 1847 2113 2222 2358 2605 2613 3047
97  494  619  746  910 1273 1661 1847 2113 2222 2358 2605 2613 3047
attr(,"class")
[1] "omit"```

We can now estimate the squared correlation between the two vectors and compare it to R2.

```> cor(yhat,WCGS\\$chol[-as.numeric(idx)])^2
[1] 0.01473513
> summary(WCGS.lm3)\\$r.squared
[1] 0.01473513```

## Problem 4.4

Recall from Chapter 4 that the effect of physical activity on glucose level has been considered as significant based on overall model fit and individual comparison (linear contrasts `beta`0 `+` `beta`2 = 0 and -`beta`3 + `beta`5 = 0).

```> HERS <- read.table("hersdata.txt",header=TRUE,na.strings=".")
> HERS\$physactf <- as.factor(HERS\$physact)
> physact.desc <- c("1=Much less active","2=Somewhat less active",
+ "3=About as active","4=Somewhat more active","5=Much more active")
> HERS.lm <- lm(glucose~physactf,data=HERS,subset=diabetes==0)
> summary(HERS.lm)
Call:
lm(formula = glucose ~ physactf, data = HERS, subset = diabetes ==
0)

Residuals:
Min       1Q   Median       3Q      Max
-48.9867  -6.9867  -0.9867   5.8056  29.8571

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  98.4206     0.9393 104.784   <2e-16 ***
physactf2    -0.8584     1.0842  -0.792   0.4286
physactf3    -1.2262     1.0111  -1.213   0.2254
physactf4    -2.4339     1.0108  -2.408   0.0161 *
physactf5    -3.2777     1.1211  -2.924   0.0035 **
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

Residual standard error: 9.716 on 2027 degrees of freedom
Multiple R-squared: 0.008668,  Adjusted R-squared: 0.006712
F-statistic: 4.431 on 4 and 2027 DF,  p-value: 0.001441
> confint(HERS.lm)
2.5 %      97.5 %
(Intercept) 96.578530 100.2625914
physactf2   -2.984617   1.2677189
physactf3   -3.209060   0.7566629
physactf4   -4.416114  -0.4515950
physactf5   -5.476291  -1.0791159
> par(mfrow=c(1,2))
> plot(glucose~physactf,data=HERS)
> with(HERS[HERS\\$diabetes==0,], stripchart(glucose~physactf,vertical=T,cex=.8,pch=1))
> abline(HERS.lm1 <- lm(glucose~physact,data=HERS,subset=diabetes==0),col="red")```
Figure 2: Glucose level and Physical activity.

The overall F test might be obtained from the ANOVA table as well.

```> anova(HERS.lm)
Analysis of Variance Table

Response: glucose
Df Sum Sq Mean Sq F value   Pr(>F)
physactf     4   1673     418   4.431 0.001441 **
Residuals 2027 191345      94
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1```

The design matrix can be obtained with model.matrix; in our case, this will be a 2032 by 5 matrix and here is a truncated output:

```> head(model.matrix(HERS.lm))
(Intercept) physactf2 physactf3 physactf4 physactf5
1           1         0         0         0         1
2           1         0         0         0         0
4           1         0         0         0         0
5           1         1         0         0         0
6           1         0         1         0         0
8           1         0         0         0         1```

If we remove the intercept but include the same 5 levels factor physactf, we will obviously get diffent results than what was obtained with the preceding model.

```> HERS.lm2 <- lm(glucose~-1+physactf,data=HERS,subset=diabetes==0)
> summary(HERS.lm2)
Call:
lm(formula = glucose ~ -1 + physactf, data = HERS, subset = diabetes ==
0)

Residuals:
Min       1Q   Median       3Q      Max
-48.9867  -6.9867  -0.9867   5.8056  29.8571

Coefficients:
Estimate Std. Error t value Pr(>|t|)
physactf1  98.4206     0.9393   104.8   <2e-16 ***
physactf2  97.5621     0.5414   180.2   <2e-16 ***
physactf3  97.1944     0.3742   259.7   <2e-16 ***
physactf4  95.9867     0.3734   257.1   <2e-16 ***
physactf5  95.1429     0.6120   155.5   <2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

Residual standard error: 9.716 on 2027 degrees of freedom
Multiple R-squared:  0.99,     Adjusted R-squared:  0.99
F-statistic: 4.023e+04 on 5 and 2027 DF,  p-value: < 2.2e-16```

At present, each `beta`i represents an estimate of the mean response in the ith category, and the design matrix now looks like

```> head(model.matrix(HERS.lm2))
physactf1 physactf2 physactf3 physactf4 physactf5
1         0         0         0         0         1
2         1         0         0         0         0
4         1         0         0         0         0
5         0         1         0         0         0
6         0         0         1         0         0
8         0         0         0         0         1```

However, t-test and global F-test are not relevant anymore since they are based on a comparison to a zero baseline (check the t-value greater than 100, for example). On the contrary, we are much more interested in between-groups comparisons, taking group 1 as the reference group (“much less active” condition) since it is the most relevant for the purpose of the study (variation of glucose level with intensity of daily activities). An alternative strategy would be to use a multiple comparison approach, like the Dunnett procedure [a], to compare each group to the control one.

## Problem 4.5

Consider the following situation: We have 5 groups who were enrolled in a randomized clinical trial. Four of them have been administered different doses of a pharmaceutical preparation (whatever the effect we might be measuring). Group 1 will be the reference group (placebo), while groups 2 to 5 will be the experimental groups. Only group 5 is suspected to have a lowered mean response rate.

```> ni <- 20
> x <- gl(5,ni)
> y <- c(rnorm(ni,10,2),replicate(3,rnorm(ni,12,2)),rnorm(ni,6,2.5))
> xy <- data.frame(x=x,y=y)
> summary(xy)
x            y
1:20   Min.   : 0.6874
2:20   1st Qu.: 8.9039
3:20   Median :10.9522
4:20   Mean   :10.2007
5:20   3rd Qu.:12.3148
Max.   :17.1221
> boxplot(y~x,xlab="Group",ylab="Response")
> lines(1:5,tapply(y,x,mean),type="b",col="red",lwd=2)
> abline(h=mean(y[x==1]),col="lightgray",lwd=2)```

Results are presented in [Fig. 3].

Figure 3: Response distribution and mean profiles in the simulated drug experiment.

Now, let's look at a simple comparison between Group 1 and Group 2:

```> xy.lm <- lm(y~x,xy)
> library(contrast)
> contrast(xy.lm,list(x='1'),list(x='2'))
lm model parameter contrast

Contrast      S.E.     Lower     Upper     t df Pr(>|t|)
1 -2.051397 0.6360268 -3.314069 -0.788724 -3.23 95   0.0017```

The contrast function allows to test for any linear contrast. The help function provides the following information:

```contrast.lm             package:contrast             R Documentation

General Contrasts of Regression Coefficients

Description:

This function computes one or more contrasts of the estimated
regression coefficients in a fit from one of the functions in
Design, along with standard errors, confidence limits, t or Z
statistics, P-values.

Usage:

contrast(fit, ...)
## S3 method for class 'lm':
contrast(fit, ...)
## S3 method for class 'gls':
contrast(fit, ...)
## S3 method for class 'lme':
contrast(fit, ...)
## S3 method for class 'geese':
contrast(fit, ...)```

The complete documentation can be found on R CRAN. There is also an R vignette which describes the detailed implementation of contrasts for linear models.

“Unfortunately”, it provides Wald test and not F test as some other statistical packages, like Statistica or SPSS [b], do propose.

 About contrast coding within R There are many routes to contrast coding with R. Among the classical way, one can set up the contrast by hand, using contrast vector or matrix notation. For example, we can build a contrast like (1,1,-1,-1) aiming at comparing two conditions with two other conditions by issuing ```> C <- c(1,1,-1,-1) > aov(cbind(x11,x12,x21,x22) %*% C ~ group, data)``` (but see the excellent tutorial written by Baron and Li). As an alternative, one can specify default contrast coding, such as `> options(contrasts = c("contr.helmert", "contr.poly"))` in the R shell. This would apply to all formula-based models (lm, aov, lme, etc.) that make use of such contrasts. There is also a function that facilitate contrast coding (see ?contrasts). Furthermore, the MBESS package also provides several functions for multiple comparisons, as well as the car package. Finally, the Design package is entirely devoted to linear modelling, and the contrast package used above relies on it. To get an idea, just check the R online help: `> help.search("contrast")`

Using a standard t-test (t.test(y[x==1],y[x==2],var.equal=TRUE)) is obviously not a good idea since it doesn't consider the right residual SS.

For the F test, it happen to be very easy once we recall that such statistics is calculated as the ratio between SS for the differences and residual SS.

The same result can be obtained by checking the output of the corresponding regression coefficients since x is here considered as a categorical variable. We can see in the following listing that the t value for x2 equals 3.225, with p=0.0017.

```> summary(xy.lm)
Call:
lm(formula = y ~ x, data = xy)

Residuals:
Min      1Q  Median      3Q     Max
-5.3143 -0.9061 -0.1659  1.1268  5.2257

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   9.8450     0.4497  21.891  < 2e-16 ***
x2            2.0514     0.6360   3.225 0.001726 **
x3            2.3090     0.6360   3.630 0.000458 ***
x4            2.2583     0.6360   3.551 0.000600 ***
x5           -4.8403     0.6360  -7.610 1.98e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.011 on 95 degrees of freedom
Multiple R-Squared: 0.6608,     Adjusted R-squared: 0.6465
F-statistic: 46.27 on 4 and 95 DF,  p-value: < 2.2e-16```

## Problem 4.10

We propose to use a simulation study to compare the two approaches (regression vs. ANOVA or t-test). We thus need to define a sample of observations whose values might vary according to the level of an exogenous factor.

Let's start with a two groups experiment. We assume that yi residuals follow a normal distribution, and a model like

in the ANOVA notation, or, equivalently,

using the regression modelling notation.

```> n <- 100
> x <- gl(2,n/2,n,labels=0:1)
> y <- c(rnorm(n/2,25,3),rnorm(n/2,30,2.5))
> tapply(y,x,mean)
0        1
25.10533 29.95281
> tapply(y,x,var)
0        1
7.848420 4.981874
> plot(x,y)```

As can be seen, there is quite the same within-group variability (though we specify a lesser sd when drawing the second sample) while we observe a 5 points (approx.) difference between the groups. The data are summarized with the help of box-and-whiskers charts in [Fig. 4].

Figure 4: Boxplots of the y response for the two groups.

A simple linear regression gives

```> y.lm <- lm(y~x)
> summary(y.lm)
Call:
lm(formula = y ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-5.5379 -1.5516 -0.1123  1.5510  5.6804

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  25.1053     0.3582  70.089  < 2e-16 ***
x1            4.8475     0.5066   9.569 1.04e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.533 on 98 degrees of freedom
Multiple R-Squared: 0.483,      Adjusted R-squared: 0.4778
F-statistic: 91.57 on 1 and 98 DF,  p-value: 1.042e-15```

and the results is highly significant, if we focus on the slope estimate. Here, x is treated as a categorical variable but the same results would have been obtained if it was first converted to a numerical one (see Problem 4.1).

Now, what gives a t-test (for independent samples), assuming homoscedasticity?

```> t.test(y~x,var.equal=TRUE)
Two Sample t-test

data:  y by x
t = -9.5694, df = 98, p-value = 1.042e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-5.852736 -3.842221
sample estimates:
mean in group 0 mean in group 1
25.10533        29.95281```

We get exactly the same result, except for the sign of the t statistics. The result is highly significant, which is not very surprising given the sample size we used and the effect size we choose (+5 points between the two groups).

Now, if we extend this simulation to a three groups design, we can check that we will get exactly the same pattern of results, whether we look at the output of regression model or of an ANOVA table.

```> x <- gl(3,n/3,n,labels=0:2)
> y <- as.numeric(x) + rnorm(n)
> by(y,x,summary)
INDICES: 0
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-0.6392  0.1236  0.8728  0.8005  1.2560  2.1900

INDICES: 1
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-0.6749  1.1830  1.6820  1.6710  2.1510  3.8590

INDICES: 2
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
1.215   2.291   2.880   2.955   3.521   5.104
> plot(x,y)
> stripchart(y~x,vertical=T)
> points(1:3,tapply(y,x,mean),pch="x",cex=2,col="red")
> abline(lm(y~as.numeric(x)),col="red")```

Evolution of y between the three groups is illustrated in [Fig. 5] (boxplots issued by the plot(x,y) command have been discarded). Straight line represents the regression line, and group means are represented as red crosses.

Figure 5: Evolution of y as a function of factor level.

Regression output (see below) indicates that y values differ according to factor level as each coefficient appears highly significant.

```> summary(lm(y~x))
Call:
lm(formula = y ~ x)

Residuals:
Min       1Q   Median       3Q      Max
-2.34550 -0.63749  0.03854  0.48981  2.18849

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.8005     0.1509   5.304 7.15e-07 ***
x1            0.8701     0.2151   4.045 0.000105 ***
x2            2.1542     0.2151  10.016  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8801 on 97 degrees of freedom
Multiple R-Squared: 0.511,      Adjusted R-squared: 0.5009
F-statistic: 50.68 on 2 and 97 DF,  p-value: 8.534e-16```

In this case, each mean is compared with baseline level (the group coded 0). The first coefficient estimate, x1, is thus estimated as the difference between group 1 and group 0.

```> mean(y[x==1])-mean(y[x==0])
[1] 0.8700501```

If we use an ANOVA, we also get significant result for the factor X. But, now, the interpretation is less detailed and we only have the MS of X compared to residual MS.

```> summary(aov(y~x))
Df Sum Sq Mean Sq F value    Pr(>F)
x            2 78.522  39.261  50.684 8.534e-16 ***
Residuals   97 75.139   0.775
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1```

Assessing between-groups comparison would imply to use orthogonal contrast, for instance treatment contrast (see the previous note).

However, we will get the same output than with the lm function by using summary.lm instead of just summary on aov object:

```> summary.lm(aov(y~x))
Call:
aov(formula = y ~ x)

Residuals:
Min       1Q   Median       3Q      Max
-2.34550 -0.63749  0.03854  0.48981  2.18849

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.8005     0.1509   5.304 7.15e-07 ***
x1            0.8701     0.2151   4.045 0.000105 ***
x2            2.1542     0.2151  10.016  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8801 on 97 degrees of freedom
Multiple R-Squared: 0.511,      Adjusted R-squared: 0.5009
F-statistic: 50.68 on 2 and 97 DF,  p-value: 8.534e-16```

Indeed, without any other covariates, ANOVA and regression on categorical variable always lead to the same results.

```> model.tables(aov(y~x))
Tables of effects

x
0       1      2
-0.998 -0.1280  1.156
rep 34.000 33.0000 33.000
> plot.design(y~x)```
Figure 6: Group means difference.

## Notes

1. [a] Charles Dunnett bio.

2. [b] Markus Brauer wrote an article on the use of contrasts in ANOVA and the way to do it with SPSS. It is available through the Persee interface. There is also a tutorial on post-hoc tests (Peter Bibby). Futher, be aware that orthogonal contrasts for repeated-measures design have to be specified using SPSS script because they are not available through standard menus and dialog box.