5  Analysis Across Sites

5.1 Example 5.1 (Pg. 68)

Example 5.1 (Pg. 68)

In 1985 species/provenance trials were laid out at six sites in Thailand as part of an ACIAR project extending over several years to investigate Australian multi-purpose tree species. The experimental design in each case was an RCB design with three replicates and the number of seedlots ranged from 30 to 42. Plots consisted of 5 \(\times\) 5 trees with a 2m \(\times\) 2m spacing. Plot summary files were constructed for the 24-month measurement according to the methods described in Chapter 3. Analyses were performed on the plot mean data at each site.

 library(car)

 library(dae)

 library(dplyr)

 library(emmeans)

 library(ggplot2)

 library(lmerTest)

 library(magrittr)

 library(predictmeans)

 data(DataExam5.1)

 # Pg.68
 fm5.4 <-
       lm(
           formula = ht ~ site*seedlot
         , data    = DataExam5.1
         )

 # Pg. 73
 anova(fm5.4)
Analysis of Variance Table

Response: ht
             Df  Sum Sq Mean Sq F value Pr(>F)
site          3  919585  306528     NaN    NaN
seedlot      26 3176289  122165     NaN    NaN
site:seedlot 78  707957    9076     NaN    NaN
Residuals     0       0     NaN               

 # Pg. 73
 emmeans(object = fm5.4, specs = ~ site)
 site       emmean  SE df lower.CL upper.CL
 Ratchaburi    462 NaN  0      NaN      NaN
 Sai Thong     628 NaN  0      NaN      NaN
 Si Sa Ket     494 NaN  0      NaN      NaN
 Sakaerat      370 NaN  0      NaN      NaN

Results are averaged over the levels of: seedlot 
Confidence level used: 0.95 

 emmeans(object = fm5.4, specs = ~ seedlot)
 seedlot emmean  SE df lower.CL upper.CL
 13877      365 NaN  0      NaN      NaN
 13866      353 NaN  0      NaN      NaN
 13689      559 NaN  0      NaN      NaN
 13688      546 NaN  0      NaN      NaN
 13861      627 NaN  0      NaN      NaN
 13854      628 NaN  0      NaN      NaN
 13684      660 NaN  0      NaN      NaN
 13864      422 NaN  0      NaN      NaN
 13863      586 NaN  0      NaN      NaN
 13683      770 NaN  0      NaN      NaN
 13681      695 NaN  0      NaN      NaN
 14175      438 NaN  0      NaN      NaN
 14660      521 NaN  0      NaN      NaN
 13653      592 NaN  0      NaN      NaN
 13846      440 NaN  0      NaN      NaN
 13621      384 NaN  0      NaN      NaN
 13871      272 NaN  0      NaN      NaN
 13519      422 NaN  0      NaN      NaN
 13514      369 NaN  0      NaN      NaN
 13148      273 NaN  0      NaN      NaN
 13990      282 NaN  0      NaN      NaN
 14537      780 NaN  0      NaN      NaN
 14106      772 NaN  0      NaN      NaN
 12013      616 NaN  0      NaN      NaN
 14130      422 NaN  0      NaN      NaN
 14485      123 NaN  0      NaN      NaN
 11935      273 NaN  0      NaN      NaN

Results are averaged over the levels of: site 
Confidence level used: 0.95 

 ANOVAfm5.4 <- anova(fm5.4)

 ANOVAfm5.4[4, 1:3] <- c(208, 208*1040, 1040)

 ANOVAfm5.4[3, 4]   <- ANOVAfm5.4[3, 3]/ANOVAfm5.4[4, 3]

 ANOVAfm5.4[3, 5]   <-
            pf(
               q        = ANOVAfm5.4[3, 4]
           , df1        = ANOVAfm5.4[3, 1]
           , df2        = ANOVAfm5.4[4, 1]
           , lower.tail = FALSE
           )

 # Pg. 73
 ANOVAfm5.4
Analysis of Variance Table

Response: ht
              Df  Sum Sq Mean Sq F value                Pr(>F)    
site           3  919585  306528     NaN                   NaN    
seedlot       26 3176289  122165     NaN                   NaN    
site:seedlot  78  707957    9076  8.7273 < 0.00000000000000022 ***
Residuals    208  216320    1040                                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 # Pg. 80
 DataExam5.1 %>%
   filter(seedlot %in% c("13653", "13871")) %>%
   ggplot(
     data = .
   , mapping = aes(
                   x     = sitemean
                 , y     = ht
                 , color = seedlot
                 , shape = seedlot
                 )
   ) +
   geom_point() +
   geom_smooth(
      method    = lm
    , se        = FALSE
    , fullrange = TRUE
    ) +
   theme_classic() +
   labs(
       x = "SiteMean"
     , y = "SeedLot Mean"
     )

 Tab5.10 <-
   DataExam5.1 %>%
   summarise(Mean = mean(ht), .by = seedlot) %>%
   left_join(
      DataExam5.1 %>%
      nest_by(seedlot) %>%
      mutate(fm1 = list(lm(ht ~ sitemean, data = data))) %>%
      summarise(Slope = coef(fm1)[2])
   , by = "seedlot"
      )

 # Pg. 81
 Tab5.10
   seedlot   Mean      Slope
1    11935 272.75 0.53017435
2    14485 123.00 0.10170020
3    14130 422.25 0.54976906
4    12013 616.25 2.06723798
5    14106 771.75 1.37751724
6    14537 779.75 0.96012145
7    13990 281.75 0.08298796
8    13148 273.25 0.05333546
9    13514 368.75 0.12307233
10   13519 422.00 0.29211648
11   13871 271.50 0.83048203
12   13621 383.75 1.20085607
13   13846 440.00 0.70691001
14   13653 591.50 1.59434380
15   14660 521.25 0.93353990
16   14175 438.00 1.33770745
17   13681 695.00 1.30937837
18   13683 769.75 1.79629735
19   13863 586.25 1.48034730
20   13864 422.50 0.61113857
21   13684 660.00 1.67860570
22   13854 628.00 1.62026853
23   13861 626.75 1.43784662
24   13688 546.50 1.72717652
25   13689 558.75 1.19475332
26   13866 352.75 0.61009734
27   13877 364.75 0.79221858

 ggplot(data = Tab5.10, mapping = aes(x = Mean, y = Slope)) +
  geom_point(size = 2) +
  theme_bw() +
  labs(
      x = "SeedLot Mean"
    , y = "Regression Coefficient"
    )

 DevSS1 <-
   DataExam5.1 %>%
   nest_by(seedlot) %>%
   mutate(fm1 = list(lm(ht ~ sitemean, data = data))) %>%
   summarise(SSE = anova(fm1)[2, 2]) %>%
   ungroup() %>%
   summarise(Dev = sum(SSE)) %>%
   as.numeric()

 ANOVAfm5.4[2, 2]
[1] 3176289

 length(levels(DataExam5.1$SeedLot))
[1] 0

 ANOVAfm5.4.1 <-
   rbind(
    ANOVAfm5.4[1:3, ]
   , c(
       ANOVAfm5.4[2, 1]
     , ANOVAfm5.4[3, 2] - DevSS1
     , (ANOVAfm5.4[3, 2] - DevSS1)/ANOVAfm5.4[2, 1]
     , NA
     , NA
     )
   , c(
       ANOVAfm5.4[3, 1]-ANOVAfm5.4[2, 1]
     , DevSS1
     , DevSS1/(ANOVAfm5.4[3, 1]-ANOVAfm5.4[2, 1])
     , DevSS1/(ANOVAfm5.4[3, 1]-ANOVAfm5.4[2, 1])/ANOVAfm5.4[4, 3]
     , pf(
             q = DevSS1/(ANOVAfm5.4[3, 1]-ANOVAfm5.4[2, 1])/ANOVAfm5.4[4, 3]
         , df1 = ANOVAfm5.4[3, 1]-ANOVAfm5.4[2, 1]
         , df2 = ANOVAfm5.4[4, 1]
         , lower.tail = FALSE
         )
     )
   , ANOVAfm5.4[4, ]
   )

 rownames(ANOVAfm5.4.1) <-
   c(
     "Site"
   , "seedlot"
   , "site:seedlot"
   , "  regressions"
   , "  deviations"
   , "Residuals"
   )

 # Pg. 82
 ANOVAfm5.4.1
Analysis of Variance Table

Response: ht
               Df  Sum Sq Mean Sq F value                Pr(>F)    
Site            3  919585  306528     NaN                   NaN    
seedlot        26 3176289  122165     NaN                   NaN    
site:seedlot   78  707957    9076  8.7273 < 0.00000000000000022 ***
  regressions  26  308503   11866                                  
  deviations   52  399454    7682  7.3863 < 0.00000000000000022 ***
Residuals     208  216320    1040                                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2 Example 5.2 (Pg. 72)

Example 5.2 (Pg. 72)

 library(car)

 library(dae)

 library(dplyr)

 library(emmeans)

 library(ggplot2)

 library(lmerTest)

 library(magrittr)

 library(predictmeans)

 data(DataExam5.2)

 # Pg. 75
 fm5.7 <-
  lm(
      formula = ht ~ site*seedlot
    , data    = DataExam5.2
    )

 # Pg. 77
 anova(fm5.7)
Analysis of Variance Table

Response: ht
              Df  Sum Sq Mean Sq F value Pr(>F)
site           5 4157543  831509     NaN    NaN
seedlot       36 4425296  122925     NaN    NaN
site:seedlot 150 1351054    9007     NaN    NaN
Residuals      0       0     NaN               

 fm5.9 <-
  lm(
      formula = ht ~ site*seedlot
    , data    = DataExam5.2
    )

 # Pg. 77
 anova(fm5.9)
Analysis of Variance Table

Response: ht
              Df  Sum Sq Mean Sq F value Pr(>F)
site           5 4157543  831509     NaN    NaN
seedlot       36 4425296  122925     NaN    NaN
site:seedlot 150 1351054    9007     NaN    NaN
Residuals      0       0     NaN               

 ANOVAfm5.9 <- anova(fm5.9)

 ANOVAfm5.9[4, 1:3] <- c(384, 384*964, 964)

 ANOVAfm5.9[3, 4]   <- ANOVAfm5.9[3, 3]/ANOVAfm5.9[4, 3]

 ANOVAfm5.9[3, 5]   <-
     pf(
         q = ANOVAfm5.9[3, 4]
     , df1 = ANOVAfm5.9[3, 1]
     , df2 = ANOVAfm5.9[4, 1]
     , lower.tail = FALSE
     )

 # Pg. 77
 ANOVAfm5.9
Analysis of Variance Table

Response: ht
              Df  Sum Sq Mean Sq F value                Pr(>F)    
site           5 4157543  831509     NaN                   NaN    
seedlot       36 4425296  122925     NaN                   NaN    
site:seedlot 150 1351054    9007  9.3434 < 0.00000000000000022 ***
Residuals    384  370176     964                                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Tab5.14 <-
   DataExam5.2 %>%
   summarise(
        Mean = round(mean(ht, na.rm = TRUE), 0)
      , .by  = seedlot
      ) %>%
    left_join(
       DataExam5.2 %>%
       nest_by(seedlot) %>%
       mutate(fm2 = list(lm(ht ~ sitemean, data = data))) %>%
       summarise(Slope = round(coef(fm2)[2], 2))
     , by = "seedlot"
      ) %>%
   as.data.frame()

 # Pg. 81
 Tab5.14
   seedlot Mean Slope
1    13877  291  0.80
2    13866  302  0.56
3    13689  463  1.07
4    13688  458  1.20
5    13861  538  1.08
6    13854  536  1.16
7    13686  726  1.77
8    13684  549  1.34
9    13864  371  0.79
10   13863  586  1.19
11   13683  673  1.78
12   13681  610  1.33
13   13680  645  2.08
14   14623  652  0.49
15   14175  378  1.22
16   14660  445  0.88
17   13691  511  1.61
18   13653  477  1.35
19   13846  381  0.93
20   13621  304  0.98
21   14176  177  0.22
22   13871  220  0.64
23   14622  355  1.65
24   13876  240  0.77
25   13519  342  0.71
26   13514  290  0.66
27   13148  230  0.35
28   13990  260  0.29
29   14537  671  1.15
30   14106  673  1.18
31   12013  529  1.30
32   14130  390  0.45
33   14485  108  0.16
34   14166  210  0.61
35   11935  211  0.61
36   14170  244  0.59
37   14152  150  0.63

 DevSS2 <-
   DataExam5.2 %>%
   nest_by(seedlot) %>%
   mutate(fm2 = list(lm(ht ~ sitemean, data = data))) %>%
   summarise(SSE = anova(fm2)[2, 2]) %>%
   ungroup() %>%
   summarise(Dev = sum(SSE)) %>%
   as.numeric()

 ANOVAfm5.9.1 <-
   rbind(
      ANOVAfm5.9[1:3, ]
    , c(
         ANOVAfm5.9[2, 1]
       , ANOVAfm5.9[3, 2] - DevSS2
       , (ANOVAfm5.9[3, 2] - DevSS2)/ANOVAfm5.9[2, 1]
       , NA
       , NA
       )
    , c(
         ANOVAfm5.9[3, 1]-ANOVAfm5.9[2, 1]
       , DevSS2
       , DevSS2/(ANOVAfm5.9[3, 1]-ANOVAfm5.9[2, 1])
       , DevSS2/(ANOVAfm5.9[3, 1]-ANOVAfm5.9[2, 1])/ANOVAfm5.9[4, 3]
       , pf(
               q = DevSS2/(ANOVAfm5.9[3, 1]-ANOVAfm5.9[2, 1])/ANOVAfm5.9[4, 3]
           , df1 = ANOVAfm5.9[3, 1]-ANOVAfm5.9[2, 1]
           , df2 = ANOVAfm5.9[4, 1]
           , lower.tail = FALSE
           )
       )
    , ANOVAfm5.9[4, ]
   )

 rownames(ANOVAfm5.9.1) <-
   c(
      "site"
    , "seedlot"
    , "site:seedlot"
    , "  regressions"
    , "  deviations"
    , "Residuals"
    )

 # Pg. 82
 ANOVAfm5.9.1
Analysis of Variance Table

Response: ht
               Df  Sum Sq Mean Sq F value                Pr(>F)    
site            5 4157543  831509     NaN                   NaN    
seedlot        36 4425296  122925     NaN                   NaN    
site:seedlot  150 1351054    9007  9.3434 < 0.00000000000000022 ***
  regressions  36  703203   19533                                  
  deviations  114  647851    5683  5.8951 < 0.00000000000000022 ***
Residuals     384  370176     964                                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Code <-
  c(
    "a","a","a","a","b","b","b","b"
  , "c","d","d","d","d","e","f","g"
  , "h","h","i","i","j","k","l","m"
  ,"n","n","n","o","p","p","q","r"
  , "s","t","t","u","v"
  )

 Tab5.14$Code <- Code

 ggplot(
    data = Tab5.14
  , mapping = aes(x = Mean, y = Slope)
  ) +
  geom_point(size = 2) +
  geom_text(
     mapping = aes(label = Code)
   , hjust   = -0.5
   , vjust   = -0.5
  ) +
  theme_bw() +
  labs(
      x = "SeedLot Mean"
    , y = "Regression Coefficient"
    )