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, 201, 202, 203, 301, 302, 303, 401, 402, 403, 501, 502, 503)), 
treatment=as.factor(c(1, 3, 2, 3, 2, 1, 2, 1, 3, 3, 2, 1, 2, 1, 3)), 
replicate=as.factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5)), 
subsample =as.factor(c(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)), 
column=as.factor(c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3)), 
blkcol=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), 
assessment1=c(340.00000000, 347.00000000, 339.00000000, 343.00000000, 333.00000000, 345.00000000, 344.00000000, 330.00000000, 349.00000000, 355.00000000, NA, 342.00000000, NA, 338.00000000, NA)
)

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.

##                      Raw       IID
## Mean          342.083333  0.000000
## Missing         3.000000  3.000000
## N              12.000000 12.000000
## ErrorDF        11.000000  9.000000
## SD              6.815201  5.412605
## SD_reml         6.815201  4.895886
## SD_ml           6.525059  4.687454
## g_1             0.005179 -0.332724
## G_1             0.005950 -0.382271
## SES             0.714598  0.714598
## SkewnessT       0.008326 -0.534945
## Skewness        0.993506  0.603333
## g_2            -0.273572 -0.844859
## G_2             0.298658 -0.609053
## SEK             1.381702  1.381702
## KurtosisT       0.216152 -0.440799
## Kurtosis        0.832825  0.667896
## MinRep          3.000000  3.000000
## MaxRep          5.000000  5.000000
## Unique         12.000000 12.000000
## Treatments      3.000000  3.000000
## ShapiroWilksW   0.988326  0.948111
## ShapiroWilks    0.999275  0.609560
## LeveneF         0.008297  0.008297
## Levenes         0.991744  0.991744
## BartlettK       0.045682  0.045682
## Bartlett        0.977418  0.977418

Recommendations

This section considers where alternative analysis may be preferred over the standard analysis.

##              ActionCode     Criteria     Value
## Levenes             IID      Levenes  0.008297
## ShapiroWilks        IID ShapiroWilks  0.948111
## Skewness            IID     Skewness -0.534945
## Kurtosis            IID     Kurtosis -0.440799
##                                                                      Comment
## Levenes      Does not fail test of homogeneity of variances among treatments
## ShapiroWilks            Does not fail general test of normality of residuals
## Skewness                         Does not fail test of skewness of residuals
## Kurtosis                  Does not fail test of excess kurtosis of residuals

Standard Analysis

The standard analysis will be an analysis of variance of the linear model determined by the design of the experiment.

#skipping standard analysis in diagnostic script

Alternate Analysis

Alternative analysis suggested by the Recommendations may be included here.

#skipping alternate analysis in diagnostic script

Plot

Analysis of Ranks

This chunk of code is taken from the R script written during Print Reports

arm.dat <- arm.dat[!is.na(arm.dat$assessment1),]
arm.dat$rank <- rank(arm.dat$assessment1, ties.method='average')
R <- tapply(arm.dat$rank, list(arm.dat$treatment), sum)
n <- tapply(arm.dat$rank, list(arm.dat$treatment), length)
k <- length(R)
N <- sum(n)
S2 <- (1/(N-1))*(sum(arm.dat$rank^2) - (N*(N+1)^2)/4)
T <- (1/S2)*(sum(R^2/n)-(N*(N+1)^2)/4)

print(RankChi2 <- T)
## [1] 5.65641
print(RankDF <- k-1)
## [1] 2
print(RankP <- 1-pchisq(T,df=k-1))
## [1] 0.05911887
print(CritChi <- qchisq(1 - alpha, df=k-1))
## [1] 5.991465
print(means <- R/n)
##         1         2         3 
##  4.800000  4.666667 10.000000
sd.rank <- sqrt(S2*((N-1-T)/(N-k)))
crit.table <- matrix(0,nrow=k,ncol=k)
post.tbl <- matrix(0,nrow=k,ncol=k)
crit.t <- qt(1-alpha/2,N-k)
for(i in 1:(k-1)) {
  for(j in (i+1):k) {
    post.tbl[i,j] <- abs(means[i]-means[j])
    crit.table [i,j] <- crit.t*sd.rank*(sqrt(1/n[i]+1/n[j]))
    post.tbl[j,i] <- post.tbl[i,j]>crit.table[i,j]
    crit.table [j,i] <- 2*(1-pt(post.tbl[i,j]/(sd.rank*(sqrt(1/n[i]+1/n[j]))),N-k))
  }
}
print(crit.table)
##            [,1]       [,2]     [,3]
## [1,] 0.00000000 4.58975088 4.215955
## [2,] 0.94904040 0.00000000 4.800076
## [3,] 0.02104812 0.03312175 0.000000
print(post.tbl)
##      [,1]      [,2]     [,3]
## [1,]    0 0.1333333 5.200000
## [2,]    0 0.0000000 5.333333
## [3,]    1 1.0000000 0.000000

Analysis in R

These code chunks are not included in ARM scripts.

In base R, we have

kruskal.test(assessment1 ~ treatment, data=arm.dat)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  assessment1 by treatment
## Kruskal-Wallis chi-squared = 5.6564, df = 2, p-value = 0.05912

agricolae also provides Kruskal-Wallis and mean comparisons

library(agricolae)
## Warning: package 'agricolae' was built under R version 3.6.1
res <- kruskal(arm.dat$assessment1, arm.dat$treatment, p.adj="none")
print(res)
## $statistics
##     Chisq Df    p.chisq
##   5.65641  2 0.05911887
## 
## $parameters
##             test p.ajusted            name.t ntr alpha
##   Kruskal-Wallis      none arm.dat$treatment   3  0.05
## 
## $means
##   arm.dat.assessment1      rank      std r Min Max Q25 Q50   Q75
## 1            339.0000  4.800000 5.656854 5 330 345 338 340 342.0
## 2            338.6667  4.666667 5.507571 3 333 344 336 339 341.5
## 3            348.5000 10.000000 5.000000 4 343 355 346 348 350.5
## 
## $comparison
## NULL
## 
## $groups
##   arm.dat$assessment1 groups
## 3           10.000000      a
## 1            4.800000      b
## 2            4.666667      b
## 
## attr(,"class")
## [1] "group"

```