Example

Bread dough experiment

Test bread rise height caused by factor time:

We test 3 values for time with 4 replicates each (12 experimental units of 3 treatments with 4 replicates)

set.seed(234789)
f <- factor(rep(c(35,40,45),each = 4))
randomized_times <- sample(f,12)
allocation <- 1:12
design <- data.frame(allocated_loaf = allocation, time = randomized_times)
design

##    allocated_loaf time
## 1               1   40
## 2               2   45
## 3               3   35
## 4               4   40
## 5               5   40
## 6               6   45
## 7               7   45
## 8               8   45
## 9               9   35
## 10             10   40
## 11             11   35
## 12             12   35

Linear model used for CRD inference

Obviously we are only making inference about the mean effect within each treatment given the replications:
Cell mean model
\(Y_{ij} = \mu_i + \epsilon_{ij}\)

A more useful representation:
Effects model
\(Y_{ij} = \mu + \tau_i + \epsilon_{ij}\)
This formulation is more useful because \(\tau_i\) is interpreted as the difference of the mean within treatment \(i\) and the overall average.

If the design is executed properly the log likelihood function will be maximised if we minimise the sum of squares of residuals (linear regression).

BUT in R the lm() function will omit the first factor \(\tau_1\) and derive coefficients relative to \(\mu + \tau_1\):

bread_example <- daewr::bread

bread_example %>% lm(formula = height~time) %>% summary

## 
## Call:
## lm(formula = height ~ time, data = .)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.812 -1.141  0.000  1.266  2.250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4375     0.7655   7.104 5.65e-05 ***
## time40        2.8125     1.0825   2.598   0.0288 *  
## time45        2.8750     1.0825   2.656   0.0262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.531 on 9 degrees of freedom
## Multiple R-squared:  0.5056, Adjusted R-squared:  0.3958 
## F-statistic: 4.602 on 2 and 9 DF,  p-value: 0.042

So we can see the lm() function did not derive a coefficient for the time = 35 factor level. In this call the intercept is interpreted as \(\mu + \tau_1\) and the other factor levels are all evaluated relative to this intercept:

Obviously from this equation it’s easy to derive e.g. \(\mu + \tau_2\) by adding the 1st and second coefficients together…

Hypothesis test for No Treatment Effect

Now let’s test if the treatment (changing the time factor) has any affect on the response variable (the height of the baked bread):

\[\begin{align} H_0:& \tau_i = \tau_j \ \forall \ {i,j} \\ H_1:& \textrm{atleast 2 } \tau_i \textrm{'s differ} \end{align}\]

bread_example %>% aov(formula = height~time) %>% summary

##             Df Sum Sq Mean Sq F value Pr(>F)  
## time         2  21.57  10.786   4.602  0.042 *
## Residuals    9  21.09   2.344                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So we can see that the p-value here is the same as the p-value in the lm() summary and since it is < 0.05 we can safely reject \(H_0\) and say different time factor levels had a significant impact on the response variable height.

Verify linearity assumption

Of course, using a lm() model assumes linearity and we should validate this assumed relationship between time and height:

Assumptions validation function

Validate_assumptions <- function(df, formula_ = NULL, model, factor, response) {
  qq_plot <- 
    df %>% 
    ggplot(aes(sample = model %>% residuals))+
    stat_qq()+
    ggtitle("QQ plot of residuals")+
    xlab("")+
    ylab("")
  
  Res_fac <- df %>% 
    ggplot()+
    geom_point(aes(y = model %>% residuals(), x = df[,factor] ))+
    ggtitle("Residuals VS factor levels") +
    xlab("Residuals")+
    ylab("Factor levels")

  Res_fit <- df %>% 
    ggplot(aes(y = model %>% residuals(), x = model %>% predict))+
    geom_point()+
    ggtitle("Residuals VS fitted")+
    xlab("Residuals")+
    ylab("fitted")
  
  Res_eu <- df %>% 
    ggplot(aes(y = model %>% residuals(), x = seq(nrow(df))))+
    geom_point()+
    ggtitle("Residuals VS experimental unit")+
    xlab("Residuals")+
    ylab("Experimental unit")
  
  multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
  
  multiplot(qq_plot, Res_fac, Res_fit, Res_eu, cols=2)
  
}

Validate_assumptions(df = bread_example,model =  bread_example %>% lm(formula = height~time), factor = "time") 

# qqnorm(y = bread_example %>% lm(formula = height~time) %>% residuals)

In this case we do seem to have a fairly normal response relationship but with some variability in the lower corner. Consider that we only have 12 points.

A more serious violation of assumption is heterogeneity of variances. In general when the linearity assumption is violated so will the homogeniety of variances assumption. To remedy this we can apply a transformation to the data to reconsile the linearity assumption.

Box-Cox Power Transformations

** When variation is relative to the mean observation within each treatment **

This transformation will change the data so that it is distributed normally. Normally distributed data has honogeneous variances

These transformations ensure that \(\sigma \propto \mu^{1 - \lambda}\)

To find \(\lambda\) we can do:

box_cox <- boxcox(aov(bread_example,formula = height~time))

lambda = box_cox$x[which(box_cox$y == box_cox$y %>% max)]
lambda

## [1] -0.5050505

Now we transform the data:

bread_example_transformed <- 
  bread_example %>% 
  mutate(height = height^lambda)

Let’s redo the normality validation on residuals:

qqnorm(y = bread_example_transformed %>% lm(formula = height~time) %>% residuals())

And if we look at the hypothesis test with the now normally distributed data:

bread_example_transformed %>% aov(formula = height~time) %>% summary

##             Df  Sum Sq  Mean Sq F value Pr(>F)  
## time         2 0.01732 0.008662   6.134 0.0209 *
## Residuals    9 0.01271 0.001412                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Our new p-value is smaller which indicates that an observation as or more extreme than what we have is even more unlikely under the null hypothesis that there is no treatment effect.

Distributional transformations

** When variation is relative to the mean observation within each treatment **

An alternative to box-cox power transformation

OK, well let’s use the lognormal transformation and see the effect on the anova:

bread_example_lognormal <- 
  bread_example %>% 
  mutate(height = log(height))

Let’s redo the normality validation on residuals:

qqnorm(y = bread_example_lognormal %>% lm(formula = height~time) %>% residuals())

In this case we seem to fit the tail distribution better…

And if we look at the hypothesis test with the now more? normally distributed errors:

bread_example_lognormal %>% aov(formula = height~time) %>% summary

##             Df Sum Sq Mean Sq F value Pr(>F)  
## time         2 0.4598 0.22992    5.64 0.0258 *
## Residuals    9 0.3669 0.04077                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The p-value is again shifted towards 0.0258 which indicates with higher certainty that we reject the null hypothesis

Weighted linear model

When we believe the heterogeneous variances are not related to the treatment mean we can use a weighted least squares

Use sd() as weights within treatments:

bread_example_Weighted <- 
    bread_example %>% 
  group_by(time) %>% 
  mutate(sd_treatment = 1/sd(height))
bread_example_Weighted

## # A tibble: 12 x 4
## # Groups:   time [3]
##     loaf time  height sd_treatment
##    <dbl> <fct>  <dbl>        <dbl>
##  1  1.00 35      4.50        1.04 
##  2  2.00 35      5.00        1.04 
##  3  3.00 35      5.50        1.04 
##  4  4.00 35      6.75        1.04 
##  5  5.00 40      6.50        0.485
##  6  6.00 40      6.50        0.485
##  7  7.00 40     10.5         0.485
##  8  8.00 40      9.50        0.485
##  9  9.00 45      9.75        0.735
## 10 10.0  45      8.75        0.735
## 11 11.0  45      6.50        0.735
## 12 12.0  45      8.25        0.735

Now that we have the weights for each treatment we can use it in a weighted OLS:

bread_example_Weighted %>% lm(formula = height~time,weights = bread_example_Weighted$sd_treatment) %>% summary

## 
## Call:
## lm(formula = height ~ time, data = ., weights = bread_example_Weighted$sd_treatment)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5543 -1.0203  0.0050  0.9611  1.5671 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4375     0.5941   9.152 7.44e-06 ***
## time40        2.8125     1.0520   2.674   0.0255 *  
## time45        2.8750     0.9220   3.118   0.0124 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.209 on 9 degrees of freedom
## Multiple R-squared:  0.5805, Adjusted R-squared:  0.4872 
## F-statistic: 6.226 on 2 and 9 DF,  p-value: 0.02006

This gives us a p-value of 0.02006 which is very similar to the other transformed methods

Creating a weighted OLS function

So to generalize this for simple case of CRD:

weighted_ols_anova <-  function(df, formula_, factor, response) {
  factor_enq <- dplyr::enquo(factor)
  response_enq <- dplyr::enquo(response)
  # factor_quo <- dplyr::quo(factor)
  # response_quo <- dplyr::quo(response)
  df <- 
    df %>% 
  group_by(!!factor_enq) %>%
  mutate(sd_treatment = 1/sd(!!response_enq)) %>% 
    ungroup
  
   # y = df %>% select(!!response_enq) %>% flatten_dbl
   # x = df %>% select(!!factor_enq) %>% flatten_dbl
  
    df %>% lm(formula =  formula_, weights = sd_treatment) %>% anova
    # df %>% lm(y = y, x = x, weights = sd_treatment) %>% anova

  # return(df)
}

weighted_ols_anova(bread,height~time,time,height)

## Analysis of Variance Table
## 
## Response: height
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## time       2 18.209  9.1047  6.2263 0.02006 *
## Residuals  9 13.161  1.4623                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Generalized linear models

Pol data using logistic link function assuming multinomial distribution

Participants were asked to rank a treatment from 1 to 5

Data:

mod_full <- MASS::polr(data = daewr::teach,
                   formula = score ~ method,
                   weight = count)

mod_reduced <- MASS::polr(data = daewr::teach,
                   formula = score ~ 1,
                   weight = count)

anova(mod_full,mod_reduced)

## Likelihood ratio tests of ordinal regression models
## 
## Response: score
##    Model Resid. df Resid. Dev   Test    Df LR stat.     Pr(Chi)
## 1      1       294   885.9465                                  
## 2 method       292   876.2986 1 vs 2     2  9.64787 0.008035108

The reduced model predicts assuming only a global mean which conceptually means:
\(y_{i,j} = \mu + \epsilon_{i,j} \ \ \ VS \ \ \ y_{i,j}=\mu + \tau_i + \epsilon_{i,j}\)

This shows that using different methods does make a difference compared to no particular method used. Because the p-value here rejects the null hypothesis that the two models are the same and therefore rejects the null hypothesis that all treatment effects are the same: \(\mu_i = \mu_j \ \forall \ {i,j}\)

How do we measure treatment effect?

Let’s visualize the scores by method used

# Predicted <- 
#   daewr::teach %>% 
#   mutate(predicted_score = mod_full %>% predict)

# Predicted %>% 
#   group_by(method) %>% 
#   nest %>% 
#   mutate(Histogram = data %>% map(~hist(.$score))) %>% 
#   .Histogram

daewr::teach %>% 
  ggplot()+
  geom_bar(aes(x = score,y = count),stat = "identity")+
  facet_wrap(~method)+
  ggtitle("Distribution of scores by method",subtitle = "Method 3 appears to outperform")

daewr::teach %>% 
  mutate_if(is.factor,as.numeric) %>% 
  group_by(method) %>% 
  summarise(mean_score = sum(score*count)/sum(count),
            Total_votes = sum(count))

## # A tibble: 3 x 3
##   method mean_score Total_votes
##    <dbl>      <dbl>       <dbl>
## 1   1.00       2.98        88.0
## 2   2.00       3.22       104  
## 3   3.00       3.47       106

We can see that on average the 3rd treatment here would outperform the rest of the treatment groups from this designed experiment.

Test for assumption violations

We can’t repeat the usual battery of tests here since we are dealing with a categorical response and the idea of residual error should be defined differently from say \(y-\hat{y}\).

In order to run the usual battery of validity tests we need the residuals of the fitted model. However, this is a categorical response experiment.

# Validate_assumptions(df = daewr::teach, model = mod_full)

Deternining the number of replicates

Bread dough experiment

The idea of practical difference is important here. If the practitioner wants to monitor any difference in cell means larger than say 3 inches he would use the Fpower1 function with \(\Delta = 3\).

Fpower1 calculates the power of the 1 way anova F distribution

Now we specify the number of replications we would like to calculate the power for given we know the variance of the experimental error \(\hat{\sigma}^2\)

If we assume the naive OLS this value is:

daewr::bread %>% lm(formula = height~time) %>% residuals() %>% var

## [1] 1.917614

Let’s assume the value \(\hat{\sigma}^2 = 1.917614\)

power <- Fpower1(alpha = rep(0.05, 6-2+1),
                 nlev = 3,
                 nreps = 2:6, 
                 Delta = 3,
                 sigma = sqrt(1.917614)
                )
power

##      alpha nlev nreps Delta    sigma     power
## [1,]  0.05    3     2     3 1.384779 0.2086501
## [2,]  0.05    3     3     3 1.384779 0.4366902
## [3,]  0.05    3     4     3 1.384779 0.6315614
## [4,]  0.05    3     5     3 1.384779 0.7730600
## [5,]  0.05    3     6     3 1.384779 0.8666144

This table then shows us that with as many as 6 replicates per treatment we can detect with ~87% accuracy a cell mean shift of 3 inches.

Based on the affordability in terms of time and money we can limit the factor under investigation using screening designs and perhaps limit the number of factors such that we get enough power out of the affordable number of replications.

Choosing optimal treatments

In general it is possible to construct hypothesis tests of comparisons before the design of experiments is carried out and responses are recorded but we will not focus on that here.

Instead we will assume the practitioner is interested in comparing treatments to one another in an unbiased way or compared to the current control treatment being used.

Unplanned Comparisons

When comparing specific treatments with one another after the experiment has been designed and carried out we must adjust for bias.

In R we can use the TukeyHSD function to compare cell mean shifts between different treatments

Sugarbeet dataset

This dataset is another CRD design. Let’s illustrate the TukeyHSD function here:

Sugerbeet_anova <- aov(yield ~ treat, data = daewr::sugarbeet)

Sugerbeet_anova

## Call:
##    aov(formula = yield ~ treat, data = daewr::sugarbeet)
## 
## Terms:
##                    treat Residuals
## Sum of Squares  291.0050   29.5865
## Deg. of Freedom        3        14
## 
## Residual standard error: 1.453727
## Estimated effects may be unbalanced

Sugerbeet_TukeyHSD <- TukeyHSD(Sugerbeet_anova, ordered = TRUE)

Sugerbeet_TukeyHSD

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = yield ~ treat, data = daewr::sugarbeet)
## 
## $treat
##     diff        lwr       upr     p adj
## B-A  6.3  3.3122236  9.287776 0.0001366
## D-A 10.0  7.1655464 12.834454 0.0000004
## C-A 10.1  7.2655464 12.934454 0.0000003
## D-B  3.7  0.8655464  6.534454 0.0094231
## C-B  3.8  0.9655464  6.634454 0.0077551
## C-D  0.1 -2.5723484  2.772348 0.9995162

Here we can see that treatments C and D do not reject the null hypothesis that they have different cell means. However, the test also shows that no other treatments are considered the same and all reject at 5%.

Based purely on the confidence interval showed we could conclude that treatments C and D appear to be the most optimal if a higher mean in the response was required.

This function does not explicitly tell us what those means are, merely their differences.

Another, more conservative, function we can use is the SNK.test function from the agricolae package:

Sugarbeet_SNK <- SNK.test(Sugerbeet_anova, "treat", alpha = 0.05)

Sugerbeet_anova

## Call:
##    aov(formula = yield ~ treat, data = daewr::sugarbeet)
## 
## Terms:
##                    treat Residuals
## Sum of Squares  291.0050   29.5865
## Deg. of Freedom        3        14
## 
## Residual standard error: 1.453727
## Estimated effects may be unbalanced

This time we can clearly see the treatment means along with their summary statistics. As for rejection of the null hypothesis we can compare the M column in the groups output. Treatments that did not reject the null would have the same letters here.

Comparison to the current control treatment

To compare different treatments to the current control treatment we can use the function glht from the multcomp package

The glht function assumes that the 1st factor level is in fact the control treatment

Sugarbeet_glht <- glht(Sugerbeet_anova,linfct = mcp(treat = "Dunnett"), alternative = "greater")

Sugarbeet_glht %>% summary

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = yield ~ treat, data = daewr::sugarbeet)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>t)    
## B - A <= 0   6.3000     1.0279   6.129 <1e-04 ***
## C - A <= 0  10.1000     0.9752  10.357 <1e-04 ***
## D - A <= 0  10.0000     0.9752  10.254 <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

By specifying alternative = "greater" we can test the null hypothesis of the cell mean being greater than that of the control treatment using the Dunnett method.

Here we can again see that treatments C and D are superior to the control treatment A.