Skip to contents

Example 6.2 from Experimental Design and Analysis for Tree Improvement

Source: R/Exam6.2.R
Exam6.2.Rd

Exam 6.2 Dbh mean, Dbh varince and number of trees per plot from 3 provinces("PNG","Sabah","Queensland") with 4 replications of 48 families.

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

DataExam6.2

Author

  1. Muhammad Yaseen (myaseen208@gmail.com)

  2. Sami Ullah (samiullahuos@gmail.com)

Examples

library(car)
library(dae)
library(dplyr)
library(emmeans)
library(ggplot2)
library(lmerTest)
library(magrittr)
library(predictmeans)

data(DataExam6.2)

DataExam6.2.1 <-
    DataExam6.2 %>%
    filter(Province == "PNG")

# Pg. 94
fm6.3 <-
     lm(
          formula = Dbh.mean ~ Replication + Family
        , data    = DataExam6.2.1
       )

b    <- anova(fm6.3)


HM      <- function(x){length(x)/sum(1/x)}
w       <- HM(DataExam6.2.1$Dbh.count)
S2      <- b[["Mean Sq"]][length(b[["Mean Sq"]])]
Sigma2t <- mean(DataExam6.2.1$Dbh.variance)
sigma2m <- S2-(Sigma2t/w)

fm6.3.1 <-
  lmer(
      formula = Dbh.mean ~ 1 + Replication + (1|Family)
    , data    = DataExam6.2.1
    , REML    = TRUE
    )

# Pg. 104
# summary(fm6.3.1)
varcomp(fm6.3.1)
#> Computing profile confidence intervals ...
#>                      vcov     SE  2.5 % 97.5 %
#> Family.(Intercept) 0.2584 0.1286 0.0538 0.5767
#> residual           1.1667 0.1506 0.8954 1.4774
sigma2f <- 0.2584
h2 <- (sigma2f/(0.3))/(Sigma2t + sigma2m + sigma2f)
cbind(hmean = w, Sigma2t, sigma2m, sigma2f, h2)
#>         hmean  Sigma2t   sigma2m sigma2f        h2
#> [1,] 4.408602 3.920732 0.2773606  0.2584 0.1932761

fm6.4 <-
  lm(
      formula = Dbh.mean ~ Replication+Family
     , data   = DataExam6.2
     )

b    <- anova(fm6.4)
HM      <- function(x){length(x)/sum(1/x)}
w       <- HM(DataExam6.2$Dbh.count)
S2      <- b[["Mean Sq"]][length(b[["Mean Sq"]])]
Sigma2t <- mean(DataExam6.2$Dbh.variance)
sigma2m <- S2-(Sigma2t/w)

fm6.4.1 <-
 lmer(
   formula = Dbh.mean ~ 1 + Replication + Province + (1|Family)
 , data    = DataExam6.2
 , REML    = TRUE
    )

# Pg. 107
varcomp(fm6.4.1)
#> Computing profile confidence intervals ...
#>                      vcov     SE  2.5 % 97.5 %
#> Family.(Intercept) 0.3514 0.1358 0.1203 0.6361
#> residual           1.0951 0.1304 0.8584 1.3634
sigma2f <- 0.3514
h2 <- (sigma2f/(0.3))/(Sigma2t+sigma2m+sigma2f)
cbind(hmean = w, Sigma2t, sigma2m, sigma2f, h2)
#>         hmean  Sigma2t  sigma2m sigma2f        h2
#> [1,] 4.451314 3.860156 0.227873  0.3514 0.2638477

fm6.7.1 <-
 lmer(
   formula = Dbh.mean ~ 1+Replication+(1|Family)
 , data    = DataExam6.2.1
 , REML = TRUE
 )

# Pg. 116
varcomp(fm6.7.1)
#> Computing profile confidence intervals ...
#>                      vcov     SE  2.5 % 97.5 %
#> Family.(Intercept) 0.2584 0.1286 0.0538 0.5767
#> residual           1.1667 0.1506 0.8954 1.4774
sigma2f[1] <- 0.2584

fm6.7.2<-
 lmer(
   formula = Ht.mean ~ 1 + Replication + (1|Family)
 , data    = DataExam6.2.1
 , REML    = TRUE
   )

# Pg. 116
varcomp(fm6.7.2)
#> Computing profile confidence intervals ...
#>                      vcov     SE  2.5 % 97.5 %
#> Family.(Intercept) 0.2711 0.1243 0.0743 0.5794
#> residual           1.0549 0.1362 0.8097 1.3359
sigma2f[2] <- 0.2711

fm6.7.3 <-
 lmer(
   formula = Sum.means ~ 1 + Replication + (1|Family)
 , data    = DataExam6.2.1
 , REML    = TRUE
 , control = lmerControl()
 )

# Pg. 116
varcomp(fm6.7.3)
#> Computing profile confidence intervals ...
#>                      vcov     SE  2.5 % 97.5 %
#> Family.(Intercept) 0.8729 0.3907 0.2553 1.8421
#> residual           3.2428 0.4186 2.4888 4.1063
sigma2f[3] <- 0.873
sigma2xy   <- 0.5*(sigma2f[3]-sigma2f[1]-sigma2f[2])
GenCorr <- sigma2xy/sqrt(sigma2f[1]*sigma2f[2])
cbind(
     S2x = sigma2f[1]
   , S2y = sigma2f[2]
   , S2.x.plus.y = sigma2f[3]
   , GenCorr
   )
#>         S2x    S2y S2.x.plus.y   GenCorr
#> [1,] 0.2584 0.2711       0.873 0.6489119