TPSS 603 - Experimental Design

Ashlyn Godbehere

1. Create a potential randomization for a CR experiment that has 3 treatments and 5 replications per treatment.

#create a randomization for a CR design 
#design.crd() #function from agricolae that will do the randomization

trt <- c("trt1", "trt2", "trt3")  # Only 3 treatments
r <- c(5, 5, 5)                   # 5 replications for each treatment

# Create randomization for CR design
design.crd(trt, r)

## $parameters
## $parameters$design
## [1] "crd"
## 
## $parameters$trt
## [1] "trt1" "trt2" "trt3"
## 
## $parameters$r
## [1] 5 5 5
## 
## $parameters$serie
## [1] 2
## 
## $parameters$seed
## [1] -6067656
## 
## $parameters$kinds
## [1] "Super-Duper"
## 
## $parameters[[7]]
## [1] TRUE
## 
## 
## $book
##    plots r  trt
## 1    101 1 trt3
## 2    102 2 trt3
## 3    103 1 trt1
## 4    104 1 trt2
## 5    105 2 trt1
## 6    106 2 trt2
## 7    107 3 trt1
## 8    108 3 trt3
## 9    109 4 trt3
## 10   110 5 trt3
## 11   111 3 trt2
## 12   112 4 trt1
## 13   113 4 trt2
## 14   114 5 trt1
## 15   115 5 trt2

2. What is the Null Hypothesis? Answer: A hypothesis where there is no difference between the variables being tested. All fertilizer treatments have the same effect on plant height. H₀: μ₁ = μ₂ = μ₃ = … = μₜ

3. What is the Alternative hypothesis? Answer: A hypothesis where there is a difference or effect between the variables being tested. At least one fertilizer treatment has a different effect on plant height. Hₐ: At least one μᵢ ≠ μⱼ (where i ≠ j) or Hₐ: Not all μᵢ are equal

4. What is the Linear additive model? # this is the abstraction of the generalized model Answer: A linear additive model is where the variable being measured is explained by adding together different effects. This model says each plant height measurement equals the overall average plant height, plus the fertilizer treatment effect, plus random experimental error. Yᵢₖ = μ + αₖ + εᵢₖ μ = grand mean aₖ = an effect of treatment for group k εᵢₖ = a person i’s residual within group k

#this is where you use the data from the file attached
setwd("/Users/ashlyngodbehere/Desktop/R/R_class")
read.csv("~/Desktop/R/R_class/CR_design.csv")

##      Trt Plant_height
## 1  Fert1           46
## 2  Fert1           42
## 3  Fert1           44
## 4  Fert1           43
## 5  Fert1           45
## 6  Fert1           49
## 7  Fert2           15
## 8  Fert2           16
## 9  Fert2            9
## 10 Fert2           10
## 11 Fert2           11
## 12 Fert2           14
## 13 Fert3           26
## 14 Fert3           24
## 15 Fert3           29
## 16 Fert3           28
## 17 Fert3           27
## 18 Fert3           31

cr=read.csv("CR_design.csv")                    
head(cr,6) #look at the first 6 lines of the data

##     Trt Plant_height
## 1 Fert1           46
## 2 Fert1           42
## 3 Fert1           44
## 4 Fert1           43
## 5 Fert1           45
## 6 Fert1           49

#create the anova model for the design using the data provided
mod <- aov(Plant_height ~ Trt, data=cr) #create a linear model for the data

5. What are the treatments? Answer: The treatments are the fertilizers: Fert1, Fert2, and Fert3.

6. How many factors are being studied? Answer: Only one, fertilizer type.

7. ANOVA table for the experiment.

summary(mod) #Print the summary of the anova table from the model abve

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Trt          2 3141.8  1570.9   231.4 5.33e-12 ***
## Residuals   15  101.8     6.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8. Is the normalilty of residuals assumption true (qqplot)? Answer:
True. The Q-Q plot shows points following the diagonal line with some minor spots not on the line.

par(mfrow=c(2,2)) #set the plot
plot(mod) #use the plot funciton to make the qq plot

9. Do we reject or fail to reject the null hypothesis? Answer:
We reject the null hypothesis because p-value (5.33e-12) < 0.05. This means the fertilizers have significantly different effects on plant height.

10. Construct a contrast between two treatments and interpret the result. Answer: The contrast between Fert1 and Fert2 shows that Fert1 produces plants that are 32.33 cm taller on average(p < 0.05), this means Fert1 is significantly better than Fert2 for plant growth.

mod # look at linear model

## Call:
##    aov(formula = Plant_height ~ Trt, data = cr)
## 
## Terms:
##                       Trt Residuals
## Sum of Squares  3141.7778  101.8333
## Deg. of Freedom         2        15
## 
## Residual standard error: 2.60555
## Estimated effects may be unbalanced

summary(mod) # look at the summary of the model

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Trt          2 3141.8  1570.9   231.4 5.33e-12 ***
## Residuals   15  101.8     6.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit.contrast(mod,"Trt", c(1,-1,0)) # create the contrast between two treatments

##                  Estimate Std. Error  t value     Pr(>|t|)
## Trt c=( 1 -1 0 ) 32.33333   1.504315 21.49373 1.108421e-12

fit.contrast(mod,"Trt", c(1,0,-1))

##                  Estimate Std. Error  t value     Pr(>|t|)
## Trt c=( 1 0 -1 ) 17.33333   1.504315 11.52241 7.515129e-09

fit.contrast(mod,"Trt", c(0,1,-1))

##                  Estimate Std. Error   t value     Pr(>|t|)
## Trt c=( 0 1 -1 )      -15   1.504315 -9.971317 5.189285e-08

Answer:

Construct all pairwise comparisons using an HSD or LSD, create a plot, and interpret the results

# Tukey's HSD (Honest Significant Difference) - more conservative
tukey_results <- TukeyHSD(mod)
print(tukey_results)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Plant_height ~ Trt, data = cr)
## 
## $Trt
##                  diff       lwr       upr p adj
## Fert2-Fert1 -32.33333 -36.24075 -28.42592 0e+00
## Fert3-Fert1 -17.33333 -21.24075 -13.42592 0e+00
## Fert3-Fert2  15.00000  11.09259  18.90741 1e-07

# Fisher's LSD (Least Significant Difference) - less conservative
library(agricolae)
lsd_results <- LSD.test(mod, "Trt", p.adj = "none")
print(lsd_results)

## $statistics
##    MSerror Df     Mean       CV t.value      LSD
##   6.788889 15 28.27778 9.214124 2.13145 3.206371
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none    Trt   3  0.05
## 
## $means
##       Plant_height      std r       se      LCL      UCL Min Max   Q25  Q50
## Fert1     44.83333 2.483277 6 1.063711 42.56609 47.10058  42  49 43.25 44.5
## Fert2     12.50000 2.880972 6 1.063711 10.23275 14.76725   9  16 10.25 12.5
## Fert3     27.50000 2.428992 6 1.063711 25.23275 29.76725  24  31 26.25 27.5
##         Q75
## Fert1 45.75
## Fert2 14.75
## Fert3 28.75
## 
## $comparison
## NULL
## 
## $groups
##       Plant_height groups
## Fert1     44.83333      a
## Fert3     27.50000      b
## Fert2     12.50000      c
## 
## attr(,"class")
## [1] "group"

plot(tukey_results)

plot(lsd_results)

Answer: All three fertilizers are significantly different. Fert1 gave the tallest plants (44.83 cm), Fert3 was middle (27.50 cm), and Fert2 gave the shortest plants (12.50 cm). All comparisons had p < 0.001, so each fertilizer works differently. Fert1 is clearly the best choice.

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.