Example 8.1.1 from Experimental Design and Analysis for Tree Improvement
Source:R/Exam8.1.1.R
Exam8.1.1.RdExam8.1.1 presents the Mixed Effects Analysis of Diameter at breast height (Dbh) of 60 SeedLots under layout of row column design with 6 rows and 10 columns in 18 countries and 59 provinces of 18 selected countries given in Example 8.1.
References
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/).
Author
Muhammad Yaseen (myaseen208@gmail.com)
Sami Ullah (samiullahuos@gmail.com)
Examples
library(car)
library(dae)
library(dplyr)
library(emmeans)
library(ggplot2)
library(lmerTest)
library(magrittr)
library(predictmeans)
data(DataExam8.1)
# Pg. 155
fm8.8 <-
lmerTest::lmer(
formula = dbh ~ 1 + repl + col + prov +
(1|repl:row) + (1|repl:col)
, data = DataExam8.1
, REML = TRUE
)
# Pg. 157
if (FALSE) { # \dontrun{
varcomp(fm8.8)
} # }
anova(fm8.8)
#> Type III Analysis of Variance Table with Satterthwaite's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> repl 2.581 0.86023 3 21.257 4.4082 0.01469 *
#> col 24.874 2.76378 9 23.705 14.1627 1.511e-07 ***
#> prov 55.433 0.95574 58 136.623 4.8976 1.306e-14 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(fm8.8, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> repl 2.580 0.86016 3 22.622 4.4078 0.0139 *
#> col 24.824 2.75827 9 22.947 14.1337 2.098e-07 ***
#> prov 54.795 0.94473 58 133.852 4.8396 2.830e-14 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predictmeans(model = fm8.8, modelterm = "repl")
#> $`Predicted Means`
#> repl
#> 1 2 3 4
#> 3.7543 3.1269 3.2128 3.5023
#>
#> $`Standard Error of Means`
#> All means have the same SE
#> 0.13629
#>
#> $`Standard Error of Differences`
#> Max.SED Min.SED Aveg.SED
#> 0.1927495 0.1927423 0.1927472
#>
#> $LSD
#> Max.LSD Min.LSD Aveg.LSD
#> 0.39910 0.39909 0.39910
#> attr(,"Significant level")
#> [1] 0.05
#> attr(,"Degree of freedom")
#> [1] 22.62
#>
#> $mean_table
#> repl Mean SE Df LL(95%) UL(95%)
#> 1 1 3.7543 0.13629 22.6214 3.4721 4.0365
#> 2 2 3.1269 0.13629 22.6214 2.8447 3.4091
#> 3 3 3.2128 0.13629 22.6214 2.9306 3.4950
#> 4 4 3.5023 0.13629 22.6214 3.2201 3.7845
#>
predictmeans(model = fm8.8, modelterm = "col")
#> $`Predicted Means`
#> col
#> 1 2 3 4 5 6 7 8 9 10
#> 3.5053 3.4996 3.8509 3.8280 3.5947 3.7829 3.3059 3.5158 3.1776 1.9301
#>
#> $`Standard Error of Means`
#> col
#> 1 2 3 4 5 6 7 8 9 10
#> 0.15496 0.15649 0.15493 0.15494 0.15490 0.15721 0.15479 0.15699 0.15485 0.15648
#>
#> $`Standard Error of Differences`
#> Max.SED Min.SED Aveg.SED
#> 0.2132295 0.2061223 0.2086161
#>
#> $LSD
#> Max.LSD Min.LSD Aveg.LSD
#> 0.43742 0.42284 0.42796
#> attr(,"Significant level")
#> [1] 0.05
#> attr(,"Degree of freedom")
#> [1] 27.12
#>
#> $mean_table
#> col Mean SE Df LL(95%) UL(95%)
#> 1 1 3.5053 0.15496 27.11736 3.1874 3.8232
#> 2 2 3.4996 0.15649 27.11736 3.1786 3.8207
#> 3 3 3.8509 0.15493 27.11736 3.5330 4.1687
#> 4 4 3.8280 0.15494 27.11736 3.5101 4.1458
#> 5 5 3.5947 0.15490 27.11736 3.2769 3.9125
#> 6 6 3.7829 0.15721 27.11736 3.4604 4.1054
#> 7 7 3.3059 0.15479 27.11736 2.9883 3.6234
#> 8 8 3.5158 0.15699 27.11736 3.1938 3.8379
#> 9 9 3.1776 0.15485 27.11736 2.8599 3.4952
#> 10 10 1.9301 0.15648 27.11736 1.6091 2.2511
#>
predictmeans(model = fm8.8, modelterm = "prov")
#> $`Predicted Means`
#> prov
#> 1 2 3 4 5 6 7 8 10 11 12
#> 2.4222 3.1425 2.8646 2.2599 3.7189 3.5445 3.9474 2.7469 2.6459 2.0497 2.7497
#> 13 14 15 16 17 18 19 20 21 22 23
#> 2.9833 2.6534 3.4459 3.9209 3.4294 3.3624 3.6540 3.3685 3.2210 3.3315 3.2453
#> 24 25 26 27 28 29 30 31 32 33 34
#> 3.6933 3.8460 3.5233 3.5271 3.1001 3.8997 3.5447 4.1294 4.2755 3.7429 3.8890
#> 35 36 37 38 39 40 41 42 45 46 47
#> 4.1179 3.8526 3.4538 3.0379 3.5306 3.7412 3.7422 3.8765 4.2126 3.4660 4.4688
#> 48 50 51 52 53 54 55 56 57 58 59
#> 3.7503 2.5571 3.4571 3.3227 3.5926 3.5337 2.7157 2.7702 3.9805 2.9968 3.3978
#> 60 61 62 63
#> 3.0115 3.2238 3.2814 3.5751
#>
#> $`Standard Error of Means`
#> prov
#> 1 2 3 4 5 6 7 8 10 11
#> 0.25362 0.25446 0.25278 0.25443 0.25423 0.25588 0.25420 0.25338 0.25417 0.25259
#> 12 13 14 15 16 17 18 19 20 21
#> 0.25306 0.25436 0.25259 0.25232 0.25502 0.25319 0.25340 0.25309 0.25294 0.25204
#> 22 23 24 25 26 27 28 29 30 31
#> 0.25350 0.25355 0.25358 0.25316 0.25289 0.25294 0.25291 0.25341 0.25254 0.25330
#> 32 33 34 35 36 37 38 39 40 41
#> 0.25377 0.25387 0.25510 0.25347 0.25351 0.25352 0.25416 0.25393 0.25450 0.25416
#> 42 45 46 47 48 50 51 52 53 54
#> 0.25303 0.25257 0.25323 0.25316 0.25356 0.25354 0.25372 0.25361 0.25237 0.25321
#> 55 56 57 58 59 60 61 62 63
#> 0.25374 0.25639 0.25173 0.25319 0.25250 0.25303 0.25378 0.25266 0.25500
#>
#> $`Standard Error of Differences`
#> Max.SED Min.SED Aveg.SED
#> 0.3591712 0.3326719 0.3482665
#>
#> $LSD
#> Max.LSD Min.LSD Aveg.LSD
#> 0.70924 0.65691 0.68770
#> attr(,"Significant level")
#> [1] 0.05
#> attr(,"Degree of freedom")
#> [1] 162.73
#>
#> $mean_table
#> prov Mean SE Df LL(95%) UL(95%)
#> 1 1 2.4222 0.25362 162.7269 1.9214 2.9230
#> 2 2 3.1425 0.25446 162.7269 2.6401 3.6450
#> 3 3 2.8646 0.25278 162.7269 2.3655 3.3638
#> 4 4 2.2599 0.25443 162.7269 1.7575 2.7623
#> 5 5 3.7189 0.25423 162.7269 3.2169 4.2209
#> 6 6 3.5445 0.25588 162.7269 3.0392 4.0498
#> 7 7 3.9474 0.25420 162.7269 3.4454 4.4493
#> 8 8 2.7469 0.25338 162.7269 2.2466 3.2473
#> 9 10 2.6459 0.25417 162.7269 2.1440 3.1478
#> 10 11 2.0497 0.25259 162.7269 1.5509 2.5485
#> 11 12 2.7497 0.25306 162.7269 2.2500 3.2494
#> 12 13 2.9833 0.25436 162.7269 2.4810 3.4855
#> 13 14 2.6534 0.25259 162.7269 2.1546 3.1522
#> 14 15 3.4459 0.25232 162.7269 2.9476 3.9441
#> 15 16 3.9209 0.25502 162.7269 3.4173 4.4245
#> 16 17 3.4294 0.25319 162.7269 2.9294 3.9293
#> 17 18 3.3624 0.25340 162.7269 2.8620 3.8628
#> 18 19 3.6540 0.25309 162.7269 3.1543 4.1538
#> 19 20 3.3685 0.25294 162.7269 2.8690 3.8680
#> 20 21 3.2210 0.25204 162.7269 2.7234 3.7187
#> 21 22 3.3315 0.25350 162.7269 2.8310 3.8321
#> 22 23 3.2453 0.25355 162.7269 2.7447 3.7460
#> 23 24 3.6933 0.25358 162.7269 3.1926 4.1941
#> 24 25 3.8460 0.25316 162.7269 3.3461 4.3459
#> 25 26 3.5233 0.25289 162.7269 3.0239 4.0227
#> 26 27 3.5271 0.25294 162.7269 3.0276 4.0265
#> 27 28 3.1001 0.25291 162.7269 2.6007 3.5996
#> 28 29 3.8997 0.25341 162.7269 3.3993 4.4001
#> 29 30 3.5447 0.25254 162.7269 3.0460 4.0434
#> 30 31 4.1294 0.25330 162.7269 3.6293 4.6296
#> 31 32 4.2755 0.25377 162.7269 3.7744 4.7766
#> 32 33 3.7429 0.25387 162.7269 3.2416 4.2442
#> 33 34 3.8890 0.25510 162.7269 3.3853 4.3928
#> 34 35 4.1179 0.25347 162.7269 3.6174 4.6185
#> 35 36 3.8526 0.25351 162.7269 3.3520 4.3532
#> 36 37 3.4538 0.25352 162.7269 2.9532 3.9544
#> 37 38 3.0379 0.25416 162.7269 2.5360 3.5398
#> 38 39 3.5306 0.25393 162.7269 3.0291 4.0320
#> 39 40 3.7412 0.25450 162.7269 3.2386 4.2437
#> 40 41 3.7422 0.25416 162.7269 3.2403 4.2440
#> 41 42 3.8765 0.25303 162.7269 3.3769 4.3761
#> 42 45 4.2126 0.25257 162.7269 3.7139 4.7114
#> 43 46 3.4660 0.25323 162.7269 2.9660 3.9661
#> 44 47 4.4688 0.25316 162.7269 3.9689 4.9687
#> 45 48 3.7503 0.25356 162.7269 3.2496 4.2510
#> 46 50 2.5571 0.25354 162.7269 2.0564 3.0577
#> 47 51 3.4571 0.25372 162.7269 2.9561 3.9582
#> 48 52 3.3227 0.25361 162.7269 2.8219 3.8235
#> 49 53 3.5926 0.25237 162.7269 3.0943 4.0910
#> 50 54 3.5337 0.25321 162.7269 3.0337 4.0337
#> 51 55 2.7157 0.25374 162.7269 2.2146 3.2167
#> 52 56 2.7702 0.25639 162.7269 2.2640 3.2765
#> 53 57 3.9805 0.25173 162.7269 3.4834 4.4776
#> 54 58 2.9968 0.25319 162.7269 2.4968 3.4968
#> 55 59 3.3978 0.25250 162.7269 2.8992 3.8963
#> 56 60 3.0115 0.25303 162.7269 2.5118 3.5112
#> 57 61 3.2238 0.25378 162.7269 2.7227 3.7250
#> 58 62 3.2814 0.25266 162.7269 2.7825 3.7803
#> 59 63 3.5751 0.25500 162.7269 3.0716 4.0787
#>
# Pg. 161
RCB1 <-
aov(dbh ~ prov + repl, data = DataExam8.1)
RCB <-
emmeans(RCB1, specs = "prov") %>%
as_tibble()
Mixed <-
emmeans(fm8.8, specs = "prov") %>%
as_tibble()
table8.9 <-
left_join(
x = RCB
, y = Mixed
, by = "prov"
, suffix = c(".RCBD", ".Mixed")
)
print(table8.9)
#> # A tibble: 59 × 11
#> prov emmean.RCBD SE.RCBD df.RCBD lower.CL.RCBD upper.CL.RCBD emmean.Mixed
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1.85 0.382 174 1.10 2.60 2.42
#> 2 2 3.24 0.382 174 2.49 3.99 3.14
#> 3 3 2.94 0.382 174 2.18 3.69 2.86
#> 4 4 2.50 0.382 174 1.74 3.25 2.26
#> 5 5 3.34 0.382 174 2.59 4.10 3.72
#> 6 6 3.37 0.382 174 2.62 4.13 3.54
#> 7 7 3.98 0.382 174 3.23 4.74 3.95
#> 8 8 2.63 0.382 174 1.88 3.39 2.75
#> 9 10 2.47 0.382 174 1.71 3.22 2.65
#> 10 11 2.40 0.382 174 1.64 3.15 2.05
#> # ℹ 49 more rows
#> # ℹ 4 more variables: SE.Mixed <dbl>, df.Mixed <dbl>, lower.CL.Mixed <dbl>,
#> # upper.CL.Mixed <dbl>