---
title: 'ARM Analysis : assessment1'
output: html_document
---
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.

```{r,echo=FALSE}
# Utility Functions
moments <- function(x, df=0, alpha=0.05) {
  missing <- sum(is.na(x))
  x <- x[!is.na(x)]
  bar_x <- mean(x)
  N <- length(x)
  if(df==0) {
    df <- N-1
  }
  SS <- sum((x- bar_x)^2)
  s_reml <- sqrt(SS/(N-1))
  s_ml <- sqrt(SS/N)
  g_1 <- (sum((x - bar_x)^3)/N)/(s_ml^3)
  G_1 <- g_1*sqrt(N*(N-1))/(N-2)
  SES <- sqrt((6*N*N-1)/((N-2)*(N+1)*(N+1)))
  skew.t <- G_1/SES
  skew.p <- 2*(1-pt(abs(skew.t),N-1))
  g_2 <- (sum((x - bar_x)^4)/N)/(s_ml^4)-3
  G_2 <- (sum((x - bar_x)^4))/(s_reml^4) * (N*(N+1))/((N-1)*(N-2)*(N-3))
  G_2 <- G_2 - (3*(N-1)^2)/((N-2)*(N-3))
  SEK <- 2*SES*(sqrt((N^2-1)/((N-3)*(N+5))))
  kurt.t <- G_2/SEK
  kurt.p <- 2*(1-pt(abs(kurt.t),N-1))
  return(list(Mean=bar_x,
    Missing = missing,
    N = N,
    ErrorDF = df,
    SD = sqrt(SS/df),
    SD_reml = s_reml,
    SD_ml = s_ml,
    g_1 = g_1,
    G_1 = G_1,
    SES = SES,
    SkewnessT = skew.t,
    Skewness = ifelse(skew.p<1e-4,0,skew.p),
    g_2 = g_2,
    G_2 = G_2,
    SEK = SEK,
    KurtosisT = kurt.t,
    Kurtosis = ifelse(kurt.p<1e-4,0,kurt.p)))
}
normality <- function(x, g, alpha=0.05,dist='pnorm') {
  g <- g[!is.na(x)]
  x <- x[!is.na(x)]
   tbl <- tapply(x,list(g),length)
  ret <- list(MinRep = min(tbl,na.rm=TRUE), MaxRep = max(tbl,na.rm=TRUE), Unique=length(unique(x)), Treatments = length(unique(g)),
       ShapiroWilksW=NA,
      ShapiroWilks=NA,
      LeveneF=NA,
      Levenes=NA,
      BartlettK=NA,
      Bartlett=NA)
  if(ret$Unique>1) { 
    sw <- shapiro.test(x)
    names(sw$statistic) <- NULL
    ret$ShapiroWilksW =sw$statistic
    ret$ShapiroWilks=ifelse(sw$p.value<1e-4,0,sw$p.value)
  }
  if(ret$MaxRep>2 & ret$Treatments>1) { 
    z.med <- tapply(x,g,median, na.rm=TRUE)
    z.dat <- z.med[g]
    lev.tbl  <- summary(aov(abs(z.dat-x) ~ g))[[1]]
    ret$LeveneF = lev.tbl[1,4]
    ret$Levenes = ifelse(lev.tbl[1,5]<1e-4,0,lev.tbl[1,5])
  }
  if(ret$MinRep>2 & ret$Treatments>1) { 
    bar <- bartlett.test(x=x, g=g)
    names(bar$statistic) <- NULL
    ret$BartlettK=ifelse(is.finite(bar$statistic),bar$statistic,NA)
    ret$Bartlett=ifelse(bar$p.value<1e-4,0,bar$p.value)
  }
  return(ret)
}
```
```{r,echo=FALSE}
# Analysis Constants
alpha <- 0.05000000
```
# Data
```{r}
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
```{r,echo=FALSE}
raw.moments <- moments(arm.dat$assessment1)
raw.normality <- normality(arm.dat$assessment1, arm.dat$treatment)
if(raw.normality$MaxRep>1 & raw.normality$Treatments>1) {
analysis.model <- lm(assessment1 ~ treatment, data = arm.dat, x = TRUE)
}
```
# Diagnostics
Before proceeding with the standard analysis, we check model assumptions.
```{r,echo=FALSE}
Criteria <- c('Levenes','ShapiroWilks','Skewness','Kurtosis')
TestColumn <- 'IID'
if(!exists('analysis.model')) {
  if(raw.normality$MaxRep==1) {
    Diagnostics <- data.frame(Raw = c(unlist(raw.moments), unlist(raw.normality)))
    Diagnostics$N <- NA
    TestColumn <- 'N'
  } else {
    Diagnostics <- data.frame(Raw = c(unlist(raw.moments), unlist(raw.normality)))
    Diagnostics$ES <- NA
    TestColumn <- 'ES'
  }
} else {
  Residuals <- data.frame(Plot = arm.dat$plot,Treatment=arm.dat$treatment)
  Residuals$IID <- resid(analysis.model)[row.names(arm.dat)]
  lm.moments <- moments(Residuals$IID, df=analysis.model$df.residual)
  lm.normality <- normality(Residuals$IID, Residuals$Treatment)
  models <- list(IID=analysis.model)
  Diagnostics <- data.frame(Raw = c(unlist(raw.moments), unlist(raw.normality)),
                            IID = c(unlist(lm.moments), unlist(lm.normality)))
  best.col<-2
  Diagnostics[1,2:dim(Diagnostics)[2]] <- 0
  normality.tests <- Diagnostics[Criteria,TestColumn]
  normality.tests[is.na(normality.tests)] <- 1
  if(Diagnostics['SD','IID']<1e-6) { 
    Diagnostics[Criteria,2] <- NA
    names(Diagnostics)[2] <- 'NM'
    TestColumn <- 'NM'
  } else {
    arm.dat$StdRes <- ifelse(is.na(Residuals$IID),0,Residuals$IID/lm.moments$SD)
    qoutlier <- qnorm(1-(alpha/(2*length(Residuals$IID))))
    arm.dat$Outlier <- abs(arm.dat$StdRes)>qoutlier
    Outliers <- arm.dat[arm.dat$Outlier,c('plot','treatment','replicate','column','assessment1','StdRes')]
    if(dim(Outliers)[1]>0) {
      print(Outliers)
      outlier.lm <- lm(assessment1 ~ treatment + Outlier, data = arm.dat)
      ex3.lm <- lm(assessment1 ~ treatment + Outlier, data = arm.dat[!arm.dat$Outlier,])
      Residuals$EX  <- resid(ex3.lm)[row.names(arm.dat)]
      Diagnostics$EX = c(unlist(moments(Residuals$EX,ex3.lm$df.residual)), unlist(normality(Residuals$EX,Residuals$Treatment)))
      models <- append(models, list(EX=ex3.lm))
    }
    if(any(normality.tests < 0.05)) { 
        if( (min(arm.dat$assessment1,na.rm=TRUE)>=0) ) {
        if( (min(arm.dat$assessment1,na.rm=TRUE)>0) ) {
          arm.dat$logx <- log(arm.dat$assessment1)
        } else {
        arm.dat$logx <- log(.Machine$double.eps + arm.dat$assessment1)
      }
      log.lm <- lm(logx ~ treatment, data = arm.dat)
      Residuals$AL  <- resid(log.lm)[row.names(arm.dat)]
      Diagnostics$AL = c(unlist(moments(Residuals$AL,log.lm$df.residual)), unlist(normality(Residuals$AL,Residuals$Treatment)))
      models <- append(models, list(AL=log.lm))
      arm.dat$sqrtx <- sqrt(0.5 + arm.dat$assessment1)
      sqrt.lm <- lm(sqrtx ~ treatment, data = arm.dat)
      Residuals$AS  <- resid(sqrt.lm)[row.names(arm.dat)]
      Diagnostics$AS = c(unlist(moments(Residuals$AS,sqrt.lm$df.residual)), unlist(normality(Residuals$AS,Residuals$Treatment)))
      models <- append(models, list(AS=sqrt.lm))
      if(max(arm.dat$assessment1,na.rm=TRUE)<=100) {
        if( max(arm.dat$assessment1,na.rm=TRUE)<=1) {
          arm.dat$asinx <- asin(sqrt((arm.dat$assessment1)))
          asin.lm <- lm(asinx ~ treatment, data = arm.dat)
        } else {
          arm.dat$asinx <- asin(sqrt((arm.dat$assessment1)/100))
          asin.lm <- lm(asinx ~ treatment, data = arm.dat)
        }
        Residuals$AA  <- resid(asin.lm)[row.names(arm.dat)]
        Diagnostics$AA = c(unlist(moments(Residuals$AA)), unlist(normality(Residuals$AA,Residuals$Treatment)))
        models <- append(models, list(AA=asin.lm))
      }
      }
    block.lm <- lm(assessment1 ~ 0, data = arm.dat)
    arm.dat$AlignedRank  <- rank(resid(block.lm))[row.names(arm.dat)]
    arank.lm <- lm(AlignedRank ~ treatment, data = arm.dat)
    models <- append(models, list(AR=arank.lm))
    Residuals$AR  <- resid(arank.lm)[row.names(arm.dat)]
    Diagnostics$AR <- c(unlist(moments(Residuals$AR,arank.lm$df.residual)), unlist(normality(Residuals$AR,Residuals$Treatment)))
    rank.col <- dim(Diagnostics)[2]
    exclude.cols <-c(1,rank.col)
    }
  }
}
Diagnostics[1,2:dim(Diagnostics)[2]] <- 0
format(Diagnostics, digits=4,scientific=FALSE)
```
# Recommendations
This section considers where alternative analysis may be preferred over the standard analysis.
```{r,echo=FALSE}
Recommendations <- data.frame(ActionCode = 'IID',Criteria=Criteria,Value=rep(NA,4),Comment = rep(' ',4), stringsAsFactors=FALSE)
row.names(Recommendations) <- Criteria
Recommendations['Levenes','Comment'] <- 'Does not fail test of homogeneity of variances among treatments'
Recommendations['ShapiroWilks','Comment'] <- 'Does not fail general test of normality of residuals'
Recommendations['Skewness','Comment'] <- 'Does not fail test of skewness of residuals'
Recommendations['Kurtosis','Comment'] <- 'Does not fail test of excess kurtosis of residuals'
if(exists('N',where=Diagnostics)) { 
  Recommendations$ActionCode = 'N'
  Recommendations$Comment = 'Unreplicated, mark as non-analyzable'
} else if(Diagnostics['SD',TestColumn]<1e-6) { 
  Recommendations$ActionCode = 'NM'
  Recommendations$Comment = 'Non-analyzable, error is effectively 0'
} else {
  Recommendations['Levenes','Value'] <- Diagnostics[c('LeveneF'),'IID']
  Recommendations['ShapiroWilks','Value'] <- Diagnostics[c('ShapiroWilksW'),'IID']
  Recommendations['Skewness','Value'] <- Diagnostics[c('SkewnessT'),'IID']
  Recommendations['Kurtosis','Value'] <- Diagnostics[c('KurtosisT'),'IID']
  if(!is.na(Diagnostics[c('Levenes'),2])) { 
    if(Diagnostics[c('Levenes'),2] < 0.05) { 
      best.col <- which.min(abs(Diagnostics[c('LeveneF'),-exclude.cols]))+1 
      if(Diagnostics['Levenes',best.col] < 0.05 | !(best.col>2)) {
        Recommendations['Levenes','Comment'] = 'Homogeneity of variances not stabilized by available transformations'
        Recommendations['Levenes','ActionCode'] = 'AR'
      } else {
        Recommendations['Levenes','ActionCode'] = names(Diagnostics)[best.col]
        if(Recommendations['Levenes','ActionCode'] == 'EX') {
          Recommendations['Levenes','Comment'] = 'Exclude outliers to improve homogeneity of variance'
        } else {
          Recommendations['Levenes','Comment'] = 'Transform to stabilize variance'
        }
      }
    }
  }
  if(Diagnostics['ShapiroWilks','IID'] < 0.05) { 
    best.col <- which.max(Diagnostics['ShapiroWilksW',-exclude.cols])+1 
    if(Diagnostics[c('ShapiroWilks'),best.col] < 0.05 | !(best.col>2)) {
      Recommendations['ShapiroWilks','ActionCode'] = 'AR'
      Recommendations['ShapiroWilks','Comment'] = 'Available transformations do not improve normality'
    } else {
      Recommendations['ShapiroWilks','ActionCode'] = names(Diagnostics)[best.col]
      if(Recommendations['ShapiroWilks','ActionCode'] == 'EX') {
        Recommendations['ShapiroWilks','Comment'] = 'Exclude outliers to improve normality of residuals'
      } else {
        Recommendations['ShapiroWilks','Comment'] = 'Analyze transformed assessments to improve normality of residuals'
      }
    }
  }
  if(Diagnostics['Skewness','IID'] < 0.05) { 
    best.col <- which.min(abs(Diagnostics[c('SkewnessT'),-exclude.cols]))+1 
    if(Diagnostics[c('Skewness'),best.col] < 0.05 | !(best.col>2)) {
      Recommendations['Skewness','ActionCode'] = 'AR'
      Recommendations['Skewness','Comment'] = 'Exclude outliers to reduce skewness of residuals'
    } else {
      Recommendations['Skewness','ActionCode'] = names(Diagnostics)[best.col]
      if(Recommendations['Skewness','ActionCode'] == 'EX') {
        Recommendations['Skewness','Comment'] = 'Exclude outliers to reduce skewness of residuals'
      } else {
        Recommendations['Skewness','Comment'] = 'Select action code to reduce skewness of residuals'
      }
    }
  }
  if(Diagnostics['Kurtosis','IID'] < 0.05) { 
    best.col <- which.min(abs(Diagnostics[c('KurtosisT'),-exclude.cols]))+1 
    if(Diagnostics['Kurtosis',best.col] < 0.05 | !(best.col>2)) {
      Recommendations['Kurtosis','ActionCode'] = 'AR'
      Recommendations['Kurtosis','Comment'] = 'Available transformations do not correct kurtosis of residuals'
    } else {
      Recommendations['Kurtosis','ActionCode'] = names(Diagnostics)[best.col]
      if(Recommendations['Kurtosis','ActionCode'] == 'EX') {
        Recommendations['Kurtosis','Comment'] = 'Exclude outliers to reduce kurtosis of residuals'
      } else {
        Recommendations['Kurtosis','Comment'] = 'Select action code to reduce kurtosis of residuals'
      }
    }
  }
}
format(Recommendations, digits=4,scientific=FALSE)
```
# Standard Analysis
The standard analysis will be an analysis of variance of the linear model determined by the design of the experiment.
```{r}
#skipping standard analysis in diagnostic script
```
# Alternate Analysis
Alternative analysis suggested by the Recommendations may be included here.
```{r}
#skipping alternate analysis in diagnostic script
```
# Plot
```{r,echo=FALSE,fig.width=8,fig.height=8}
plot.columns <- 2
values1 <- arm.dat$assessment1
xlabA1 <- 'Assessment Value'
ylabB1 <- 'Assessment Values'
ylabD1 <- 'Assessment Value'
xlabD1 <- 'Treatment Mean'
mainA1 <- 'Histogram of Assessments'
treatments1 <- arm.dat$treatment
means <- tapply(values1,list(treatments1),mean,na.rm=TRUE)
means1 <- means[treatments1]
main1 <- 'Assessment 1 Box-Whisker'
code1 <- 'Raw'
not.IID <- which(Recommendations$ActionCode != 'IID')
recommended <- 1
if(length(not.IID)>0) {recommended <- which.min(not.IID)}
code2 <- Recommendations$ActionCode[recommended]
if(code2 %in% c('NM', 'ES', 'N')) plot.columns <- 1
if(plot.columns>1) {
  values2 <- Residuals$IID
  treatments2 <- Residuals$Treatment
  means2 <- predict(analysis.model)[Residuals$Plot]
  main2 <- 'AOV Residual Box-Whisker'
  xlabA2 <- 'Residual Value'
  ylabB2 <- 'Residual Values'
  ylabD2 <- 'Standardized Residuals'
  xlabD2 <- 'AOV Predicted Value'
  mainA2 <- 'Histogram of AOV Residuals'
  if(code2 == 'AR') {
    code2 <- 'IID'
  }
  if(!(code2 %in% c('IID','AR'))) {
    code1 <- 'IID'
    values1 <- values2
    means1 <- means2
    main1 <- main2
    mainA1 <- mainA2
    xlabA1 <- xlabA2
    ylabB1 <- ylabB2
    ylabD1 <- ylabD2
    xlabD1 <- xlabD2
      main2 <- paste(code2,'Residual Box-Whisker')
      mainA2 <- paste('Histogram of', code2, 'Residuals')
      ylabB2 <- paste(code2,'Residuals')
      ylabD2 <- paste(code2,'Std Residuals')
      xlabD2 <- paste(code2,'Predicted Value')
      treatments2 <- models[[code2]]$model$treatment
      means2 <- predict(models[[code2]])
      values2 <- residuals(models[[code2]])
  }
}
sigma.rng <- c(-3,-2,0,2,3)
line.col <- c('red','blue','darkgray','blue','red')
lbl1.skew <- paste('Skewness p = ', format(Diagnostics[c('Skewness'),code1],digits=3))
lbl1.kur <- paste('Kurtosis p = ', format(Diagnostics[c('Kurtosis'),code1],digits=3))
lbl1.sw <- paste('Shapiro-Wilks p = ', format(Diagnostics[c('ShapiroWilks'),code1],digits=3))
lbl1.lev <- paste('Levene p = ', format(Diagnostics[c('Levenes'),code1],digits=3))
if(plot.columns>1) {
  lbl2.skew <- paste('Skewness p = ', format(Diagnostics[c('Skewness'),code2],digits=3))
  lbl2.kur <- paste('Kurtosis p = ', format(Diagnostics[c('Kurtosis'),code2],digits=3))
  lbl2.sw <- paste('Shapiro-Wilks p = ', format(Diagnostics[c('ShapiroWilks'),code2],digits=3))
  lbl2.lev <- paste('Levene p = ', format(Diagnostics[c('Levenes'),code2],digits=3))
}
error.lwd <- 1
error.lty <- 2
stat.cex <- 0.9
par(mfrow=c(4,plot.columns))
plot(values1 ~ treatments1, main = main1,xlab='Treatment Number',ylab=ylabB1)
mtext(lbl1.lev,cex=stat.cex)
if(plot.columns>1) {
  plot(values2 ~ treatments2, main = main2, xlab='Treatment Number',ylab=ylabB2)
  mtext(lbl2.lev,cex=stat.cex)
}
qqnorm.plot <- function(values) {
  error.lwd <- 1
  error.lty <- 2
  good<-!is.na(values)
  qqnorm(values[good])
  #from https://stats.stackexchange.com/questions/111288/confidence-bands-for-qq-line
  ord<-order(values[good])
  ord.x<-values[good][ord]
  n <- length(values[good])
  P<-ppoints(n)
  z<-qnorm(P)
  coef<-coef(lm(ord.x~z))
  a<-coef[1]
  b<-coef[2]
  conf<- 1-alpha
  zz<-qnorm(1-(1-conf)/2)
  SE<-(b/dnorm(z))*sqrt(P*(1-P)/n)
  fit.value<-a+b*z
  abline(a,b,lty=1,lwd=2,col='gray')
  lines(z,fit.value+zz*SE,lty=error.lty,lwd=error.lwd,col='blue')
  lines(z,fit.value-zz*SE,lty=error.lty,lwd=error.lwd,col='blue')
  conf<- 1-alpha/2
  zz<-qnorm(1-(1-conf)/2)
  lines(z,fit.value+zz*SE,lty=error.lty,lwd=error.lwd,col='red')
  lines(z,fit.value-zz*SE,lty=error.lty,lwd=error.lwd,col='red')
}
qqnorm.plot(values1)
mtext(lbl1.sw, cex=stat.cex)
legend('topleft',inset=0.01, legend = c('2 SD', '3 SD'),col = c('blue', 'red'), lty = error.lty, box.lty=0)
if(plot.columns>1) {
  qqnorm.plot(values2)
  legend('topleft',inset=0.01, legend = c('2 SD', '3 SD'),col = c('blue', 'red'), lty = error.lty, box.lty=0)
  mtext(lbl2.sw, cex=stat.cex)
}
plot(x= means1,y=values1,xlab=xlabD1,ylab=ylabD1, main = 'Spread-Location')
abline(h=Diagnostics[c('Mean'),code1] + sigma.rng*Diagnostics[c('SD'),code1], col=line.col,lty=error.lty,lwd=error.lwd)
mtext(lbl1.skew, cex=stat.cex)
if(plot.columns>1) {
  plot(x=means2, y=values2,xlab=xlabD2,ylab=ylabD2, main = 'Spread-Location')
  abline(h=Diagnostics[c('Mean'),code2] + sigma.rng*Diagnostics[c('SD'),code2], col=line.col,lty=error.lty,lwd=error.lwd)
  mtext(lbl2.skew,cex=stat.cex)
}
hist(values1, xlab=xlabA1, col='lightgrey',main=mainA1)
abline(v=Diagnostics[c('Mean'),code1] + sigma.rng*Diagnostics[c('SD'),code1], col=line.col, lty=error.lty,lwd=error.lwd)
mtext(lbl1.kur,cex=stat.cex)
if(plot.columns>1) {
  hist(values2 , xlab=xlabA2, col='lightgrey',main=mainA2)
  abline(v=Diagnostics[c('Mean'),code2] + sigma.rng*Diagnostics[c('SD'),code2], col=line.col, lty=error.lty,lwd=error.lwd)
  mtext(lbl2.kur,cex=stat.cex)
}

```

# Analysis of Ranks

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

```{r}
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)
print(RankDF <- k-1)
print(RankP <- 1-pchisq(T,df=k-1))
print(CritChi <- qchisq(1 - alpha, df=k-1))

print(means <- R/n)

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)
print(post.tbl)
```

# Analysis in R
These code chunks are not included in ARM scripts.

In base R, we have

```{r}
kruskal.test(assessment1 ~ treatment, data=arm.dat)
```

`agricolae` also provides Kruskal-Wallis and mean comparisons
```{r}
library(agricolae)
res <- kruskal(arm.dat$assessment1, arm.dat$treatment, p.adj="none")
print(res)
```
