Bonferroni-Welch/Welch t

Gastric Example

See SAS Welch t for the description of the data.

> x <- matrix(scan(),ncol=3,byrow=T)

1: 1 2 3 1 2 5 1 1 1 1 3 4 1 2 2 2 0 1 2 1 1

22: 2 1 1 2 0 0 2 2 0 3 0 0 3 0 0 3 1 1 3 1 1 3 1 0

46:

```Read 45 items
```

> gastric <- data.frame(gastric=x[,2], duodenal=x[,3],

+ drug=factor(x[,1],labels=c("N","S","C")))

> gastric

```   gastric duodenal drug
1        2        3    N
2        2        5    N
3        1        1    N
4        3        4    N
5        2        2    N
6        0        1    S
7        1        1    S
8        1        1    S
9        0        0    S
10       2        0    S
11       0        0    C
12       0        0    C
13       1        1    C
14       1        1    C
15       1        0    C
```

> car::leveneTest(gastric ~ drug, data=gastric))

```Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group  2  0.2222  0.804
12
```

> car::leveneTest(duodenal ~ drug, data=gastric))

```Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group  2  2.4615 0.1271
12
```

> oneway.test(duodenal~drug, data=gastric) # Welch F

```        One-way analysis of means (not assuming equal variances)

data:  duodenal and drug
F = 5.5597, num df = 2.000, denom df = 7.365, p-value = 0.03375
```

> oneway.test(duodenal~drug, data=gastric, var.equal=T) # F (incorrect)

```        One-way analysis of means

data:  duodenal and drug
F = 10.129, num df = 2, denom df = 12, p-value = 0.00265
```

> for(i in 1:2)for(j in (i+1):3)

+   print(t.test(duodenal~drug,data=gastric,

+                subset=(drug %in% levels(drug)[c(i,j)]),

+                conf.level=0.90)

+        )

```        Welch Two Sample t-test

data:  duodenal by drug
t = 3.2071, df = 4.946, p-value = 0.02417
alternative hypothesis: true difference in means is not equal to 0
90 percent confidence interval:
0.888457 3.911543
sample estimates:
mean in group N mean in group S
3.0             0.6

Welch Two Sample t-test

data:  duodenal by drug
t = 3.4744, df = 4.946, p-value = 0.01808
alternative hypothesis: true difference in means is not equal to 0
90 percent confidence interval:
1.088457 4.111543
sample estimates:
mean in group N mean in group C
3.0             0.4

Welch Two Sample t-test

data:  duodenal by drug
t = 0.5774, df = 8, p-value = 0.5796
alternative hypothesis: true difference in means is not equal to 0
90 percent confidence interval:
-0.4441663  0.8441663
sample estimates:
mean in group S mean in group C
0.6             0.4
```

The aove gives the individual comparisons of means. See the discussion in SAS Welch t. I've constructed several R functions to obtain Bonferroni-Welch/Welch tests/confidence intervals. Run the following R command (assuming that your computer is connected to the internet) to obtain these functions. The `welch` function assumes the existence of R package:multcomp. However, you don't have to load the package when you run this function.

> #

> # source('http://www.stat.wmich.edu/wang/R/codes/welch.R')

> #

After this, make sure that you save the image by either log-off R and answer 'yes' to exit dialog; or do FileSave Workspace and respond with .RData (case-sensitive) for the file in the default location (type in for the first time, click on .RData icon for subsequent sessions). The documentations were included as comments in the R functions.

> with(gastric, welch(duodenal,drug,linfct="Tukey",

```*******************************
* Welch Inference About Means *
*  All Pairwise Comparisons   *
*******************************

Response: duodenal

==> Linear Functions <==
N  S C
S - N -1  1 0
C - N -1  0 1
C - S  0 -1 1

==> Hypothesis Tests <==
(row names show the alternative hypotheses)

tW      DFw          P
S - N != 0 -3.2071349 4.946372 0.02417052
C - N != 0 -3.4743961 4.946372 0.01807831
C - S != 0 -0.5773503 8.000000 0.57958400

==> Confidence Intervals <==
Individual Confidence Level = 90%
Estimate       S.E   t.crit        lwr        upr
S - N     -2.4 0.7483315 2.019884 -3.9115430 -0.8884570
C - N     -2.6 0.7483315 2.019884 -4.1115430 -1.0884570
C - S     -0.2 0.3464102 1.859548 -0.8441663  0.4441663
```

The result is identical to the above t tests (only reverse signs since the contrasts in welch are of negative sign of the respective contrasts in t tests). Now, suppose Bonferroni-Welch of all pairwise comparsions are desired:

> with(gastric, welch(duodenal,drug,linfct="Tukey",

+ int.df=FALSE, conf.level=0.90))

```******************************************
* Bonferroni-Welch Inference About Means *
*        All Pairwise Comparisons        *
******************************************

Response: duodenal

==> Linear Functions <==
N  S C
S - N -1  1 0
C - N -1  0 1
C - S  0 -1 1

==> Hypothesis Tests <==
(row names show the alternative hypotheses)

tW      DFw          P
S - N != 0 -3.2071349 4.946372 0.07251157
C - N != 0 -3.4743961 4.946372 0.05423492
C - S != 0 -0.5773503 8.000000 1.00000000

==> Confidence Intervals <==
Family Confidence Level = 90%
Estimate       S.E   t.crit       lwr        upr
S - N     -2.4 0.7483315 2.922830 -4.587246 -0.2127541
C - N     -2.6 0.7483315 2.922830 -4.787246 -0.4127541
C - S     -0.2 0.3464102 2.566019 -1.088895  0.6888952
```

> with(gastric, welch(duodenal,drug,linfct="Tukey",

+ conf.level=0.90)) -> gastric.welch

> gastric.welch # use integer d.o.f.'s

```******************************************
* Bonferroni-Welch Inference About Means *
*        All Pairwise Comparisons        *
******************************************

Response: duodenal

Note: The degrees of freedom used for the calculations
of P values and critical values are integers (from
rounding down the reported degrees of freedom).

==> Linear Functions <==
N  S C
S - N -1  1 0
C - N -1  0 1
C - S  0 -1 1

==> Hypothesis Tests <==
(row names show the alternative hypotheses)

tW      DFw          P
S - N != 0 -3.2071349 4.946372 0.09803377
C - N != 0 -3.4743961 4.946372 0.07644444
C - S != 0 -0.5773503 8.000000 1.00000000

==> Confidence Intervals <==
Family Confidence Level = 90%
Estimate       S.E   t.crit       lwr         upr
S - N     -2.4 0.7483315 3.186316 -4.784420 -0.01557978
C - N     -2.6 0.7483315 3.186316 -4.984420 -0.21557978
C - S     -0.2 0.3464102 2.566019 -1.088895  0.68889520
```

Now produce plot of the confidence intervals.

> windows(width=5.5,height=3,pointsize=10)

> plot(gastric.welch,plotmath=T,

+ sub="note: integer d.o.f.'s used")

Currently, some warning messages were produced (annoying!). Will fix the problem in the next version. The plot is, however, produced.

Suppose now, we wish to compare treatments with control (Standard drug):

> levels(gastric\$drug)

```[1] "N" "S" "C"
```

Note that level 2 is the control.

> windows(width=5.5,height=2.5,pointsize=10)

> plot(with(gastric,welch(duodenal,drug,base=2,conf.level=0.9)))

The confidence intervals are

> contr <- list(psi.1=c(N=1,S=-1),

+ psi.2=c(C=1,S=-1), psi.3=c(N=1,C=1,S=-2))

> contr

```\$psi.1
N  S
1 -1

\$psi.2
C  S
1 -1

\$psi.3
N  C  S
1  1 -2
```

> x <- with(data=gastric,

+ welch(duodenal,drug,linfct=contr, conf.level=0.9,

+ divisor=c(1,1,2)))

> x

```******************************************
* Bonferroni-Welch Inference About Means *
* User-Defined Linear Functions of Means *
******************************************

Response: duodenal

Note: The degrees of freedom used for the calculations
of P values and critical values are integers (from
rounding down the reported degrees of freedom).

==> Linear Functions <==
N  S   C
psi.1 1.0 -1 0.0
psi.2 0.0 -1 1.0
psi.3 0.5 -1 0.5

==> Hypothesis Tests <==
(row names show the alternative hypotheses)

tW      DFw          P
psi.1 != 0  3.2071349 4.946372 0.09803377
psi.2 != 0 -0.5773503 8.000000 1.00000000
psi.3 != 0  2.4596748 8.226221 0.11801611

==> Confidence Intervals <==
Family Confidence Level = 90%
Estimate       S.E   t.crit         lwr       upr
psi.1      2.4 0.7483315 3.186316  0.01557978 4.7844202
psi.2     -0.2 0.3464102 2.566019 -1.08889520 0.6888952
psi.3      1.1 0.4472136 2.566019 -0.04755877 2.2475588
```

> windows(width=5.5,height=3,pointsize=10)

> plot(x,plotmath=T,left.margin=8)

The confidence intervals are plotted: