setwd("E:/Methods")
library(asbio)
## Loading required package: tcltk
library(DescTools)
##
## Attaching package: 'DescTools'
## The following object is masked from 'package:asbio':
##
## Mode
library(ggplot2)
library(gplots)
## Registered S3 method overwritten by 'gplots':
## method from
## reorder.factor DescTools
##
## Attaching package: 'gplots'
## The following object is masked from 'package:DescTools':
##
## reorder.factor
## The following object is masked from 'package:stats':
##
## lowess
Dat <- read.table("http://users.stat.ufl.edu/~rrandles/sta4210/Rclassnotes/data/textdatasets/KutnerData/Chapter%2028%20Data%20Sets/CH28PR17.txt")
colnames(Dat) <- c("Help", "Exec", "Month", "Report")
attach(Dat)
Exec <- as.factor(Exec)
levels(Exec) = c("Harrison", "Smith", "Carmichael", "Loeb",
"Month")
Month <- as.factor(Month)
levels(Month) = c("March", "April", "May", "June", "July")
Report <- as.factor(Report)
levels(Report) = c("A", "B", "C", "D", "E")
Anova1 <- aov(Dat$Help ~ Exec + Month + Report)
Anova1
## Call:
## aov(formula = Dat$Help ~ Exec + Month + Report)
##
## Terms:
## Exec Month Report Residuals
## Sum of Squares 220.16 10.96 527.36 19.28
## Deg. of Freedom 4 4 4 12
##
## Residual standard error: 1.267544
## Estimated effects may be unbalanced
Plot1 <- plot(Anova1$residuals, Anova1$fitted.values, xlab = "Residuals", ylab = "Fitted Values", main = "Residuals Vs. Fitted Values")
Plot1
## NULL
Standardized1 <- rstandard(Anova1)
Standardized1
## 1 2 3 4 5 6 7
## -1.0020725 -0.7743288 1.0476213 0.3643900 0.3643900 0.3643900 0.1366463
## 8 9 10 11 12 13 14
## -0.3188413 1.0476213 -1.2298163 0.5921338 -0.7743288 -1.2298163 0.1366463
## 15 16 17 18 19 20 21
## 1.2753651 -0.7743288 2.1863401 0.5921338 -0.0910975 -1.9130476 0.8198775
## 22 23 24 25
## -0.7743288 -0.0910975 -1.4575601 1.5031088
qqnorm(Standardized1, datax = TRUE)
qqline(Standardized1, datax = TRUE)
Normal1 <- shapiro.test(Standardized1)
Normal1
##
## Shapiro-Wilk normality test
##
## data: Standardized1
## W = 0.98256, p-value = 0.9307
source("https://raw.githubusercontent.com/athienit/STA4210material/main/check.R")
check(Anova1,tests= TRUE)
## Loading required package: lawstat
## Loading required package: car
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:lawstat':
##
## levene.test
## The following object is masked from 'package:DescTools':
##
## Recode
## $Independence
## $Independence[[1]]
##
## Runs Test - Two sided
##
## data: re
## Standardized Runs Statistic = 0.21283, p-value = 0.8315
##
##
## $Independence[[2]]
## lag Autocorrelation D-W Statistic p-value
## 1 -0.2675519 2.404564 0.734
## Alternative hypothesis: rho != 0
##
##
## $Normality
##
## Shapiro-Wilk normality test
##
## data: re
## W = 0.98256, p-value = 0.9307
##
##
## [[3]]
## [1] "Constant Variance only valid if data are in groups"
##
## $ConstantVar
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.0136 0.9996
## 20
### It appears that no assumptions were violated. Both blocks should be considered random as we are not necessarily just interested in these few executives or these few month but rather the underlying population. Additionally this is a small sample of months and execs. ###
Plot2 <- plotmeans(Dat$Help ~ Report, xlab = "Summary Report", ylab = "Help Score", main = "Main Effects")
Plot2
## NULL
### This plot implies that there is such thing as too much detail in a summary report. The most detailed reports A amd B seem to provide too much detial and information to be helpful to the reader. It seems that D and to a lesser extent C and E provide enough information with out bogging the reader down in too much information. ###
### Ho: None of the reports are different from one antother ###
### H1: At least one of the reports is differt from one other. ###
Decsion_Rule <- qf(0.01, 4, 12, lower.tail = FALSE) ### Decision Rule ###
Decsion_Rule
## [1] 5.411951
summary(Anova1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Exec 4 220.2 55.04 34.257 1.78e-06 ***
## Month 4 11.0 2.74 1.705 0.213
## Report 4 527.4 131.84 82.058 1.31e-08 ***
## Residuals 12 19.3 1.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### With an F-Value larger than the decision rule and a P-value smaller than .01/.05/.1 we reject Ho: that none of the reports are different from one another. Instead we have Support for H1: that there is at least one report differt from one other. ###
Tukey1 <- PostHocTest(Anova1, method = "hsd", conf.level = .95, which = "Report")
Tukey1
##
## Posthoc multiple comparisons of means : Tukey HSD
## 95% family-wise confidence level
##
## $Report
## diff lwr.ci upr.ci pval
## B-A 0.4 -2.155253 2.9552527 0.98592
## C-A 8.0 5.444747 10.5552527 3.0e-06 ***
## D-A 12.0 9.444747 14.5552527 3.3e-08 ***
## E-A 6.4 3.844747 8.9552527 3.1e-05 ***
## C-B 7.6 5.044747 10.1552527 5.2e-06 ***
## D-B 11.6 9.044747 14.1552527 4.8e-08 ***
## E-B 6.0 3.444747 8.5552527 6.0e-05 ***
## D-C 4.0 1.444747 6.5552527 0.00237 **
## E-C -1.6 -4.155253 0.9552527 0.32396
## E-D -5.6 -8.155253 -3.0447473 0.00012 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### The results from the Pairwise comparisons are consitent with previous results. The only comparisons that were not signifigant were between A and B which is asgain consistent. ###
MSCOL <- 2.74
MSROW <- 55.04
r <- 5
r_1 <- 5 - 1
r_2 <- 5 + 1
MSE <- 1.61
E1 = ((MSCOL + MSROW + (r_1 * MSE)) / (r_2 * MSE))
E1
## [1] 6.648033
E2 = ((MSCOL + (r_1 * MSE)) / (r * MSE))
E2
## [1] 1.140373
E3 = ((MSROW + (r_1 * MSE)) / (r * MSE))
E3
## [1] 7.637267
df = anova(Anova1)[, "Df"]
MS = anova(Anova1)[, "Mean Sq"]
(MS[2] + MS[3] + df[1] * MS[4])/((df[1] + 2) * MS[4])
## [1] 14.62725
### The Latin Squares design method had the same power as a CRD but was more efficent as it was able to achiece the same power with 14.63 ~ 15 or 6.65 ~ 7 times fewer observations depending on which formula you use. ###