Exam5.1: Example 5.1 from Experimental Design and Analysis for Tree...

In eda4treeR: Experimental Design and Analysis for Tree Improvement

Example 5.1 from Experimental Design and Analysis for Tree Improvement

Description

Exam5.1 presents the height of 27 seedlots from 4 sites.

Author(s)

1. Muhammad Yaseen ([e-mail])

2. Sami Ullah ([e-mail])

References

1. E.R. Williams, C.E. Harwood and A.C. Matheson (2023). Experimental Design and Analysis for Tree Improvement. CSIRO Publishing (https://www.publish.csiro.au/book/3145/).

See Also

DataExam5.1

Examples

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)
# Pg. 73
emmeans(object = fm5.4, specs = ~ site)
emmeans(object = fm5.4, specs = ~ seedlot)

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

# 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

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]

length(levels(DataExam5.1$SeedLot))

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

eda4treeR documentation built on Sept. 14, 2024, 1:08 a.m.