In-Class-Activity-1
Liora Yehuday
2024-09-10
Dataset: fish_transport.csv
Source: K. Thongprajukaew, S. Takaeh, N. Esor, S. Saekhow, S. Malawa, N. Nuntapong, W.Hahor, A. Choodum (2023). “Optimal water volume for transportation of male Siamese fighting fish (Betta splendens),” Aquaculture Reports, Vol. 28, https://doi.org/10.1016/j.aqrep.2022.101430.
Description: Skin coloration scores for Siamese fighting fish when shipped in various water volumes. Responses: L* (0=darkness, 100=brightness) and a(+a=red, -a=green). Treatments: Control (1,not transported), 40mL(2), 60mL(3), 80mL(4), 100mL(5), 120mL(6), 140mL(7) n=15 fish per treatment in Completely Randomized Design (each fish in different container). Data generated to match means and SDs (paper reported SEMs: SEM=SD/sqrt(n) => SD=sqrt(n)SEM).
Variable Names: Y= shipTrt (labels given in Description) trt= trtID (Fish ID within Trt) L.star a.star
Read in data, print out top and bottom of data
Response = L.star
library(effectsize)
library(agricolae)
ft <- read.csv("http://www.stat.ufl.edu/~winner/data/fish_transport.csv")
head(ft); tail(ft)## shipTrt trtID L.star a.star
## 1 1 1 17.7575 24.4971
## 2 1 2 17.9925 26.9582
## 3 1 3 18.4705 27.6317
## 4 1 4 13.4241 18.0226
## 5 1 5 19.7282 25.7558
## 6 1 6 13.1565 28.1487## shipTrt trtID L.star a.star
## 100 7 10 28.8744 15.1266
## 101 7 11 26.1500 14.0599
## 102 7 12 29.7907 16.2409
## 103 7 13 29.9627 14.3175
## 104 7 14 21.6483 14.0810
## 105 7 15 28.4363 13.6530attach(ft)Create factor variable for Water Volume Treatments and plot data, Trt Means, Overall Mean, Trt SD’s and n’s
shipTrt.f <- factor(shipTrt, levels=1:7,
labels=c("C","40mL", "60mL","80mL","100mL","120mL","140mL"))
all.mean <- mean(L.star)
trt.mean <- as.vector(tapply(L.star, shipTrt.f, mean))
trt.sd <- as.vector(tapply(L.star, shipTrt.f, sd))
trt.n <- as.vector(tapply(L.star, shipTrt.f, length))
n_T <- sum(trt.n)
r <- length(trt.n)
plot(L.star ~ shipTrt, pch=16, col="blue", cex=0.5, xaxt="n",
xlab="Treatment", ylab="L*", xlim=c(0.5,7.5))
axis(1,1:7,labels=c("Control","40mL","60mL","80mL","100mL","120mL","140mL"))
arrows(0.8, mean(L.star), 7.2, mean(L.star), code=3, col="red",
angle=0, lwd=1.5)
for(i1 in 1:7) {
arrows(i1-0.2,mean(L.star[shipTrt==i1]),i1+0.2,mean(L.star[shipTrt==i1]),
code=3, col="purple", angle=0, lwd=1.5)
}
## Obtain the Analysis of Variance directly and with aov function
SSTR <- sum(trt.n *(trt.mean-all.mean)^2); df_TR <- r-1
SSE <- sum((trt.n-1)*trt.sd^2); df_E <- n_T-r
SSTO <- SSTR + SSE; df_TO <- n_T-1
s2 <- SSE/df_E
aov.out1 <- rbind(df_TR,df_E,df_TO)
aov.out2 <- rbind(SSTR,SSE,SSTO)
aov.out3 <- rbind(SSTR/df_TR,SSE/df_E,NA)
aov.out4 <- rbind(aov.out3[1]/aov.out3[2],NA,NA)
aov.out5 <- rbind(qf(.95,df_TR,df_E),NA,NA)
aov.out6 <- rbind(1-pf(aov.out4[1],df_TR,df_E),NA,NA)
aov.out7 <- rbind(SSTR/SSTO,NA,NA)
aov.out <- cbind(aov.out1,aov.out2,aov.out3,aov.out4,aov.out5,aov.out6,aov.out7)
colnames(aov.out) <- c("df","Sum Sq", "Mean Sq", "F*", "F(.95)", "P(>F*)", "eta^2")
rownames(aov.out) <- c("Treatment", "Error", "Total")
round(aov.out,4)## df Sum Sq Mean Sq F* F(.95) P(>F*) eta^2
## Treatment 6 1256.654 209.4423 13.3707 2.1925 0 0.4501
## Error 98 1535.101 15.6643 NA NA NA NA
## Total 104 2791.755 NA NA NA NA NAoptions(contrasts=c("contr.sum", "contr.poly"))
mod1 <- aov(L.star ~ shipTrt.f)
summary.lm(mod1)##
## Call:
## aov(formula = L.star ~ shipTrt.f)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0693 -2.2446 0.1312 2.1744 16.2785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.7571 0.3862 61.508 < 2e-16 ***
## shipTrt.f1 -7.0571 0.9461 -7.459 3.55e-11 ***
## shipTrt.f2 -2.9571 0.9461 -3.126 0.002335 **
## shipTrt.f3 1.4428 0.9461 1.525 0.130468
## shipTrt.f4 0.7429 0.9461 0.785 0.434243
## shipTrt.f5 1.4429 0.9461 1.525 0.130463
## shipTrt.f6 3.4428 0.9461 3.639 0.000439 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.958 on 98 degrees of freedom
## Multiple R-squared: 0.4501, Adjusted R-squared: 0.4165
## F-statistic: 13.37 on 6 and 98 DF, p-value: 5.083e-11anova(mod1)## Analysis of Variance Table
##
## Response: L.star
## Df Sum Sq Mean Sq F value Pr(>F)
## shipTrt.f 6 1256.7 209.442 13.371 5.083e-11 ***
## Residuals 98 1535.1 15.664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1options(es.use_symbols=TRUE)
eta_squared(mod1)## For one-way between subjects designs, partial eta squared is equivalent
## to eta squared. Returning eta squared.## # Effect Size for ANOVA
##
## Parameter | η² | 95% CI
## -------------------------------
## shipTrt.f | 0.45 | [0.31, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].Contrast between Control and all other treatments c(6,-1,-1,-1,-1,-1,-1)
t_.975 <- qt(.975,df_E)
cont <- c(6, rep(-1,6))
sum(cont)## [1] 0Lhat <- sum(cont*trt.mean)
SE_Lhat <- sqrt(s2*sum(cont^2/trt.n))
SSC <- Lhat^2 / sum(cont^2/trt.n)
cont.out <- cbind(Lhat,SE_Lhat,Lhat/SE_Lhat,t_.975,2*(1-pt(abs(Lhat/SE_Lhat),df_E)),
Lhat-t_.975*SE_Lhat,Lhat+t_.975*SE_Lhat, SSC/s2, qf(.95,1,df_E),1-pf(SSC/s2,1,df_E))
colnames(cont.out) <- c("Lhat", "Std Err", "t_C", "t(.975)", "P(>|t_C|)", "LB", "UB", "F_C", "F(.95)", "P(>F_C)")
rownames(cont.out) <- c("Control vs Others")
round(cont.out,4)## Lhat Std Err t_C t(.975) P(>|t_C|) LB UB
## Control vs Others -49.4 6.6227 -7.4592 1.9845 0 -62.5425 -36.2574
## F_C F(.95) P(>F_C)
## Control vs Others 55.6396 3.9381 0Tukey Honest Significant Difference (Balanced n)
HSD <- qtukey(.95,r,n_T-r) * sqrt(s2/trt.n[1])
num.pairs <- r*(r-1)/2
p.count <- 0
tukey.out <- matrix(rep(0, num.pairs*9), ncol=9)
for (i1 in 1:(r-1)) {
for (i2 in (i1+1):r) {
p.count <- p.count + 1
tukey.out[p.count,1] <- i2
tukey.out[p.count,2] <- i1
tukey.out[p.count,3] <- trt.mean[i2]
tukey.out[p.count,4] <- trt.mean[i1]
tukey.out[p.count,5] <- trt.mean[i2]-trt.mean[i1]
tukey.out[p.count,6] <- HSD
tukey.out[p.count,7] <- tukey.out[p.count,5]-HSD
tukey.out[p.count,8] <- tukey.out[p.count,5]+HSD
tukey.out[p.count,9] <-
(1-ptukey(abs(tukey.out[p.count,5]/sqrt(s2/trt.n[1])),r,df_E))
}
}
colnames(tukey.out) <- c("Trt1","Trt2","Mean1","Mean2","Diff","HSD",
"Lower","Upper", "Adj P")
round(tukey.out, 6)## Trt1 Trt2 Mean1 Mean2 Diff HSD Lower Upper
## [1,] 2 1 20.80001 16.70001 4.100000 4.350696 -0.250696 8.450696
## [2,] 3 1 25.19999 16.70001 8.499980 4.350696 4.149284 12.850676
## [3,] 4 1 24.50000 16.70001 7.799993 4.350696 3.449297 12.150690
## [4,] 5 1 25.20001 16.70001 8.500000 4.350696 4.149304 12.850696
## [5,] 6 1 27.19999 16.70001 10.499980 4.350696 6.149284 14.850676
## [6,] 7 1 26.70001 16.70001 10.000000 4.350696 5.649304 14.350696
## [7,] 3 2 25.19999 20.80001 4.399980 4.350696 0.049284 8.750676
## [8,] 4 2 24.50000 20.80001 3.699993 4.350696 -0.650703 8.050690
## [9,] 5 2 25.20001 20.80001 4.400000 4.350696 0.049304 8.750696
## [10,] 6 2 27.19999 20.80001 6.399980 4.350696 2.049284 10.750676
## [11,] 7 2 26.70001 20.80001 5.900000 4.350696 1.549304 10.250696
## [12,] 4 3 24.50000 25.19999 -0.699987 4.350696 -5.050683 3.650710
## [13,] 5 3 25.20001 25.19999 0.000020 4.350696 -4.350676 4.350716
## [14,] 6 3 27.19999 25.19999 2.000000 4.350696 -2.350696 6.350696
## [15,] 7 3 26.70001 25.19999 1.500020 4.350696 -2.850676 5.850716
## [16,] 5 4 25.20001 24.50000 0.700007 4.350696 -3.650690 5.050703
## [17,] 6 4 27.19999 24.50000 2.699987 4.350696 -1.650710 7.050683
## [18,] 7 4 26.70001 24.50000 2.200007 4.350696 -2.150690 6.550703
## [19,] 6 5 27.19999 25.20001 1.999980 4.350696 -2.350716 6.350676
## [20,] 7 5 26.70001 25.20001 1.500000 4.350696 -2.850696 5.850696
## [21,] 7 6 26.70001 27.19999 -0.499980 4.350696 -4.850676 3.850716
## Adj P
## [1,] 0.078330
## [2,] 0.000001
## [3,] 0.000010
## [4,] 0.000001
## [5,] 0.000000
## [6,] 0.000000
## [7,] 0.045613
## [8,] 0.149958
## [9,] 0.045611
## [10,] 0.000481
## [11,] 0.001714
## [12,] 0.999001
## [13,] 1.000000
## [14,] 0.809083
## [15,] 0.943839
## [16,] 0.999001
## [17,] 0.506084
## [18,] 0.730839
## [19,] 0.809090
## [20,] 0.943842
## [21,] 0.999857TukeyHSD(mod1)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = L.star ~ shipTrt.f)
##
## $shipTrt.f
## diff lwr upr p adj
## 40mL-C 4.1000000 -0.2506965 8.450696 0.0783302
## 60mL-C 8.4999800 4.1492835 12.850676 0.0000012
## 80mL-C 7.7999933 3.4492968 12.150690 0.0000096
## 100mL-C 8.5000000 4.1493035 12.850696 0.0000012
## 120mL-C 10.4999800 6.1492835 14.850676 0.0000000
## 140mL-C 10.0000000 5.6493035 14.350696 0.0000000
## 60mL-40mL 4.3999800 0.0492835 8.750676 0.0456128
## 80mL-40mL 3.6999933 -0.6507032 8.050690 0.1499582
## 100mL-40mL 4.4000000 0.0493035 8.750696 0.0456111
## 120mL-40mL 6.3999800 2.0492835 10.750676 0.0004811
## 140mL-40mL 5.9000000 1.5493035 10.250696 0.0017137
## 80mL-60mL -0.6999867 -5.0506832 3.650710 0.9990013
## 100mL-60mL 0.0000200 -4.3506765 4.350716 1.0000000
## 120mL-60mL 2.0000000 -2.3506965 6.350696 0.8090827
## 140mL-60mL 1.5000200 -2.8506765 5.850716 0.9438387
## 100mL-80mL 0.7000067 -3.6506898 5.050703 0.9990012
## 120mL-80mL 2.6999867 -1.6507098 7.050683 0.5060844
## 140mL-80mL 2.2000067 -2.1506898 6.550703 0.7308392
## 120mL-100mL 1.9999800 -2.3507165 6.350676 0.8090899
## 140mL-100mL 1.5000000 -2.8506965 5.850696 0.9438422
## 140mL-120mL -0.4999800 -4.8506765 3.850716 0.9998567Bonferroni Simultaneous CI’s
num.pairs <- r*(r-1)/2
p.count <- 0
t.crit <- qt(1-.05/(2*num.pairs),n_T-r)
bonf.out <- matrix(rep(0, num.pairs*9), ncol=9)
for (i1 in 1:(r-1)) {
for (i2 in (i1+1):r) {
p.count <- p.count + 1
bonf.out[p.count,1] <- i2
bonf.out[p.count,2] <- i1
bonf.out[p.count,3] <- trt.mean[i2]
bonf.out[p.count,4] <- trt.mean[i1]
bonf.out[p.count,5] <- trt.mean[i2]-trt.mean[i1]
bonf.out[p.count,6] <- sqrt(s2*(1/trt.n[i2]+1/trt.n[i1]))
bonf.out[p.count,7] <- bonf.out[p.count,5]-t.crit*bonf.out[p.count,6]
bonf.out[p.count,8] <- bonf.out[p.count,5]+t.crit*bonf.out[p.count,6]
bonf.out[p.count,9] <- min(1,num.pairs*2*
(1-pt(abs(bonf.out[p.count,5]/bonf.out[p.count,6]),df_E)))
}
}
colnames(bonf.out) <- c("Trt1","Trt2","Mean1","Mean2","Diff","StdErr",
"Lower","Upper", "Adj P")
round(bonf.out, 4)## Trt1 Trt2 Mean1 Mean2 Diff StdErr Lower Upper Adj P
## [1,] 2 1 20.8 16.7 4.1 1.4452 -0.4080 8.6080 0.1162
## [2,] 3 1 25.2 16.7 8.5 1.4452 3.9920 13.0079 0.0000
## [3,] 4 1 24.5 16.7 7.8 1.4452 3.2920 12.3080 0.0000
## [4,] 5 1 25.2 16.7 8.5 1.4452 3.9920 13.0080 0.0000
## [5,] 6 1 27.2 16.7 10.5 1.4452 5.9920 15.0079 0.0000
## [6,] 7 1 26.7 16.7 10.0 1.4452 5.4920 14.5080 0.0000
## [7,] 3 2 25.2 20.8 4.4 1.4452 -0.1080 8.9079 0.0628
## [8,] 4 2 24.5 20.8 3.7 1.4452 -0.8080 8.2080 0.2517
## [9,] 5 2 25.2 20.8 4.4 1.4452 -0.1080 8.9080 0.0628
## [10,] 6 2 27.2 20.8 6.4 1.4452 1.8920 10.9079 0.0005
## [11,] 7 2 26.7 20.8 5.9 1.4452 1.3920 10.4080 0.0019
## [12,] 4 3 24.5 25.2 -0.7 1.4452 -5.2080 3.8080 1.0000
## [13,] 5 3 25.2 25.2 0.0 1.4452 -4.5079 4.5080 1.0000
## [14,] 6 3 27.2 25.2 2.0 1.4452 -2.5080 6.5080 1.0000
## [15,] 7 3 26.7 25.2 1.5 1.4452 -3.0079 6.0080 1.0000
## [16,] 5 4 25.2 24.5 0.7 1.4452 -3.8080 5.2080 1.0000
## [17,] 6 4 27.2 24.5 2.7 1.4452 -1.8080 7.2080 1.0000
## [18,] 7 4 26.7 24.5 2.2 1.4452 -2.3080 6.7080 1.0000
## [19,] 6 5 27.2 25.2 2.0 1.4452 -2.5080 6.5079 1.0000
## [20,] 7 5 26.7 25.2 1.5 1.4452 -3.0080 6.0080 1.0000
## [21,] 7 6 26.7 27.2 -0.5 1.4452 -5.0079 4.0080 1.0000pairwise.t.test(L.star, shipTrt.f, p.adjust.method="bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: L.star and shipTrt.f
##
## C 40mL 60mL 80mL 100mL 120mL
## 40mL 0.11622 - - - - -
## 60mL 1.2e-06 0.06284 - - - -
## 80mL 9.9e-06 0.25174 1.00000 - - -
## 100mL 1.2e-06 0.06284 1.00000 1.00000 - -
## 120mL 1.9e-09 0.00052 1.00000 1.00000 1.00000 -
## 140mL 9.9e-09 0.00191 1.00000 1.00000 1.00000 1.00000
##
## P value adjustment method: bonferroniSNK and Duncan’s Multiple Range Tests
(SNK.crit <- qtukey(.95,2:7,df_E) * sqrt(s2/trt.n[1]))## [1] 2.867932 3.439325 3.777241 4.016499 4.200997 4.350696(Duncan.crit <- qtukey(.95^(1:6),2:7,df_E) * sqrt(s2/trt.n[1]))## [1] 2.867932 3.018065 3.117729 3.190725 3.247375 3.293059SNK.test(mod1, "shipTrt.f", console=TRUE)##
## Study: mod1 ~ "shipTrt.f"
##
## Student Newman Keuls Test
## for L.star
##
## Mean Square Error: 15.6643
##
## shipTrt.f, means
##
## L.star std r se Min Max Q25 Q50 Q75
## 100mL 25.20001 3.098391 15 1.021903 20.3186 32.8077 23.38860 24.4906 27.31270
## 120mL 27.19999 3.485686 15 1.021903 19.7150 31.9919 25.40440 28.6501 29.26755
## 140mL 26.70001 2.711081 15 1.021903 21.6483 29.9627 24.87850 26.8312 28.91415
## 40mL 20.80001 6.971377 15 1.021903 10.7307 37.0785 16.24455 21.0834 24.60190
## 60mL 25.19999 3.872979 15 1.021903 17.9319 32.2451 23.24740 25.4399 27.65170
## 80mL 24.50000 2.711088 15 1.021903 20.4410 30.6687 22.91485 24.1533 25.11805
## C 16.70001 3.098390 15 1.021903 10.1485 21.5381 14.37980 17.7575 18.31885
##
## Alpha: 0.05 ; DF Error: 98
##
## Critical Range
## 2 3 4 5 6 7
## 2.867932 3.439325 3.777241 4.016499 4.200997 4.350696
##
## Means with the same letter are not significantly different.
##
## L.star groups
## 120mL 27.19999 a
## 140mL 26.70001 a
## 100mL 25.20001 a
## 60mL 25.19999 a
## 80mL 24.50000 a
## 40mL 20.80001 b
## C 16.70001 cduncan.test(mod1, "shipTrt.f", console=TRUE)##
## Study: mod1 ~ "shipTrt.f"
##
## Duncan's new multiple range test
## for L.star
##
## Mean Square Error: 15.6643
##
## shipTrt.f, means
##
## L.star std r se Min Max Q25 Q50 Q75
## 100mL 25.20001 3.098391 15 1.021903 20.3186 32.8077 23.38860 24.4906 27.31270
## 120mL 27.19999 3.485686 15 1.021903 19.7150 31.9919 25.40440 28.6501 29.26755
## 140mL 26.70001 2.711081 15 1.021903 21.6483 29.9627 24.87850 26.8312 28.91415
## 40mL 20.80001 6.971377 15 1.021903 10.7307 37.0785 16.24455 21.0834 24.60190
## 60mL 25.19999 3.872979 15 1.021903 17.9319 32.2451 23.24740 25.4399 27.65170
## 80mL 24.50000 2.711088 15 1.021903 20.4410 30.6687 22.91485 24.1533 25.11805
## C 16.70001 3.098390 15 1.021903 10.1485 21.5381 14.37980 17.7575 18.31885
##
## Alpha: 0.05 ; DF Error: 98
##
## Critical Range
## 2 3 4 5 6 7
## 2.867932 3.018065 3.117729 3.190725 3.247375 3.293059
##
## Means with the same letter are not significantly different.
##
## L.star groups
## 120mL 27.19999 a
## 140mL 26.70001 a
## 100mL 25.20001 a
## 60mL 25.19999 a
## 80mL 24.50000 a
## 40mL 20.80001 b
## C 16.70001 cNo differences compared to what original data found. We agree with the authors original findings.
Bartlett’s Test for equal variances
Welch’s F-test
Kruskal-Wallis Test
bartlett.test(L.star ~ shipTrt.f)##
## Bartlett test of homogeneity of variances
##
## data: L.star by shipTrt.f
## Bartlett's K-squared = 22.594, df = 6, p-value = 0.0009446oneway.test(L.star ~ shipTrt.f, var.equal=FALSE)##
## One-way analysis of means (not assuming equal variances)
##
## data: L.star and shipTrt.f
## F = 18.264, num df = 6.000, denom df = 43.356, p-value = 1.934e-10kruskal.test(L.star ~ shipTrt.f)##
## Kruskal-Wallis rank sum test
##
## data: L.star by shipTrt.f
## Kruskal-Wallis chi-squared = 45.252, df = 6, p-value = 4.17e-08