DEVHIDE
Login Or Sign up

In R, how to analyze Strip-Plot Design when there are three factors?

114 Views Asked by At

I have one dataset like below

dataA= structure(list(sowing_date= c("Early", "Early", "Early", 
                                      "Early", "Early", "Early", 
                                      "Early", "Early", "Early", 
                                      "Early", "Early", "Early", 
                                      "Late", "Late", "Late", 
                                      "Late", "Late", "Late", 
                                      "Late", "Late", "Late", 
                                      "Late", "Late", "Late"), 
                      herbicide= c("H0", "H0", "H0", "H0", "H0", 
                                    "H0", "H1", "H1", "H1", "H1", 
                                    "H1", "H1", "H0", "H0", "H0", 
                                    "H0", "H0", "H0", "H1", "H1",
                                    "H1", "H1", "H1", "H1"), 
                      nitrogen= c("N0", "N0", "N0", "N1", "N1", 
                                   "N1", "N0", "N0", "N0", "N1", 
                                   "N1", "N1", "N0", "N0", "N0", 
                                   "N1", "N1", "N1", "N0", "N0", 
                                   "N0", "N1", "N1", "N1"), 
                      Block= c("Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3",
                               "Block 1", "Block 2", "Block 3"), 
                      Yield= c(30, 27, 25, 40, 41, 42, 37, 38, 
                               40, 48, 47, 46, 25, 27, 26, 41, 
                               41, 42,38, 39, 42, 57, 59, 60)), 
                 class= "data.frame", 
                 row.names= c(NA, -24L))

   sowing_date herbicide nitrogen   Block Yield
1        Early        H0       N0 Block 1    30
2        Early        H0       N0 Block 2    27
3        Early        H0       N0 Block 3    25
4        Early        H0       N1 Block 1    40
5        Early        H0       N1 Block 2    41
6        Early        H0       N1 Block 3    42
7        Early        H1       N0 Block 1    37
8        Early        H1       N0 Block 2    38
9        Early        H1       N0 Block 3    40
10       Early        H1       N1 Block 1    48
11       Early        H1       N1 Block 2    47
12       Early        H1       N1 Block 3    46
13        Late        H0       N0 Block 1    25
14        Late        H0       N0 Block 2    27
15        Late        H0       N0 Block 3    26
16        Late        H0       N1 Block 1    41
17        Late        H0       N1 Block 2    41
18        Late        H0       N1 Block 3    42
19        Late        H1       N0 Block 1    38
20        Late        H1       N0 Block 2    39
21        Late        H1       N0 Block 3    42
22        Late        H1       N1 Block 1    57
23        Late        H1       N1 Block 2    59
24        Late        H1       N1 Block 3    60

Now, my experimental design is strip-plot design, and the corresponding R code is below

library(agricolae)
model= with(dataA, strip.plot(BLOCK=Block,COL=herbicide, ROW=sowing_date, Yield))

However, this model is when there are only two factors like below layout

enter image description here

I have three factors; sowing_date, herbicide and nitrogen. I'd like to randomize nitrogen within herbicide and sowing_date like below layout.

enter image description here

Could you let me know how to analyze strip-plot design when there are three factors? Maybe, I can use split-split plot design, but the layout is designed as strip-plot, so I think it would be more accurate to use the model for strip-plot.

Simply, I added nitrogen in the model, but it doesn't work

    model= with(dataA, strip.plot(BLOCK=Block,COL=herbicide, ROW=sowing_date, nitrogen, Yield))

    Error in strip.plot(BLOCK = Block, COL = herbicide, ROW = sowing_date,  : 
    unused argument (Yield)

Could you answer my questions?

Thanks,

r plot statistics experimental-design
Original Q&A
1

There are 1 best solutions below

3
jared_mamrot On

Original answer

As @Dave2e commented, it sounds like you're looking for a split-split plot design. You can implement it using the ssp.plot() function from the agricolae package:

Example data:

dataA= structure(list(sowing_date= c("Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late"), 
                      herbicide= c("H0", "H0", "H0", "H0", "H0", 
                                   "H0", "H1", "H1", "H1", "H1", 
                                   "H1", "H1", "H0", "H0", "H0", 
                                   "H0", "H0", "H0", "H1", "H1",
                                   "H1", "H1", "H1", "H1"), 
                      nitrogen= c("N0", "N0", "N0", "N1", "N1", 
                                  "N1", "N0", "N0", "N0", "N1", 
                                  "N1", "N1", "N0", "N0", "N0", 
                                  "N1", "N1", "N1", "N0", "N0", 
                                  "N0", "N1", "N1", "N1"), 
                      Block= c("Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3",
                               "Block 1", "Block 2", "Block 3"), 
                      Yield= c(30, 27, 25, 40, 41, 42, 37, 38, 
                               40, 48, 47, 46, 25, 27, 26, 41, 
                               41, 42,38, 39, 42, 57, 59, 60)), 
                 class= "data.frame", 
                 row.names= c(NA, -24L))

dataA
#>    sowing_date herbicide nitrogen   Block Yield
#> 1        Early        H0       N0 Block 1    30
#> 2        Early        H0       N0 Block 2    27
#> 3        Early        H0       N0 Block 3    25
#> 4        Early        H0       N1 Block 1    40
#> 5        Early        H0       N1 Block 2    41
#> 6        Early        H0       N1 Block 3    42
#> 7        Early        H1       N0 Block 1    37
#> 8        Early        H1       N0 Block 2    38
#> 9        Early        H1       N0 Block 3    40
#> 10       Early        H1       N1 Block 1    48
#> 11       Early        H1       N1 Block 2    47
#> 12       Early        H1       N1 Block 3    46
#> 13        Late        H0       N0 Block 1    25
#> 14        Late        H0       N0 Block 2    27
#> 15        Late        H0       N0 Block 3    26
#> 16        Late        H0       N1 Block 1    41
#> 17        Late        H0       N1 Block 2    41
#> 18        Late        H0       N1 Block 3    42
#> 19        Late        H1       N0 Block 1    38
#> 20        Late        H1       N0 Block 2    39
#> 21        Late        H1       N0 Block 3    42
#> 22        Late        H1       N1 Block 1    57
#> 23        Late        H1       N1 Block 2    59
#> 24        Late        H1       N1 Block 3    60

Example code:

# install.packages("agricolae")
library(agricolae)
model= with(dataA, ssp.plot(block = Block, pplot = herbicide, splot = sowing_date, ssplot = nitrogen, Y = Yield))
#> 
#> ANALYSIS SPLIT-SPLIT PLOT:  Yield 
#> Class level information
#> 
#> herbicide    :  H0 H1 
#> sowing_date  :  Early Late 
#> nitrogen     :  N0 N1 
#> Block    :  Block 1 Block 2 Block 3 
#> 
#> Number of observations:  24 
#> 
#> Analysis of Variance Table
#> 
#> Response: Yield
#>                                Df  Sum Sq Mean Sq  F value    Pr(>F)    
#> Block                           2    3.08    1.54   0.5873  0.578110    
#> herbicide                       1  864.00  864.00 329.1429  0.003024 ** 
#> Ea                              2    5.25    2.63                       
#> sowing_date                     1   54.00   54.00  27.0000  0.006533 ** 
#> herbicide:sowing_date           1   73.50   73.50  36.7500  0.003738 ** 
#> Eb                              4    8.00    2.00                       
#> nitrogen                        1 1204.17 1204.17 458.7302 2.377e-08 ***
#> nitrogen:herbicide              1    0.67    0.67   0.2540  0.627877    
#> nitrogen:sowing_date            1   54.00   54.00  20.5714  0.001910 ** 
#> nitrogen:herbicide:sowing_date  1   28.17   28.17  10.7302  0.011260 *  
#> Ec                              8   21.00    2.63                       
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> cv(a) = 4.1 %, cv(b) = 3.5 %, cv(c) = 4.1 %, Mean = 39.91667

Created on 2023-09-04 with reprex v2.0.2

Edit 1

This is outside my area of expertise, so I'm not sure if this helps or not, but if you don't want to use a split-split plot design, you could combine herbicide and nitrogen 'status' into a single factor and get 'reasonable' looking results from a two-factor split-plot design:

dataA= structure(list(sowing_date= c("Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Early", "Early", "Early", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late", 
                                     "Late", "Late", "Late"), 
                      herbicide= c("H0", "H0", "H0", "H0", "H0", 
                                   "H0", "H1", "H1", "H1", "H1", 
                                   "H1", "H1", "H0", "H0", "H0", 
                                   "H0", "H0", "H0", "H1", "H1",
                                   "H1", "H1", "H1", "H1"), 
                      nitrogen= c("N0", "N0", "N0", "N1", "N1", 
                                  "N1", "N0", "N0", "N0", "N1", 
                                  "N1", "N1", "N0", "N0", "N0", 
                                  "N1", "N1", "N1", "N0", "N0", 
                                  "N0", "N1", "N1", "N1"), 
                      Block= c("Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3", 
                               "Block 1", "Block 2", "Block 3",
                               "Block 1", "Block 2", "Block 3"), 
                      Yield= c(30, 27, 25, 40, 41, 42, 37, 38, 
                               40, 48, 47, 46, 25, 27, 26, 41, 
                               41, 42,38, 39, 42, 57, 59, 60)), 
                 class= "data.frame", 
                 row.names= c(NA, -24L))

# install.packages("agricolae")
library(agricolae)

# original model
model= with(dataA, strip.plot(BLOCK=Block,COL=herbicide, ROW=sowing_date, Yield))
#> 
#> ANALYSIS STRIP PLOT:  Yield 
#> Class level information
#> 
#> herbicide    :  H0 H1 
#> sowing_date  :  Early Late 
#> Block    :  Block 1 Block 2 Block 3 
#> 
#> Number of observations:  24 
#> 
#> model Y: Yield ~ Block + herbicide + Ea + sowing_date + Eb + sowing_date:herbicide + Ec 
#> 
#> Analysis of Variance Table
#> 
#> Response: Yield
#>                       Df  Sum Sq Mean Sq  F value   Pr(>F)   
#> Block                  2    3.08    1.54                     
#> herbicide              1  864.00  864.00 329.1429 0.003024 **
#> Ea                     2    5.25    2.63                     
#> sowing_date            1   54.00   54.00  13.9355 0.064856 . 
#> Eb                     2    7.75    3.88                     
#> sowing_date:herbicide  1   73.50   73.50   0.7865 0.390127   
#> Ec                    14 1308.25   93.45                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> cv(a) = 4.1 %, cv(b) = 4.9 %, cv(c) = 24.2 %, Mean = 39.91667

# split-split plot
model= with(dataA, ssp.plot(block = Block, pplot = herbicide, splot = sowing_date, ssplot = nitrogen, Y = Yield))
#> 
#> ANALYSIS SPLIT-SPLIT PLOT:  Yield 
#> Class level information
#> 
#> herbicide    :  H0 H1 
#> sowing_date  :  Early Late 
#> nitrogen     :  N0 N1 
#> Block    :  Block 1 Block 2 Block 3 
#> 
#> Number of observations:  24 
#> 
#> Analysis of Variance Table
#> 
#> Response: Yield
#>                                Df  Sum Sq Mean Sq  F value    Pr(>F)    
#> Block                           2    3.08    1.54   0.5873  0.578110    
#> herbicide                       1  864.00  864.00 329.1429  0.003024 ** 
#> Ea                              2    5.25    2.63                       
#> sowing_date                     1   54.00   54.00  27.0000  0.006533 ** 
#> herbicide:sowing_date           1   73.50   73.50  36.7500  0.003738 ** 
#> Eb                              4    8.00    2.00                       
#> nitrogen                        1 1204.17 1204.17 458.7302 2.377e-08 ***
#> nitrogen:herbicide              1    0.67    0.67   0.2540  0.627877    
#> nitrogen:sowing_date            1   54.00   54.00  20.5714  0.001910 ** 
#> nitrogen:herbicide:sowing_date  1   28.17   28.17  10.7302  0.011260 *  
#> Ec                              8   21.00    2.63                       
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> cv(a) = 4.1 %, cv(b) = 3.5 %, cv(c) = 4.1 %, Mean = 39.91667

# 2-factor design with nitrogen instead of herbicide
model= with(dataA, strip.plot(BLOCK=Block,COL=nitrogen, ROW=sowing_date, Yield))
#> 
#> ANALYSIS STRIP PLOT:  Yield 
#> Class level information
#> 
#> nitrogen     :  N0 N1 
#> sowing_date  :  Early Late 
#> Block    :  Block 1 Block 2 Block 3 
#> 
#> Number of observations:  24 
#> 
#> model Y: Yield ~ Block + nitrogen + Ea + sowing_date + Eb + sowing_date:nitrogen + Ec 
#> 
#> Analysis of Variance Table
#> 
#> Response: Yield
#>                      Df  Sum Sq Mean Sq    F value   Pr(>F)    
#> Block                 2    3.08    1.54                        
#> nitrogen              1 1204.17 1204.17 28900.0000 3.46e-05 ***
#> Ea                    2    0.08    0.04                        
#> sowing_date           1   54.00   54.00    13.9355  0.06486 .  
#> Eb                    2    7.75    3.88                        
#> sowing_date:nitrogen  1   54.00   54.00     0.7615  0.39757    
#> Ec                   14  992.75   70.91                        
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> cv(a) = 0.5 %, cv(b) = 4.9 %, cv(c) = 21.1 %, Mean = 39.91667

# create a new variable (herb_nitro)
dataA$herb_nitro <- ifelse(dataA$herbicide == "H0" & dataA$nitrogen == "N0", "H0N0",
                           ifelse(dataA$herbicide == "H0" & dataA$nitrogen == "N1", "H0N1",
                                  ifelse(dataA$herbicide == "H1" & dataA$nitrogen == "N0", "H1N0",
                                         ifelse(dataA$herbicide == "H1" & dataA$nitrogen == "N1", "H1N1", "NA"))))
dataA
#>    sowing_date herbicide nitrogen   Block Yield herb_nitro
#> 1        Early        H0       N0 Block 1    30       H0N0
#> 2        Early        H0       N0 Block 2    27       H0N0
#> 3        Early        H0       N0 Block 3    25       H0N0
#> 4        Early        H0       N1 Block 1    40       H0N1
#> 5        Early        H0       N1 Block 2    41       H0N1
#> 6        Early        H0       N1 Block 3    42       H0N1
#> 7        Early        H1       N0 Block 1    37       H1N0
#> 8        Early        H1       N0 Block 2    38       H1N0
#> 9        Early        H1       N0 Block 3    40       H1N0
#> 10       Early        H1       N1 Block 1    48       H1N1
#> 11       Early        H1       N1 Block 2    47       H1N1
#> 12       Early        H1       N1 Block 3    46       H1N1
#> 13        Late        H0       N0 Block 1    25       H0N0
#> 14        Late        H0       N0 Block 2    27       H0N0
#> 15        Late        H0       N0 Block 3    26       H0N0
#> 16        Late        H0       N1 Block 1    41       H0N1
#> 17        Late        H0       N1 Block 2    41       H0N1
#> 18        Late        H0       N1 Block 3    42       H0N1
#> 19        Late        H1       N0 Block 1    38       H1N0
#> 20        Late        H1       N0 Block 2    39       H1N0
#> 21        Late        H1       N0 Block 3    42       H1N0
#> 22        Late        H1       N1 Block 1    57       H1N1
#> 23        Late        H1       N1 Block 2    59       H1N1
#> 24        Late        H1       N1 Block 3    60       H1N1

model2 <- with(dataA, strip.plot(BLOCK=Block,COL=herb_nitro, ROW=sowing_date, Yield))
#> 
#> ANALYSIS STRIP PLOT:  Yield 
#> Class level information
#> 
#> herb_nitro   :  H0N0 H0N1 H1N0 H1N1 
#> sowing_date  :  Early Late 
#> Block    :  Block 1 Block 2 Block 3 
#> 
#> Number of observations:  24 
#> 
#> model Y: Yield ~ Block + herb_nitro + Ea + sowing_date + Eb + sowing_date:herb_nitro + Ec 
#> 
#> Analysis of Variance Table
#> 
#> Response: Yield
#>                        Df  Sum Sq Mean Sq F value    Pr(>F)    
#> Block                   2    3.08    1.54                      
#> herb_nitro              3 2068.83  689.61 244.591 1.164e-06 ***
#> Ea                      6   16.92    2.82                      
#> sowing_date             1   54.00   54.00  13.935 0.0648562 .  
#> Eb                      2    7.75    3.88                      
#> sowing_date:herb_nitro  3  155.67   51.89  32.487 0.0004173 ***
#> Ec                      6    9.58    1.60                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> cv(a) = 4.2 %, cv(b) = 4.9 %, cv(c) = 3.2 %, Mean = 39.91667

Created on 2023-09-04 with reprex v2.0.2

Related Questions in R

Related Questions in PLOT

Related Questions in STATISTICS

Related Questions in EXPERIMENTAL-DESIGN

Trending Questions

Popular # Hahtags

javascript python java c# php android html jquery c++ css ios sql mysql r reactjs

Popular Questions

.

Copyright © 2021 Jogjafile Inc.