code below demos One-Way ANOVA
library(car)## Loading required package: carDatalibrary(doBy)
library(emmeans)
library(xtable)
ds$trt<-as.factor(ds$trt)
ds## obs trt y
## 1 1 1 19.4
## 2 2 1 32.6
## 3 3 1 27.0
## 4 4 1 32.1
## 5 5 1 33.0
## 6 6 2 17.7
## 7 7 2 24.8
## 8 8 2 27.9
## 9 9 2 25.2
## 10 10 3 17.0
## 11 11 3 19.4
## 12 12 3 9.1
## 13 13 3 11.9
## 14 14 4 20.7
## 15 15 4 21.0
## 16 16 4 20.5
## 17 17 4 18.8
## 18 18 4 18.6
## 19 19 5 14.3
## 20 20 5 14.4
## 21 21 5 11.8
## 22 22 5 11.6
## 23 23 5 14.2
## 24 24 6 17.3
## 25 25 6 19.4
## 26 26 6 19.1
## 27 27 6 16.9
## 28 28 6 20.8summary(ds)## obs trt y
## Min. : 1.00 1:5 Min. : 9.10
## 1st Qu.: 7.75 2:4 1st Qu.:16.27
## Median :14.50 3:4 Median :19.25
## Mean :14.50 4:5 Mean :19.88
## 3rd Qu.:21.25 5:5 3rd Qu.:21.95
## Max. :28.00 6:5 Max. :33.00lm1 <-lm(y~trt,data=ds)
summary(lm1)##
## Call:
## lm(formula = y ~ trt, data = ds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.42 -1.51 0.74 1.50 5.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.820 1.599 18.028 1.15e-14 ***
## trt2 -4.920 2.398 -2.052 0.052284 .
## trt3 -14.470 2.398 -6.034 4.50e-06 ***
## trt4 -8.900 2.261 -3.937 0.000704 ***
## trt5 -15.560 2.261 -6.883 6.52e-07 ***
## trt6 -10.120 2.261 -4.476 0.000188 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.575 on 22 degrees of freedom
## Multiple R-squared: 0.743, Adjusted R-squared: 0.6846
## F-statistic: 12.72 on 5 and 22 DF, p-value: 6.965e-06anova(lm1)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 5 812.67 162.535 12.72 6.965e-06 ***
## Residuals 22 281.12 12.778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(lm1,pairwise~trt)## $emmeans
## trt emmean SE df lower.CL upper.CL
## 1 28.8 1.60 22 25.50 32.1
## 2 23.9 1.79 22 20.19 27.6
## 3 14.3 1.79 22 10.64 18.1
## 4 19.9 1.60 22 16.60 23.2
## 5 13.3 1.60 22 9.94 16.6
## 6 18.7 1.60 22 15.38 22.0
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 4.92 2.40 22 2.052 0.3470
## 1 - 3 14.47 2.40 22 6.034 0.0001
## 1 - 4 8.90 2.26 22 3.937 0.0081
## 1 - 5 15.56 2.26 22 6.883 <.0001
## 1 - 6 10.12 2.26 22 4.476 0.0023
## 2 - 3 9.55 2.53 22 3.778 0.0116
## 2 - 4 3.98 2.40 22 1.660 0.5705
## 2 - 5 10.64 2.40 22 4.437 0.0025
## 2 - 6 5.20 2.40 22 2.169 0.2912
## 3 - 4 -5.57 2.40 22 -2.323 0.2271
## 3 - 5 1.09 2.40 22 0.455 0.9972
## 3 - 6 -4.35 2.40 22 -1.814 0.4776
## 4 - 5 6.66 2.26 22 2.946 0.0710
## 4 - 6 1.22 2.26 22 0.540 0.9938
## 5 - 6 -5.44 2.26 22 -2.406 0.1971
##
## P value adjustment: tukey method for comparing a family of 6 estimates## use esticon to estimate things, to show how to use
## estimable linear fuction k'b
# this estimates (mu + trt 1)
lambda1 <- c(1,0,0,0,0,0)
esticon(lm1,lambda1)## beta0 Estimate Std.Error t.value DF Pr(>|t|)
## [1,] 0.0000e+00 2.8820e+01 1.5986e+00 1.8028e+01 2.2000e+01 1.1546e-14
## Lower Upper
## [1,] 2.5505e+01 32.135# this estimates (mu + trt 2)
lambda1 <- c(1,1,0,0,0,0)
esticon(lm1,lambda1)## beta0 Estimate Std.Error t.value DF Pr(>|t|)
## [1,] 0.0000e+00 2.3900e+01 1.7873e+00 1.3372e+01 2.2000e+01 4.8397e-12
## Lower Upper
## [1,] 2.0193e+01 27.607# this estimates (mu + trt 2) - (mu + trt 1)
lambda1 <- c(0,1,0,0,0,0)
esticon(lm1,lambda1)## beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower
## [1,] 0.000000 -4.920000 2.397945 -2.051757 22.000000 0.052284 -9.893034
## Upper
## [1,] 0.053