library(lme4)
## Loading required package: Matrix
library(DescTools)
library(metRology)
## Warning: package 'metRology' was built under R version 4.1.3
##
## Attaching package: 'metRology'
## The following objects are masked from 'package:base':
##
## cbind, rbind
setwd("C:/Users/lhomm/OneDrive/Documents/R")
Dat <- read.table('https://users.stat.ufl.edu/~rrandles/sta4210/Rclassnotes/data/textdatasets/KutnerData/Chapter%2025%20Data%20Sets/CH25PR17.txt')
colnames(Dat) <- c("Value", "Coat", "Batch", "Count")
attach(Dat)
Coat <- as.factor(Dat$Coat)
Batch <- as.factor(Dat$Batch)
levels(Coat) <- c("6Coats", "8Coats", "10Coats")
batch = as.factor(Batch)
levels(Batch) <- c("Batch1", "Batch2", "Batch3", "Batch ")
View(Dat)
Decsion_Rule <- qf(0.05, 5, 36, lower.tail=FALSE) ### Decision Rule ###
Decsion_Rule
## [1] 2.477169
Anova1 <- aov(Value ~ Coat + Batch + Coat*Batch)
summary(Anova1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Coat 2 150.39 75.19 15.591 1.33e-05 ***
## Batch 3 152.85 50.95 10.564 3.98e-05 ***
## Coat:Batch 6 1.85 0.31 0.064 0.999
## Residuals 36 173.62 4.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Interaction is not significant which makes sense since batch is a random effect. The P-value is large and the F-value is smaller than our decision rule. ###
Anova2 <- aov(Value ~ Coat + Batch)
summary(Anova2)
## Df Sum Sq Mean Sq F value Pr(>F)
## Coat 2 150.4 75.19 18.0 2.26e-06 ***
## Batch 3 152.8 50.95 12.2 7.12e-06 ***
## Residuals 42 175.5 4.18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Decsion_Rule <- qf(0.05, 2, 42, lower.tail=FALSE) ### Decision Rule ###
Decsion_Rule
## [1] 3.219942
### As expected Coat is significant as it has a small P-value and a large F-value. What is surprising is that our random effect batch is signifigant. We will try a different model. ###
Anova3 <- lmer(Value ~ Coat + (1|Batch))
summary(Anova3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Value ~ Coat + (1 | Batch)
##
## REML criterion at convergence: 207.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1991 -0.6281 0.1033 0.6555 1.3638
##
## Random effects:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 3.898 1.974
## Residual 4.178 2.044
## Number of obs: 48, groups: Batch, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 73.1062 1.1116 65.769
## Coat8Coats 3.6875 0.7227 5.103
## Coat10Coats 3.8188 0.7227 5.284
##
## Correlation of Fixed Effects:
## (Intr) Ct8Cts
## Coat8Coats -0.325
## Coat10Coats -0.325 0.500
### This suggests a strong relationship between Coat and Value. ###
Bonferroni1 <- PostHocTest(Anova2, method = "bonferroni", conf.level = .9, which = "Coat")
Bonferroni1
##
## Posthoc multiple comparisons of means : Bonferroni
## 90% family-wise confidence level
##
## $Coat
## diff lwr.ci upr.ci pval
## 8Coats-6Coats 3.68750 2.097358 5.277642 2.3e-05 ***
## 10Coats-6Coats 3.81875 2.228608 5.408892 1.3e-05 ***
## 10Coats-8Coats 0.13125 -1.458892 1.721392 1.0000
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### There is a signifigant difference between a pearl with an 8 layer coat and those pearls with 10 layers and 6 layers but no difference between a pearl with a 6 and 8 layer coat. ###
Bonferroni2 <- PostHocTest(Anova2, method = "bonferroni", conf.level = .9, which = "Batch")
Bonferroni2
##
## Posthoc multiple comparisons of means : Bonferroni
## 90% family-wise confidence level
##
## $Batch
## diff lwr.ci upr.ci pval
## Batch2-Batch1 2.7083333 0.6274465 4.789220 0.01383 *
## Batch3-Batch1 3.4250000 1.3441132 5.505887 0.00110 **
## Batch -Batch1 -0.8333333 -2.9142201 1.247553 1.00000
## Batch3-Batch2 0.7166667 -1.3642201 2.797553 1.00000
## Batch -Batch2 -3.5416667 -5.6225535 -1.460780 0.00071 ***
## Batch -Batch3 -4.2583333 -6.3392201 -2.177447 4.6e-05 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### There are siginifigant differnces between Batches: 2 and 1, 3 and 1, 4 and 2, and 4 and 3.
Bonferroni3 <- PostHocTest(Anova2, method = "bonferroni", conf.level = .9)
Bonferroni3
##
## Posthoc multiple comparisons of means : Bonferroni
## 90% family-wise confidence level
##
## $Coat
## diff lwr.ci upr.ci pval
## 8Coats-6Coats 3.68750 2.097358 5.277642 2.3e-05 ***
## 10Coats-6Coats 3.81875 2.228608 5.408892 1.3e-05 ***
## 10Coats-8Coats 0.13125 -1.458892 1.721392 1.0000
##
## $Batch
## diff lwr.ci upr.ci pval
## Batch2-Batch1 2.7083333 0.6274465 4.789220 0.01383 *
## Batch3-Batch1 3.4250000 1.3441132 5.505887 0.00110 **
## Batch -Batch1 -0.8333333 -2.9142201 1.247553 1.00000
## Batch3-Batch2 0.7166667 -1.3642201 2.797553 1.00000
## Batch -Batch2 -3.5416667 -5.6225535 -1.460780 0.00071 ***
## Batch -Batch3 -4.2583333 -6.3392201 -2.177447 4.6e-05 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Summary. ###
Mean1 <- (Anova1$fitted.values[1] + Anova1$fitted.values[5] + Anova1$fitted.values[9] +
Anova1$fitted.values[13])/4
Mean2 <- (Anova1$fitted.values[17] + Anova1$fitted.values[21] + Anova1$fitted.values[25] +
Anova1$fitted.values[29])/4
Mean3 <- (Anova1$fitted.values[33] + Anova1$fitted.values[37] + Anova1$fitted.values[41] +
Anova1$fitted.values[45])/4
Mean1
## 1
## 73.10625
Mean2
## 17
## 76.79375
Mean3
## 33
## 76.925
### T-critical value with 6 degrees of freedom and alpha = .975 is 2.44. And the SE is .1968 ###
### D1: 76.79375 - 73.10625 +/- 2.44 * .1968 = (3.2, 4.17)
### D2: 76.92500 − 76.79375 +/- 2.44 * .1968 = (-.35, .61)
### Conclusion: The difference between Mean2 and Mean1 was significant as the interval does not contain zero. The difference between Mean3 and Mean2, however, was not significant as the interval does contain zero