setwd("C:/Users/lhomm/OneDrive/Documents/R")
library(asbio)
## Loading required package: tcltk
library(DescTools)
##
## Attaching package: 'DescTools'
## The following object is masked from 'package:asbio':
##
## Mode
Dat <- read.table ("http://users.stat.ufl.edu/~rrandles/sta4210/Rclassnotes/data/textdatasets/KutnerData/Chapter%2027%20Data%20Sets/CH27PR03.txt")
colnames (Dat) <- c("Pressure", "Rabbit", "Dose" )
Dose<- factor(Dat$Dose)
levels(Dose) <- c("0.1", "0.3", "0.5", "1.0", "1.5", "3.0")
Rabbit <- factor(Dat$Rabbit, levels = 1:12)
Anova1 <- aov(Dat$Pressure ~ Dose + Rabbit)
Anova1
## Call:
## aov(formula = Dat$Pressure ~ Dose + Rabbit)
##
## Terms:
## Dose Rabbit Residuals
## Sum of Squares 5826.278 1197.444 467.389
## Deg. of Freedom 5 11 55
##
## Residual standard error: 2.915129
## Estimated effects may be unbalanced
summary(Anova1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Dose 5 5826 1165.3 137.12 < 2e-16 ***
## Rabbit 11 1197 108.9 12.81 1.43e-11 ***
## Residuals 55 467 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Anova1$residuals, Anova1$fitted.values)
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.23738, p-value = 0.8124
##
##
## $Independence[[2]]
## lag Autocorrelation D-W Statistic p-value
## 1 -0.0363987 1.981086 0.184
## Alternative hypothesis: rho != 0
##
##
## $Normality
##
## Shapiro-Wilk normality test
##
## data: re
## W = 0.98682, p-value = 0.6578
##
##
## [[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 5 0.2533 0.9367
## 66
### It seems that the block dose has an effect on blood pressure but not the individual rabbit. This is fine as we are not particularly interested in the rabbits but rather the underlying population. It also appears that the data is roughly normal although there may be some problems in the tails. ###
source("https://raw.githubusercontent.com/athienit/STA4210material/main/check.R")
check(Anova1,tests=TRUE)
## $Independence
## $Independence[[1]]
##
## Runs Test - Two sided
##
## data: re
## Standardized Runs Statistic = 0.23738, p-value = 0.8124
##
##
## $Independence[[2]]
## lag Autocorrelation D-W Statistic p-value
## 1 -0.0363987 1.981086 0.164
## Alternative hypothesis: rho != 0
##
##
## $Normality
##
## Shapiro-Wilk normality test
##
## data: re
## W = 0.98682, p-value = 0.6578
##
##
## [[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 5 0.2533 0.9367
## 66
### The plot of the residuals and the Levees test suggest that the assumption of equal variances has not been violated. ###
interaction.plot(Dose, Rabbit, Dat$Pressure, xlab = "Dose", ylab = "Pressure",
main = "Interaction Plot")
### The assumption of no interaction seems okay as most of the lines are parallel and overlapping not crossing or diverging. ###
### Ho: D the difference is equal to 0. ###
### H1: D the difference is not equal to 0. ###
Decsion_Rule <- qf(0.005, 1, 54, lower.tail=FALSE) ### Decision Rule ###
Decsion_Rule
## [1] 8.567124
tukey.add.test(Dat$Pressure, Dose, Rabbit)
##
## Tukey's one df test for additivity
## F = 1.0319384 Denom df = 54 p-value = 0.3142338
### With an F-Value smaller than our decision rule and a P-value larger than .0/.05/.1 we reject Ho: that the differences is not equal to zero. ###
### Yes we should treat Rabbit as random. We are not interested in the individual rabbits. It is also not entirely clear what we are interested. Is the underlying population all rabbits, all of one type of rabbit, or are the rabbits just meant to act as test subjects in which case the rabbits are almost eniterly inconsequential. All of this is to say that it is unlikely that we would care at all about the rabbits themselves. ###
Anova1
## Call:
## aov(formula = Dat$Pressure ~ Dose + Rabbit)
##
## Terms:
## Dose Rabbit Residuals
## Sum of Squares 5826.278 1197.444 467.389
## Deg. of Freedom 5 11 55
##
## Residual standard error: 2.915129
## Estimated effects may be unbalanced
summary(Anova1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Dose 5 5826 1165.3 137.12 < 2e-16 ***
## Rabbit 11 1197 108.9 12.81 1.43e-11 ***
## Residuals 55 467 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Decsion_Rule <- qf(0.01, 5, 55, lower.tail=FALSE) ### Decision Rule ###
Decsion_Rule
## [1] 3.369962
summary(Anova1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Dose 5 5826 1165.3 137.12 < 2e-16 ***
## Rabbit 11 1197 108.9 12.81 1.43e-11 ***
## Residuals 55 467 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### with an F-Value larger than our decision rule and a P-values smaller than .01/.05/.1. ###
Dose1 <- 21 + 19 + 12 + 9 + 7 + 18 + 9 + 20 + 18 + 8 +18 + 17
Dose1
## [1] 176
Y.hat_1 <- Dose1/ 12
Y.hat_1
## [1] 14.66667
Dose2 <- 21 + 24 + 25 + 17 +10 + 26 + 12 +20 + 18 + 12 + 22 + 23
Dose2
## [1] 230
Y.hat_2 <- Dose2 / 12
Y.hat_2
## [1] 19.16667
Dose3 <- 23 + 27 + 27 + 18 + 19 + 26 + 17 + 30 + 27 + 11 +25 + 26
Dose3
## [1] 276
Y.hat_3 <- Dose3 / 12
Y.hat_3
## [1] 23
Dose4 <- 35 + 36 + 26 + 27 + 25 + 29 + 22 + 30 + 31 + 24 + 32 + 28
Dose4
## [1] 345
Y.hat_4 <- Dose4 / 12
Y.hat_4
## [1] 28.75
Dose5 <- 36 + 36 + 33 + 34 + 31 + 39 + 33 + 38 + 42 + 26 + 38 + 34
Dose5
## [1] 420
Y.hat_5 <- Dose5 / 12
Y.hat_5
## [1] 35
Dose6 <- 48 + 46 + 40 + 39 + 38 + 44 + 40 + 41 + 49 + 31 + 38 + 35
Dose6
## [1] 489
Y.hat_6 <- Dose6 / 12
Y.hat_6
## [1] 40.75
L1 <- Y.hat_1 - Y.hat_2
L1
## [1] -4.5
L2 <- Y.hat_2 - Y.hat_3
L2
## [1] -3.833333
L3 <- Y.hat_3 - Y.hat_4
L3
## [1] -5.75
L4 <- Y.hat_4 - Y.hat_5
L4
## [1] -6.25
L5 <- Y.hat_5 - Y.hat_6
L5
## [1] -5.75
### The contrasts of interest from the list are 0.3-0.1, 0.5-0.3, 1.0-.0.5, 1.5 -1.0, and 3.0 -1.5 ###
### The contrasts are positives due to the differences in the way the function sets up the contrasts. ###
Bonferroni <- PostHocTest(Anova1, method = "bonferroni", conf.level = .95, which = "Dose")
Bonferroni
##
## Posthoc multiple comparisons of means : Bonferroni
## 95% family-wise confidence level
##
## $Dose
## diff lwr.ci upr.ci pval
## 0.3-0.1 4.500000 0.8478034 8.152197 0.00580 **
## 0.5-0.1 8.333333 4.6811367 11.985530 5.6e-08 ***
## 1.0-0.1 14.083333 10.4311367 17.735530 1.4e-15 ***
## 1.5-0.1 20.333333 16.6811367 23.985530 < 2e-16 ***
## 3.0-0.1 26.083333 22.4311367 29.735530 < 2e-16 ***
## 0.5-0.3 3.833333 0.1811367 7.485530 0.03220 *
## 1.0-0.3 9.583333 5.9311367 13.235530 1.1e-09 ***
## 1.5-0.3 15.833333 12.1811367 19.485530 < 2e-16 ***
## 3.0-0.3 21.583333 17.9311367 25.235530 < 2e-16 ***
## 1.0-0.5 5.750000 2.0978034 9.402197 0.00017 ***
## 1.5-0.5 12.000000 8.3478034 15.652197 6.2e-13 ***
## 3.0-0.5 17.750000 14.0978034 21.402197 < 2e-16 ***
## 1.5-1.0 6.250000 2.5978034 9.902197 3.8e-05 ***
## 3.0-1.0 12.000000 8.3478034 15.652197 6.2e-13 ***
## 3.0-1.5 5.750000 2.0978034 9.402197 0.00017 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Because this function does all of the contrasts the individual error wise rate is different then if I did 6 contrasts with confidence level .95. I adjusted the family confidence to .85 to get interval closer to what they would be if the function did only 6 at alpha level .95. ###
Bonferroni1 <- PostHocTest(Anova1, method = "bonferroni", conf.level = .85, which = "Dose")
Bonferroni1
##
## Posthoc multiple comparisons of means : Bonferroni
## 85% family-wise confidence level
##
## $Dose
## diff lwr.ci upr.ci pval
## 0.3-0.1 4.500000 1.3245652 7.675435 0.00580 **
## 0.5-0.1 8.333333 5.1578985 11.508768 5.6e-08 ***
## 1.0-0.1 14.083333 10.9078985 17.258768 1.4e-15 ***
## 1.5-0.1 20.333333 17.1578985 23.508768 < 2e-16 ***
## 3.0-0.1 26.083333 22.9078985 29.258768 < 2e-16 ***
## 0.5-0.3 3.833333 0.6578985 7.008768 0.03220 *
## 1.0-0.3 9.583333 6.4078985 12.758768 1.1e-09 ***
## 1.5-0.3 15.833333 12.6578985 19.008768 < 2e-16 ***
## 3.0-0.3 21.583333 18.4078985 24.758768 < 2e-16 ***
## 1.0-0.5 5.750000 2.5745652 8.925435 0.00017 ***
## 1.5-0.5 12.000000 8.8245652 15.175435 6.2e-13 ***
## 3.0-0.5 17.750000 14.5745652 20.925435 < 2e-16 ***
## 1.5-1.0 6.250000 3.0745652 9.425435 3.8e-05 ***
## 3.0-1.0 12.000000 8.8245652 15.175435 6.2e-13 ***
## 3.0-1.5 5.750000 2.5745652 8.925435 0.00017 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Bonferroni)
plot(Bonferroni1)
### All of the contrasts are significant so it seems that there significant differences between the doses. ###
Anova1
## Call:
## aov(formula = Dat$Pressure ~ Dose + Rabbit)
##
## Terms:
## Dose Rabbit Residuals
## Sum of Squares 5826.278 1197.444 467.389
## Deg. of Freedom 5 11 55
##
## Residual standard error: 2.915129
## Estimated effects may be unbalanced
anova(Anova1)
## Analysis of Variance Table
##
## Response: Dat$Pressure
## Df Sum Sq Mean Sq F value Pr(>F)
## Dose 5 5826.3 1165.26 137.12 < 2.2e-16 ***
## Rabbit 11 1197.4 108.86 12.81 1.427e-11 ***
## Residuals 55 467.4 8.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DF <- anova(Anova1)[,"Df"]
DF
## [1] 5 11 55
MS <- anova(Anova1)[, "Mean Sq"]
MS
## [1] 1165.25556 108.85859 8.49798
(DF[2]*MS[2]+(DF[2]+1)*DF[1]*MS[3])/(sum(DF)*MS[3])
## [1] 2.829709
### Because we got a value larger than 1 it suggest that the CRBD is more efficient than the CRD. ###