fruitfly <- read.table("C:/Rexamples/fruitflydata.txt", header=T) View(fruitfly) summary(fruitfly) fruitfly$line2 <- fruitfly$line fruitfly$line2 <- sub("SS", "S", fruitfly$line2, ignore.case =FALSE, fixed=FALSE) fruitfly$line2 <- sub("RS", "S", fruitfly$line2, ignore.case =FALSE, fixed=FALSE) fruitfly attach(fruitfly) # get summary info and check assumptions by(fecund, line, function(x) summary(x)) by(fecund, line, function(x) sd(x)) boxplot(fecund~line) by(fecund, line, function(x) stem(x)) qqnorm(fecund[line=='NS'], main="Normal Q-Q Plot, NS line") qqline(fecund[line=='NS']) shapiro.test(fecund[line=='NS']) qqnorm(fecund[line=='RS'], main="Normal Q-Q Plot, RS line") qqline(fecund[line=='RS']) shapiro.test(fecund[line=='RS']) qqnorm(fecund[line=='SS'], main="Normal Q-Q Plot, SS line") qqline(fecund[line=='SS']) shapiro.test(fecund[line=='SS']) # model comparison F tests anova_model1 <- aov(fecund~line) summary(anova_model1) anova_model2 <- aov(fecund~line2) summary(anova_model2) anova(anova_model1,anova_model2) # additional commands to estimate contrasts # create indicator variables for NS, RS, NS, S fruitfly$xNS <- as.double(fruitfly$line == "NS") fruitfly$xRS <- as.double(fruitfly$line == "RS") fruitfly$xSS <- as.double(fruitfly$line == "SS") fruitfly$xS <- as.double(fruitfly$line2 == "S") fruitfly$x <- fruitfly$xRS + fruitfly$xSS attach(fruitfly) # model1 cell means model 3 means: NS, RS, SS (-1 requests no intercept) fit_model1 <- lm(fecund ~ xNS + xRS + xSS -1) predict(fit_model1, se.fit=TRUE, interval="confidence") # estimate the contrast RS vs SS out_model1 <- data.frame(xNS=c(0), xRS=c(1), xSS=c(-1)) estimate_RSvsSS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_RSvsSS # model2 cell means model 2 means: NS and S (-1 requests no intercept) fit_model2 <- lm(fecund ~ xNS + xS -1) predict(fit_model2, se.fit=TRUE, interval="confidence") # estimate the contrast NS vs S out_model2 <- data.frame(xNS=c(1), xS=c(-1)) estimate_NSvsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") estimate_NSvsS # --------------------------------------------------------------------------------------- # alternate approach: model1 cell means model 3 means: NS, RS, SS (-1 requests no intercept) fit_model1 <- lm(fecund ~ xNS + xRS + xSS -1) predict(fit_model1, se.fit=TRUE, interval="confidence") # estimate the contrast RS vs SS out_model1 <- data.frame(xNS=c(0), xRS=c(1), xSS=c(-1)) estimate_RSvsSS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_RSvsSS # estimate the contrast RS vs NS out_model1 <- data.frame(xNS=c(-1), xRS=c(1), xSS=c(0)) estimate_RSvsNS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_RSvsNS # estimate the contrast SS vs NS out_model1 <- data.frame(xNS=c(-1), xRS=c(0), xSS=c(1)) estimate_SSvsNS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_SSvsNS # estimate the contrast NS vs others out_model1 <- data.frame(xNS=c(1), xRS=c(-.5), xSS=c(-.5)) estimate_NSvsothers <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_NSvsothers # ------------------------------------------------- # Scheffe simultaneous confidence intervals # ------------------------------------------------- estimate <- as.data.frame(rbind(estimate_RSvsSS$fit,estimate_RSvsNS$fit,estimate_SSvsNS$fit,estimate_NSvsothers$fit)) stderr <- as.data.frame(rbind(estimate_RSvsSS$se.fit,estimate_RSvsNS$se.fit,estimate_SSvsNS$se.fit,estimate_NSvsothers$se.fit)) multi <- sqrt(2*qf(.95,2,72)) lowerCL <- estimate$fit - multi*stderr upperCL <- estimate$fit + multi*stderr scheffe_intervals <- cbind(estimate$fit,stderr,lowerCL,upperCL) colnames(scheffe_intervals) <- c("estimate", "stderr", "lowerCL", "upperCL") rownames(scheffe_intervals) <- c("RSvsSS", "RSvsNS", "SSvsNS", "NSvsothers") scheffe_intervals --------------------- NOTE ------------------------ The first line above which reads the data from the file "fruitflydata.txt" assumes that this file is in the directory "C:/Rexamples". Obey this convention or modify the command. -------------------------- results ---------------------------------- The boxplots and normal probability plots (qqplots) requested in this program will appear in the plot window of RStudio. They can be saved as pdfs or images. I have not posted them here. These plots are somewhat crude since I did not specify titles, labels, or such. The remainder of this file is copied and pasted from the RStudio command window. -------------------------------------------------------- > fruitfly <- read.table("C:/Rexamples/fruitflydata.txt", header=T) > View(fruitfly) > summary(fruitfly) line fecund NS:25 Min. :10.80 RS:25 1st Qu.:20.30 SS:25 Median :26.70 Mean :27.42 3rd Qu.:34.50 Max. :51.80 > fruitfly$line2 <- fruitfly$line > fruitfly$line2 <- sub("SS", "S", fruitfly$line2, ignore.case =FALSE, fixed=FALSE) > fruitfly$line2 <- sub("RS", "S", fruitfly$line2, ignore.case =FALSE, fixed=FALSE) > fruitfly line fecund line2 1 RS 12.8 S 2 RS 21.6 S 3 RS 14.8 S 4 RS 23.1 S 5 RS 34.6 S 6 RS 19.7 S 7 RS 22.6 S 8 RS 29.6 S 9 RS 16.4 S 10 RS 20.3 S 11 RS 29.3 S 12 RS 14.9 S 13 RS 27.3 S 14 RS 22.4 S 15 RS 27.5 S 16 RS 20.3 S 17 RS 38.7 S 18 RS 26.4 S 19 RS 23.7 S 20 RS 26.1 S 21 RS 29.5 S 22 RS 38.6 S 23 RS 44.4 S 24 RS 23.2 S 25 RS 23.6 S 26 SS 38.4 S 27 SS 32.9 S 28 SS 48.5 S 29 SS 20.9 S 30 SS 11.6 S 31 SS 22.3 S 32 SS 30.2 S 33 SS 33.4 S 34 SS 26.7 S 35 SS 39.0 S 36 SS 12.8 S 37 SS 14.6 S 38 SS 12.2 S 39 SS 23.1 S 40 SS 29.4 S 41 SS 16.0 S 42 SS 20.1 S 43 SS 23.3 S 44 SS 22.9 S 45 SS 22.5 S 46 SS 15.1 S 47 SS 31.0 S 48 SS 16.9 S 49 SS 16.1 S 50 SS 10.8 S 51 NS 35.4 NS 52 NS 27.4 NS 53 NS 19.3 NS 54 NS 41.8 NS 55 NS 20.3 NS 56 NS 37.6 NS 57 NS 36.9 NS 58 NS 37.3 NS 59 NS 28.2 NS 60 NS 23.4 NS 61 NS 33.7 NS 62 NS 29.2 NS 63 NS 41.7 NS 64 NS 22.6 NS 65 NS 40.4 NS 66 NS 34.4 NS 67 NS 30.4 NS 68 NS 14.9 NS 69 NS 51.8 NS 70 NS 33.8 NS 71 NS 37.9 NS 72 NS 29.5 NS 73 NS 42.4 NS 74 NS 36.6 NS 75 NS 47.4 NS > attach(fruitfly) > # get summary info and check assumptions > by(fecund, line, function(x) summary(x)) line: NS Min. 1st Qu. Median Mean 3rd Qu. Max. 14.90 28.20 34.40 33.37 37.90 51.80 ----------------------------------------------- line: RS Min. 1st Qu. Median Mean 3rd Qu. Max. 12.80 20.30 23.60 25.26 29.30 44.40 ----------------------------------------------- line: SS Min. 1st Qu. Median Mean 3rd Qu. Max. 10.80 16.00 22.50 23.63 30.20 48.50 > by(fecund, line, function(x) sd(x)) line: NS [1] 8.942013 ----------------------------------------------- line: RS [1] 7.772391 ----------------------------------------------- line: SS [1] 9.768466 > boxplot(fecund~line) > by(fecund, line, function(x) stem(x)) The decimal point is 1 digit(s) to the right of the | 1 | 59 2 | 033789 3 | 00444577788 4 | 02227 5 | 2 The decimal point is 1 digit(s) to the right of the | 1 | 3556 2 | 000223334466789 3 | 00599 4 | 4 The decimal point is 1 digit(s) to the right of the | 1 | 122355667 2 | 012333379 3 | 013389 4 | 9 line: NS NULL ----------------------------------------------- line: RS NULL ----------------------------------------------- line: SS NULL > qqnorm(fecund[line=='NS'], main="Normal Q-Q Plot, NS line") > qqline(fecund[line=='NS']) > shapiro.test(fecund[line=='NS']) Shapiro-Wilk normality test data: fecund[line == "NS"] W = 0.9839, p-value = 0.9498 > qqnorm(fecund[line=='RS'], main="Normal Q-Q Plot, RS line") > qqline(fecund[line=='RS']) > shapiro.test(fecund[line=='RS']) Shapiro-Wilk normality test data: fecund[line == "RS"] W = 0.9496, p-value = 0.245 > qqnorm(fecund[line=='SS'], main="Normal Q-Q Plot, SS line") > qqline(fecund[line=='SS']) > shapiro.test(fecund[line=='SS']) Shapiro-Wilk normality test data: fecund[line == "SS"] W = 0.9396, p-value = 0.1446 > # model comparison F tests > anova_model1 <- aov(fecund~line) > summary(anova_model1) Df Sum Sq Mean Sq F value Pr(>F) line 2 1362 681.1 8.666 0.000424 *** Residuals 72 5659 78.6 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova_model2 <- aov(fecund~line2) > summary(anova_model2) Df Sum Sq Mean Sq F value Pr(>F) line2 1 1329 1329 17.05 9.59e-05 *** Residuals 73 5692 78 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(anova_model1,anova_model2) Analysis of Variance Table Model 1: fecund ~ line Model 2: fecund ~ line2 Res.Df RSS Df Sum of Sq F Pr(>F) 1 72 5659.0 2 73 5692.2 -1 -33.13 0.4215 0.5182 > # additional commands to estimate contrasts > # create indicator variables for NS, RS, NS, S > fruitfly$xNS <- as.double(fruitfly$line == "NS") > fruitfly$xRS <- as.double(fruitfly$line == "RS") > fruitfly$xSS <- as.double(fruitfly$line == "SS") > fruitfly$xS <- as.double(fruitfly$line2 == "S") > fruitfly$x <- fruitfly$xRS + fruitfly$xSS > attach(fruitfly) The following objects are masked from fruitfly (position 3): fecund, line, line2 > # model1 cell means model 3 means: NS, RS, SS (-1 requests no intercept) > fit_model1 <- lm(fecund ~ xNS + xRS + xSS -1) > predict(fit_model1, se.fit=TRUE, interval="confidence") $fit fit lwr upr 1 25.256 21.72138 28.79062 2 25.256 21.72138 28.79062 3 25.256 21.72138 28.79062 4 25.256 21.72138 28.79062 5 25.256 21.72138 28.79062 6 25.256 21.72138 28.79062 7 25.256 21.72138 28.79062 8 25.256 21.72138 28.79062 9 25.256 21.72138 28.79062 10 25.256 21.72138 28.79062 11 25.256 21.72138 28.79062 12 25.256 21.72138 28.79062 13 25.256 21.72138 28.79062 14 25.256 21.72138 28.79062 15 25.256 21.72138 28.79062 16 25.256 21.72138 28.79062 17 25.256 21.72138 28.79062 18 25.256 21.72138 28.79062 19 25.256 21.72138 28.79062 20 25.256 21.72138 28.79062 21 25.256 21.72138 28.79062 22 25.256 21.72138 28.79062 23 25.256 21.72138 28.79062 24 25.256 21.72138 28.79062 25 25.256 21.72138 28.79062 26 23.628 20.09338 27.16262 27 23.628 20.09338 27.16262 28 23.628 20.09338 27.16262 29 23.628 20.09338 27.16262 30 23.628 20.09338 27.16262 31 23.628 20.09338 27.16262 32 23.628 20.09338 27.16262 33 23.628 20.09338 27.16262 34 23.628 20.09338 27.16262 35 23.628 20.09338 27.16262 36 23.628 20.09338 27.16262 37 23.628 20.09338 27.16262 38 23.628 20.09338 27.16262 39 23.628 20.09338 27.16262 40 23.628 20.09338 27.16262 41 23.628 20.09338 27.16262 42 23.628 20.09338 27.16262 43 23.628 20.09338 27.16262 44 23.628 20.09338 27.16262 45 23.628 20.09338 27.16262 46 23.628 20.09338 27.16262 47 23.628 20.09338 27.16262 48 23.628 20.09338 27.16262 49 23.628 20.09338 27.16262 50 23.628 20.09338 27.16262 51 33.372 29.83738 36.90662 52 33.372 29.83738 36.90662 53 33.372 29.83738 36.90662 54 33.372 29.83738 36.90662 55 33.372 29.83738 36.90662 56 33.372 29.83738 36.90662 57 33.372 29.83738 36.90662 58 33.372 29.83738 36.90662 59 33.372 29.83738 36.90662 60 33.372 29.83738 36.90662 61 33.372 29.83738 36.90662 62 33.372 29.83738 36.90662 63 33.372 29.83738 36.90662 64 33.372 29.83738 36.90662 65 33.372 29.83738 36.90662 66 33.372 29.83738 36.90662 67 33.372 29.83738 36.90662 68 33.372 29.83738 36.90662 69 33.372 29.83738 36.90662 70 33.372 29.83738 36.90662 71 33.372 29.83738 36.90662 72 33.372 29.83738 36.90662 73 33.372 29.83738 36.90662 74 33.372 29.83738 36.90662 75 33.372 29.83738 36.90662 $se.fit [1] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [7] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [13] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [19] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [25] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [31] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [37] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [43] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [49] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [55] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [61] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [67] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [73] 1.773105 1.773105 1.773105 $df [1] 72 $residual.scale [1] 8.865525 > # estimate the contrast RS vs SS > out_model1 <- data.frame(xNS=c(0), xRS=c(1), xSS=c(-1)) > estimate_RSvsSS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_RSvsSS $fit fit lwr upr 1 1.628 -3.370708 6.626708 $se.fit [1] 2.507549 $df [1] 72 $residual.scale [1] 8.865525 > # model2 cell means model 2 means: NS and S (-1 requests no intercept) > fit_model2 <- lm(fecund ~ xNS + xS -1) > predict(fit_model2, se.fit=TRUE, interval="confidence") $fit fit lwr upr 1 24.442 21.95315 26.93085 2 24.442 21.95315 26.93085 3 24.442 21.95315 26.93085 4 24.442 21.95315 26.93085 5 24.442 21.95315 26.93085 6 24.442 21.95315 26.93085 7 24.442 21.95315 26.93085 8 24.442 21.95315 26.93085 9 24.442 21.95315 26.93085 10 24.442 21.95315 26.93085 11 24.442 21.95315 26.93085 12 24.442 21.95315 26.93085 13 24.442 21.95315 26.93085 14 24.442 21.95315 26.93085 15 24.442 21.95315 26.93085 16 24.442 21.95315 26.93085 17 24.442 21.95315 26.93085 18 24.442 21.95315 26.93085 19 24.442 21.95315 26.93085 20 24.442 21.95315 26.93085 21 24.442 21.95315 26.93085 22 24.442 21.95315 26.93085 23 24.442 21.95315 26.93085 24 24.442 21.95315 26.93085 25 24.442 21.95315 26.93085 26 24.442 21.95315 26.93085 27 24.442 21.95315 26.93085 28 24.442 21.95315 26.93085 29 24.442 21.95315 26.93085 30 24.442 21.95315 26.93085 31 24.442 21.95315 26.93085 32 24.442 21.95315 26.93085 33 24.442 21.95315 26.93085 34 24.442 21.95315 26.93085 35 24.442 21.95315 26.93085 36 24.442 21.95315 26.93085 37 24.442 21.95315 26.93085 38 24.442 21.95315 26.93085 39 24.442 21.95315 26.93085 40 24.442 21.95315 26.93085 41 24.442 21.95315 26.93085 42 24.442 21.95315 26.93085 43 24.442 21.95315 26.93085 44 24.442 21.95315 26.93085 45 24.442 21.95315 26.93085 46 24.442 21.95315 26.93085 47 24.442 21.95315 26.93085 48 24.442 21.95315 26.93085 49 24.442 21.95315 26.93085 50 24.442 21.95315 26.93085 51 33.372 29.85224 36.89176 52 33.372 29.85224 36.89176 53 33.372 29.85224 36.89176 54 33.372 29.85224 36.89176 55 33.372 29.85224 36.89176 56 33.372 29.85224 36.89176 57 33.372 29.85224 36.89176 58 33.372 29.85224 36.89176 59 33.372 29.85224 36.89176 60 33.372 29.85224 36.89176 61 33.372 29.85224 36.89176 62 33.372 29.85224 36.89176 63 33.372 29.85224 36.89176 64 33.372 29.85224 36.89176 65 33.372 29.85224 36.89176 66 33.372 29.85224 36.89176 67 33.372 29.85224 36.89176 68 33.372 29.85224 36.89176 69 33.372 29.85224 36.89176 70 33.372 29.85224 36.89176 71 33.372 29.85224 36.89176 72 33.372 29.85224 36.89176 73 33.372 29.85224 36.89176 74 33.372 29.85224 36.89176 75 33.372 29.85224 36.89176 $se.fit [1] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [7] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [13] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [19] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [25] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [31] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [37] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [43] 1.248797 1.248797 1.248797 1.248797 1.248797 1.248797 [49] 1.248797 1.248797 1.766066 1.766066 1.766066 1.766066 [55] 1.766066 1.766066 1.766066 1.766066 1.766066 1.766066 [61] 1.766066 1.766066 1.766066 1.766066 1.766066 1.766066 [67] 1.766066 1.766066 1.766066 1.766066 1.766066 1.766066 [73] 1.766066 1.766066 1.766066 $df [1] 73 $residual.scale [1] 8.830328 > # estimate the contrast NS vs S > out_model2 <- data.frame(xNS=c(1), xS=c(-1)) > estimate_NSvsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") > estimate_NSvsS $fit fit lwr upr 1 8.93 4.619188 13.24081 $se.fit [1] 2.16298 $df [1] 73 $residual.scale [1] 8.830328 > # --------------------------------------------------------------------------------------- > # alternate approach: model1 cell means model 3 means: NS, RS, SS (-1 requests no intercept) > fit_model1 <- lm(fecund ~ xNS + xRS + xSS -1) > predict(fit_model1, se.fit=TRUE, interval="confidence") $fit fit lwr upr 1 25.256 21.72138 28.79062 2 25.256 21.72138 28.79062 3 25.256 21.72138 28.79062 4 25.256 21.72138 28.79062 5 25.256 21.72138 28.79062 6 25.256 21.72138 28.79062 7 25.256 21.72138 28.79062 8 25.256 21.72138 28.79062 9 25.256 21.72138 28.79062 10 25.256 21.72138 28.79062 11 25.256 21.72138 28.79062 12 25.256 21.72138 28.79062 13 25.256 21.72138 28.79062 14 25.256 21.72138 28.79062 15 25.256 21.72138 28.79062 16 25.256 21.72138 28.79062 17 25.256 21.72138 28.79062 18 25.256 21.72138 28.79062 19 25.256 21.72138 28.79062 20 25.256 21.72138 28.79062 21 25.256 21.72138 28.79062 22 25.256 21.72138 28.79062 23 25.256 21.72138 28.79062 24 25.256 21.72138 28.79062 25 25.256 21.72138 28.79062 26 23.628 20.09338 27.16262 27 23.628 20.09338 27.16262 28 23.628 20.09338 27.16262 29 23.628 20.09338 27.16262 30 23.628 20.09338 27.16262 31 23.628 20.09338 27.16262 32 23.628 20.09338 27.16262 33 23.628 20.09338 27.16262 34 23.628 20.09338 27.16262 35 23.628 20.09338 27.16262 36 23.628 20.09338 27.16262 37 23.628 20.09338 27.16262 38 23.628 20.09338 27.16262 39 23.628 20.09338 27.16262 40 23.628 20.09338 27.16262 41 23.628 20.09338 27.16262 42 23.628 20.09338 27.16262 43 23.628 20.09338 27.16262 44 23.628 20.09338 27.16262 45 23.628 20.09338 27.16262 46 23.628 20.09338 27.16262 47 23.628 20.09338 27.16262 48 23.628 20.09338 27.16262 49 23.628 20.09338 27.16262 50 23.628 20.09338 27.16262 51 33.372 29.83738 36.90662 52 33.372 29.83738 36.90662 53 33.372 29.83738 36.90662 54 33.372 29.83738 36.90662 55 33.372 29.83738 36.90662 56 33.372 29.83738 36.90662 57 33.372 29.83738 36.90662 58 33.372 29.83738 36.90662 59 33.372 29.83738 36.90662 60 33.372 29.83738 36.90662 61 33.372 29.83738 36.90662 62 33.372 29.83738 36.90662 63 33.372 29.83738 36.90662 64 33.372 29.83738 36.90662 65 33.372 29.83738 36.90662 66 33.372 29.83738 36.90662 67 33.372 29.83738 36.90662 68 33.372 29.83738 36.90662 69 33.372 29.83738 36.90662 70 33.372 29.83738 36.90662 71 33.372 29.83738 36.90662 72 33.372 29.83738 36.90662 73 33.372 29.83738 36.90662 74 33.372 29.83738 36.90662 75 33.372 29.83738 36.90662 $se.fit [1] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [7] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [13] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [19] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [25] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [31] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [37] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [43] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [49] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [55] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [61] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [67] 1.773105 1.773105 1.773105 1.773105 1.773105 1.773105 [73] 1.773105 1.773105 1.773105 $df [1] 72 $residual.scale [1] 8.865525 > # estimate the contrast RS vs SS > out_model1 <- data.frame(xNS=c(0), xRS=c(1), xSS=c(-1)) > estimate_RSvsSS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_RSvsSS $fit fit lwr upr 1 1.628 -3.370708 6.626708 $se.fit [1] 2.507549 $df [1] 72 $residual.scale [1] 8.865525 > # estimate the contrast RS vs NS > out_model1 <- data.frame(xNS=c(-1), xRS=c(1), xSS=c(0)) > estimate_RSvsNS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_RSvsNS $fit fit lwr upr 1 -8.116 -13.11471 -3.117292 $se.fit [1] 2.507549 $df [1] 72 $residual.scale [1] 8.865525 > # estimate the contrast SS vs NS > out_model1 <- data.frame(xNS=c(-1), xRS=c(0), xSS=c(1)) > estimate_SSvsNS <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_SSvsNS $fit fit lwr upr 1 -9.744 -14.74271 -4.745292 $se.fit [1] 2.507549 $df [1] 72 $residual.scale [1] 8.865525 > # estimate the contrast NS vs others > out_model1 <- data.frame(xNS=c(1), xRS=c(-.5), xSS=c(-.5)) > estimate_NSvsothers <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_NSvsothers $fit fit lwr upr 1 8.93 4.600992 13.25901 $se.fit [1] 2.171601 $df [1] 72 $residual.scale [1] 8.865525 > # ------------------------------------------------- > # Scheffe simultaneous confidence intervals > # ------------------------------------------------- > estimate <- as.data.frame(rbind(estimate_RSvsSS$fit,estimate_RSvsNS$fit,estimate_SSvsNS$fit,estimate_NSvsothers$fit)) > stderr <- as.data.frame(rbind(estimate_RSvsSS$se.fit,estimate_RSvsNS$se.fit,estimate_SSvsNS$se.fit,estimate_NSvsothers$se.fit)) > multi <- sqrt(2*qf(.95,2,72)) > lowerCL <- estimate$fit - multi*stderr > upperCL <- estimate$fit + multi*stderr > scheffe_intervals <- cbind(estimate$fit,stderr,lowerCL,upperCL) > colnames(scheffe_intervals) <- c("estimate", "stderr", "lowerCL", "upperCL") > rownames(scheffe_intervals) <- c("RSvsSS", "RSvsNS", "SSvsNS", "NSvsothers") > scheffe_intervals estimate stderr lowerCL upperCL RSvsSS 1.628 2.507549 -4.639777 7.895777 RSvsNS -8.116 2.507549 -14.383777 -1.848223 SSvsNS -9.744 2.507549 -16.011777 -3.476223 NSvsothers 8.930 2.171601 3.501946 14.358054 >