--- title: "STAT 516 Lec 06" format: html --- Tensile strength data ```{r} y <- c(68,63,65,71,66,66,126,128,133,107,110,116, 93,101,98,63,60,59,56,59,57,40,41,44) agg <- as.factor(rep(c(rep("B",3),rep("S",3)),4)) comp <- as.factor(c(rep("st",6),rep("r",6), rep("l",6),rep("vl",6))) ``` ```{r} cbind(y,agg,comp) ``` ```{r} lm_out <- lm(y ~ agg + comp + agg:comp) lm_out aggregate(y ~ agg + comp, FUN = mean) ``` ```{r} summary(lm_out) ``` ```{r} plot(lm_out, which = 2) plot(lm_out,which = 1, add.smooth=F) ``` Get the ANOVA table quantities: ```{r} a <- 2 b <- 4 n <- 3 y... <- mean(y) yij. <- predict(lm_out) # returns fitted value, in this case the trt group means SStot <- sum((y - y...)^2) SSerror <- sum((y - yij.)^2) SStrt <- SStot - SSerror MSerror <- SSerror/( a*b*(n-1)) MStrt <- SStrt / (a*b - 1) Ftest <- MStrt / MSerror alpha <- 0.05 Fcrit <- qf(1-alpha,a*b - 1, a*b*(n-1)) pval <- 1 - pf(Ftest, a*b - 1, a*b*(n-1)) ``` The `summary()` output prints the F test statistic and p-value: ```{r} summary(lm_out) ``` Use the `anova()` function on the `lm()` output: ```{r} anova(lm_out) ``` ```{r} interaction.plot(agg,comp,y) interaction.plot(comp,agg,y) ``` Construct Tukey's CIs for all differences in pairs of means: ```{r} TukeyHSD(aov(lm_out)) ``` Here is how we can do Dunnett's to compare all treatments with a baseline treatment. Say the baseline is basalt aggregate and static compaction. ```{r} agg_comp <- as.factor(paste(agg,comp,sep="x")) library(DescTools) DunnettTest(y ~ agg_comp, control="Bxst", conf.level=0.95) ``` Serum glucose data: ```{r} y <- c(42.5,43.3,42.9,138.4,144.4,142.7,180.9,180.5,183.0, 39.8,40.3,41.2,132.4,132.4,130.3,176.8,173.6,174.9) method <- as.factor(c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2)) glc <- as.factor(c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3)) ``` ```{r} stripchart(y ~ method + glc, method = "jitter", vertical = T) ``` Look at interaction plots: ```{r} interaction.plot(method, glc, log(y)) interaction.plot(glc,method, log(y)) ```