ARM Analysis : assessment1
What follows are code fragments used to generate ARM reports. Code to redirect R messages and to write results to tables are not included.
This document has been modified to from Save to R Markdown to hide code used for diagnostics.
Data
arm.dat <- data.frame(
plot=as.factor(c(101, 102, 103, 104, 105, 106, 107, 201, 202, 203, 204, 205, 206, 207, 301, 302, 303, 304, 305, 306, 307, 401, 402, 403, 404, 405, 406, 407, 501, 502, 503, 504, 505, 506, 507, 601, 602, 603, 604, 605, 606, 607, 701, 702, 703, 704, 705, 706, 707, 801, 802, 803, 804, 805, 806, 807, 901, 902, 903, 904, 905, 906, 907, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1201, 1202, 1203, 1204, 1205, 1206,
1207, 1301, 1302, 1303, 1304, 1305, 1306, 1307, 1401, 1402, 1403, 1404, 1405, 1406, 1407)),
treatment=as.factor(c(1, 2, 3, 4, 5, 6, 7, 6, 7, 4, 3, 1, 2, 5, 1, 4, 3, 2, 6, 5, 7, 7, 5, 4, 6, 1, 3, 2, 1, 7, 6, 4, 2, 5, 3, 4, 1, 7, 5, 3, 6, 2, 3, 4, 2, 6, 1, 5, 7, 1, 3, 5, 4, 6, 7, 2, 4, 2, 3, 1, 7, 6, 5, 6, 7, 1, 5, 4, 2, 3, 2, 1, 4, 6, 3, 7, 5, 4, 6, 2, 1, 5, 3, 7, 3, 7, 5, 4, 6, 1, 2, 2, 3, 1, 7, 4, 6, 5)),
replicate=as.factor(c(1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14)),
subsample =as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),
block=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),
column=as.factor(c(1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7)),
blkcol=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),
assessment1=c(103.30000000, 68.42000000, 89.53000000, 77.94000000, 100.00000000, 108.20000000, 184.90000000, 61.74000000, 102.30000000, 73.90000000, 60.91000000, 42.19000000, 44.31000000, 43.87000000, 71.27000000, 86.52000000, 100.71000000, 81.88000000, 90.75000000, 100.30000000, 100.60000000, 117.10000000, 71.82000000, 60.79000000, 83.04000000, 37.59000000, 56.97000000, 60.05000000, 58.31000000, 85.77000000, 89.78000000, 60.42000000, 52.73000000,
104.33000000, 96.04000000, 181.00000000, 46.27000000, 246.80000000, 172.33000000, 129.80000000, 164.80000000, 82.18000000, 38.70000000, 45.81000000, 23.07000000, 50.69000000, 19.00000000, 59.03000000, 55.18000000, 8.31000000, 9.16000000, 10.63000000, 14.28000000, 15.84000000, 12.50000000, 8.43000000, 69.35000000, 33.48000000, 60.08000000, 20.15000000, 101.60000000, 45.28000000, 114.34000000, 41.52000000, 66.33000000, 3.16000000, 8.89000000, 18.95000000,
4.12000000, 12.73000000, 18.87000000, 4.12000000, 12.92000000, 19.85000000, 8.54000000, 16.76000000, 25.30000000, 10.95000000, 13.96000000, 11.18000000, 7.68000000, 10.54000000, 10.44000000, 14.39000000, 11.22000000, 69.38000000, 42.25000000, 25.36000000, 48.80000000, 5.29000000, 10.91000000, 5.57000000, 22.23000000, 6.00000000, 4.00000000, 2.45000000, 1.00000000, 6.24000000)
)Design Model
The standard ARM analysis is based on linear model fit using lm. This model is only valid if there are a sufficient number of assessments
Diagnostics
Before proceeding with the standard analysis, we check model assumptions.
## plot treatment replicate column assessment1 StdRes
## 37 602 1 6 2 46.27 -3.62831## Raw IID EX
## Mean 53.75786 0.0000000 0.000000
## Missing 0.00000 0.0000000 1.000000
## N 98.00000 98.0000000 97.000000
## ErrorDF 97.00000 78.0000000 77.000000
## SD 47.45880 21.2340112 18.970655
## SD_reml 47.45880 19.0411633 16.989949
## SD_ml 47.21605 18.9437654 16.902145
## g_1 1.33257 -0.1184116 0.153632
## G_1 1.35337 -0.1202602 0.156055
## SES 0.24747 0.2474726 0.248746
## SkewnessT 5.46877 -0.4859534 0.627369
## Skewness 0.00000 0.6280962 0.531907
## g_2 2.26408 3.5832181 2.184234
## G_2 2.44781 3.8368030 2.365652
## SEK 0.49032 0.4903206 0.492800
## KurtosisT 4.99225 7.8250898 4.800435
## Kurtosis 0.00000 0.0000000 0.000000
## MinRep 14.00000 14.0000000 13.000000
## MaxRep 14.00000 14.0000000 14.000000
## Unique 97.00000 98.0000000 97.000000
## Treatments 7.00000 7.0000000 7.000000
## ShapiroWilksW 0.87351 0.9422911 0.957409
## ShapiroWilks 0.00000 0.0003106 0.003180
## LeveneF 1.44304 2.9223896 4.093056
## Levenes 0.20696 0.0118554 0.001121
## BartlettK 13.45753 24.9023972 23.438979
## Bartlett 0.03632 0.0003560 0.000662Recommendations
This section considers where alternative analysis may be preferred over the standard analysis.
## ActionCode Criteria Value
## Levenes IID Levenes NA
## ShapiroWilks IID ShapiroWilks NA
## Skewness IID Skewness NA
## Kurtosis IID Kurtosis NA
## Comment
## Levenes Transformation specified, skipping recommendations
## ShapiroWilks Transformation specified, skipping recommendations
## Skewness Transformation specified, skipping recommendations
## Kurtosis Transformation specified, skipping recommendationsStandard Analysis
The standard analysis will be an analysis of variance of the linear model determined by the design of the experiment.
Alternate Analysis
Alternative analysis suggested by the Recommendations may be included here.
Plot
Analysis of Ranks
This chunk of code is taken from the R script written during Print Reports
arm.dat <- subset(arm.dat,!is.na(arm.dat$assessment1))
UniqueTreatments <- sort(unique(arm.dat$treatment))
UniqueReplicates <- unique(arm.dat$replicate)
for(rep in UniqueReplicates) {
mask <- arm.dat$replicate == rep
arm.dat$rank[mask] <- rank(arm.dat$assessment1[mask], ties.method='average')
}
b <- length(UniqueReplicates)
k <- length(UniqueTreatments)
s <- tapply(arm.dat$rank, list(arm.dat$treatment), sum)[UniqueTreatments]
S_r <- sum(arm.dat$rank^2)
S_k <- sum(s^2)/b
C <- (b*k*(k+1)^2)/4
T <- (b*(k-1)*(S_k - C))/(S_r-C)
print(RankChi2 <- T)## [1] 40.10204print(RankDF <- k-1)## [1] 6print(RankP <- 1-pchisq(T,df=k-1))## [1] 4.349584e-07print(CritChi <- qchisq(1 - alpha, df=k-1))## [1] 12.59159print(RankSum <- tapply(arm.dat$rank, list(arm.dat$treatment), sum)[UniqueTreatments])## 1 2 3 4 5 6 7
## 23 36 53 57 70 70 83se_R <- sqrt(2*b*(S_r - S_k)/((b-1)*(k-1)))
err.df <- (b-1)*(k-1)
crit.t <-qt(1 - 0.05 / 2, err.df)
Rmin <- se_R*crit.t
crit.table <- matrix(Rmin,nrow=k,ncol=k)
for (i in 1:k) {
crit.table[i,i] <- 0
}
post.tbl <- matrix(0,nrow=k,ncol=k)
for(i in 1:(k-1)) {
for(j in (i+1):k) {
post.tbl[i,j] <- abs(RankSum[i]-RankSum[j])
post.tbl[j,i] <- post.tbl[i,j]>crit.table[i,j]
crit.table[j,i] <- 2*(1-pt(post.tbl[i,j]/se_R,err.df))
}
}
print(post.tbl)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 13 30 34 47 47 60
## [2,] 0 0 17 21 34 34 47
## [3,] 1 0 0 4 17 17 30
## [4,] 1 1 0 0 13 13 26
## [5,] 1 1 0 0 0 0 13
## [6,] 1 1 0 0 0 0 13
## [7,] 1 1 1 1 0 0 0print(crit.table)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000e+00 1.707242e+01 1.707242e+01 17.072423291 17.0724233
## [2,] 1.335733e-01 0.000000e+00 1.707242e+01 17.072423291 17.0724233
## [3,] 7.765195e-04 5.095592e-02 0.000000e+00 17.072423291 17.0724233
## [4,] 1.617159e-04 1.657505e-02 6.421954e-01 0.000000000 17.0724233
## [5,] 5.028270e-07 1.617159e-04 5.095592e-02 0.133573283 0.0000000
## [6,] 5.028270e-07 1.617159e-04 5.095592e-02 0.133573283 1.0000000
## [7,] 7.991856e-10 5.028270e-07 7.765195e-04 0.003297825 0.1335733
## [,6] [,7]
## [1,] 17.0724233 17.07242
## [2,] 17.0724233 17.07242
## [3,] 17.0724233 17.07242
## [4,] 17.0724233 17.07242
## [5,] 17.0724233 17.07242
## [6,] 0.0000000 17.07242
## [7,] 0.1335733 0.00000Analysis in R
These code chunks are not included in ARM scripts.
In base R, we have
friedman.test(assessment1 ~ treatment | replicate, data = arm.dat)##
## Friedman rank sum test
##
## data: assessment1 and treatment and replicate
## Friedman chi-squared = 40.102, df = 6, p-value = 4.35e-07agricolae also provides Friedman’s test and mean comparisons
library(agricolae)## Warning: package 'agricolae' was built under R version 3.6.1res <- friedman(arm.dat$replicate, arm.dat$treatment, arm.dat$assessment1)
print(res)## $statistics
## Chisq Df p.chisq F DFerror p.F t.value LSD
## 40.10204 6 4.349584e-07 11.87587 78 2.012652e-09 1.990847 17.07242
##
## $parameters
## test name.t ntr alpha
## Friedman arm.dat$treatment 7 0.05
##
## $means
## arm.dat.assessment1 rankSum std r Min Max Q25 Q50
## 1 30.90286 23 30.50910 14 3.16 103.30 6.4200 19.575
## 2 36.08571 36 28.59641 14 4.12 82.18 10.9775 28.275
## 3 50.50429 53 40.79134 14 8.54 129.80 11.5975 47.835
## 4 52.90286 57 46.70163 14 2.45 181.00 15.4475 53.115
## 5 62.13357 70 50.38499 14 6.24 172.33 14.2975 51.450
## 6 59.66071 70 44.37954 14 1.00 164.80 25.2675 49.745
## 7 84.11500 83 68.46616 14 4.00 246.80 26.3650 77.575
## Q75
## 1 45.2500
## 2 58.2200
## 3 82.3750
## 4 72.7625
## 5 100.2250
## 6 88.0950
## 7 102.1250
##
## $comparison
## NULL
##
## $groups
## Sum of ranks groups
## 7 83 a
## 5 70 ab
## 6 70 ab
## 4 57 b
## 3 53 bc
## 2 36 cd
## 1 23 d
##
## attr(,"class")
## [1] "group"