potato <- read.table("C:/Rexamples/potatodata.txt", header=T) potato$treat2 <- potato$treat potato$treat2 <- sub("fructose", "sixcarb", potato$treat2, ignore.case =FALSE, fixed=FALSE) potato$treat2 <- sub("glucose", "sixcarb", potato$treat2, ignore.case =FALSE, fixed=FALSE) potato$treat3 <- potato$treat2 potato$treat3 <- sub("sucrose", "sugar", potato$treat3, ignore.case =FALSE, fixed=FALSE) potato$treat3 <- sub("sixcarb", "sugar", potato$treat3, ignore.case =FALSE, fixed=FALSE) View(potato) summary(potato) attach(potato) potato anova_model1 <- aov(time~treat) summary(anova_model1) anova_model2 <- aov(time~treat2) summary(anova_model2) anova_model3 <- aov(time~treat3) summary(anova_model3) anova(anova_model1,anova_model2) anova(anova_model2,anova_model3) # additional commands to estimate contrasts # create indicator variables for c,f,g,s, and 6(6carb) potato$xc <- as.double(potato$treat == "control") potato$xf <- as.double(potato$treat == "fructose") potato$xg <- as.double(potato$treat == "glucose") potato$xs <- as.double(potato$treat == "sucrose") potato$x6 <- potato$xf + potato$xg attach(potato) # model1 cell means model 4 means: c,f,g,s (-1 requests no intercept) fit_model1 <- lm(time ~ xc + xf + xg + xs -1) predict(fit_model1, se.fit=TRUE, interval="confidence") # estimate the contrast F vs G out_model1 <- data.frame(xc=c(0), xf=c(1), xg=c(-1), xs=c(0)) estimate_FvsG <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") estimate_FvsG # model2 cell means model 3 means: c,6,s (-1 requests no intercept) fit_model2 <- lm(time ~ xc + x6 + xs -1) predict(fit_model2, se.fit=TRUE, interval="confidence") # estimate the contrast 6 carb vs S out_model2 <- data.frame(xc=c(0), x6=c(1), xs=c(-1)) estimate_6vsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") estimate_6vsS # ------------------------------------------------- # more contrasts and Scheffe simultaneous confidence intervals # ------------------------------------------------- # estimate the contrast 6 carb vs C out_model2 <- data.frame(xc=c(-1), x6=c(1), xs=c(0)) estimate_6vsC <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") estimate_6vsC # estimate the contrast C vs S out_model2 <- data.frame(xc=c(1), x6=c(0), xs=c(-1)) estimate_CvsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") estimate_CvsS # estimate the contrast S vs (C and 6 carb) out_model2 <- data.frame(xc=c(-.5), x6=c(-.5), xs=c(1)) estimate_Svs6andC <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") estimate_Svs6andC estimate <- as.data.frame(rbind(estimate_CvsS$fit,estimate_6vsS$fit,estimate_6vsC$fit,estimate_Svs6andC$fit)) stderr <- as.data.frame(rbind(estimate_CvsS$se.fit,estimate_6vsS$se.fit,estimate_6vsC$se.fit,estimate_Svs6andC$se.fit)) multi <- sqrt(2*qf(.95,2,5)) 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("CvsS", "6vsS", "6vsC", "Svs6andC") scheffe_intervals --------------------- NOTE ------------------------ The first line above which reads the data from the file "potatodata.txt" assumes that this file is in the directory "C:/Rexamples". Obey this convention or modify the command. -------------------------- results ---------------------------------- The remainder of this file is copied and pasted from the RStudio command window. -------------------------------------------------------- > potato <- read.table("C:/Rexamples/potatodata.txt", header=T) > potato$treat2 <- potato$treat > potato$treat2 <- sub("fructose", "sixcarb", potato$treat2, ignore.case =FALSE, fixed=FALSE) > potato$treat2 <- sub("glucose", "sixcarb", potato$treat2, ignore.case =FALSE, fixed=FALSE) > potato$treat3 <- potato$treat2 > potato$treat3 <- sub("sucrose", "sugar", potato$treat3, ignore.case =FALSE, fixed=FALSE) > potato$treat3 <- sub("sixcarb", "sugar", potato$treat3, ignore.case =FALSE, fixed=FALSE) > View(potato) > summary(potato) treat time treat2 control :2 Min. :1.700 Length:8 fructose:2 1st Qu.:2.250 Class :character glucose :2 Median :2.550 Mode :character sucrose :2 Mean :2.725 3rd Qu.:3.150 Max. :4.000 treat3 Length:8 Class :character Mode :character > attach(potato) > potato treat time treat2 treat3 1 control 2.3 control control 2 control 1.7 control control 3 sucrose 3.6 sucrose sugar 4 sucrose 4.0 sucrose sugar 5 glucose 3.0 sixcarb sugar 6 glucose 2.8 sixcarb sugar 7 fructose 2.1 sixcarb sugar 8 fructose 2.3 sixcarb sugar > anova_model1 <- aov(time~treat) > summary(anova_model1) Df Sum Sq Mean Sq F value Pr(>F) treat 3 3.975 1.325 17.67 0.00901 ** Residuals 4 0.300 0.075 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova_model2 <- aov(time~treat2) > summary(anova_model2) Df Sum Sq Mean Sq F value Pr(>F) treat2 2 3.485 1.742 11.03 0.0147 * Residuals 5 0.790 0.158 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova_model3 <- aov(time~treat3) > summary(anova_model3) Df Sum Sq Mean Sq F value Pr(>F) treat3 1 1.402 1.4017 2.927 0.138 Residuals 6 2.873 0.4789 > anova(anova_model1,anova_model2) Analysis of Variance Table Model 1: time ~ treat Model 2: time ~ treat2 Res.Df RSS Df Sum of Sq F Pr(>F) 1 4 0.30 2 5 0.79 -1 -0.49 6.5333 0.0629 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(anova_model2,anova_model3) Analysis of Variance Table Model 1: time ~ treat2 Model 2: time ~ treat3 Res.Df RSS Df Sum of Sq F Pr(>F) 1 5 0.7900 2 6 2.8733 -1 -2.0833 13.186 0.01504 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > # additional commands to estimate contrasts > # create indicator variables for c,f,g,s, and 6(6carb) > potato$xc <- as.double(potato$treat == "control") > potato$xf <- as.double(potato$treat == "fructose") > potato$xg <- as.double(potato$treat == "glucose") > potato$xs <- as.double(potato$treat == "sucrose") > potato$x6 <- potato$xf + potato$xg > attach(potato) The following objects are masked from potato (position 3): time, treat, treat2, treat3, x6, xc, xf, xg, xs The following objects are masked from potato (position 4): time, treat, treat2, treat3 > # model1 cell means model 4 means: c,f,g,s (-1 requests no intercept) > fit_model1 <- lm(time ~ xc + xf + xg + xs -1) > predict(fit_model1, se.fit=TRUE, interval="confidence") $fit fit lwr upr 1 2.0 1.462344 2.537656 2 2.0 1.462344 2.537656 3 3.8 3.262344 4.337656 4 3.8 3.262344 4.337656 5 2.9 2.362344 3.437656 6 2.9 2.362344 3.437656 7 2.2 1.662344 2.737656 8 2.2 1.662344 2.737656 $se.fit [1] 0.1936492 0.1936492 0.1936492 0.1936492 0.1936492 0.1936492 [7] 0.1936492 0.1936492 $df [1] 4 $residual.scale [1] 0.2738613 > # estimate the contrast F vs G > out_model1 <- data.frame(xc=c(0), xf=c(1), xg=c(-1), xs=c(0)) > estimate_FvsG <- predict(fit_model1, out_model1, se.fit=TRUE, interval="confidence") > estimate_FvsG $fit fit lwr upr 1 -0.7 -1.460361 0.06036081 $se.fit [1] 0.2738613 $df [1] 4 $residual.scale [1] 0.2738613 > # model2 cell means model 3 means: c,6,s (-1 requests no intercept) > fit_model2 <- lm(time ~ xc + x6 + xs -1) > predict(fit_model2, se.fit=TRUE, interval="confidence") $fit fit lwr upr 1 2.00 1.277488 2.722512 2 2.00 1.277488 2.722512 3 3.80 3.077488 4.522512 4 3.80 3.077488 4.522512 5 2.55 2.039107 3.060893 6 2.55 2.039107 3.060893 7 2.55 2.039107 3.060893 8 2.55 2.039107 3.060893 $se.fit [1] 0.2810694 0.2810694 0.2810694 0.2810694 0.1987461 0.1987461 [7] 0.1987461 0.1987461 $df [1] 5 $residual.scale [1] 0.3974921 > # estimate the contrast 6 carb vs S > out_model2 <- data.frame(xc=c(0), x6=c(1), xs=c(-1)) > estimate_6vsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") > estimate_6vsS $fit fit lwr upr 1 -1.25 -2.134893 -0.3651073 $se.fit [1] 0.3442383 $df [1] 5 $residual.scale [1] 0.3974921 > > # ------------------------------------------------- > # more contrasts and Scheffe simultaneous confidence intervals > # ------------------------------------------------- > > # estimate the contrast 6 carb vs C > out_model2 <- data.frame(xc=c(-1), x6=c(1), xs=c(0)) > estimate_6vsC <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") > estimate_6vsC $fit fit lwr upr 1 0.55 -0.3348927 1.434893 $se.fit [1] 0.3442383 $df [1] 5 $residual.scale [1] 0.3974921 > # estimate the contrast C vs S > out_model2 <- data.frame(xc=c(1), x6=c(0), xs=c(-1)) > estimate_CvsS <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") > estimate_CvsS $fit fit lwr upr 1 -1.8 -2.821786 -0.7782139 $se.fit [1] 0.3974921 $df [1] 5 $residual.scale [1] 0.3974921 > # estimate the contrast S vs (C and 6 carb) > out_model2 <- data.frame(xc=c(-.5), x6=c(-.5), xs=c(1)) > estimate_Svs6andC <- predict(fit_model2, out_model2, se.fit=TRUE, interval="confidence") > estimate_Svs6andC $fit fit lwr upr 1 1.525 0.6777797 2.37222 $se.fit [1] 0.3295831 $df [1] 5 $residual.scale [1] 0.3974921 > estimate <- as.data.frame(rbind(estimate_CvsS$fit,estimate_6vsS$fit,estimate_6vsC$fit,estimate_Svs6andC$fit)) > stderr <- as.data.frame(rbind(estimate_CvsS$se.fit,estimate_6vsS$se.fit,estimate_6vsC$se.fit,estimate_Svs6andC$se.fit)) > multi <- sqrt(2*qf(.95,2,5)) > 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("CvsS", "6vsS", "6vsC", "Svs6andC") > scheffe_intervals estimate stderr lowerCL upperCL CvsS -1.800 0.3974921 -3.1521903 -0.44780968 6vsS -1.250 0.3442383 -2.4210312 -0.07896883 6vsC 0.550 0.3442383 -0.6210312 1.72103117 Svs6andC 1.525 0.3295831 0.4038230 2.64617699 > > >