Example 8.1.1 from Experimental Design and Analysis for Tree Improvement

Exam8.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)
library(supernova)

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
# \dontrun{
varcomp(fm8.8)
#> Computing profile confidence intervals ...
#>                        vcov     SE  2.5 % 97.5 %
#> repl:col.(Intercept) 0.0459 0.0262 0.0000 0.0565
#> repl:row.(Intercept) 0.0640 0.0294 0.0210 0.1161
#> residual             0.1951 0.0253 0.1126 0.1782
# }

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.62064  3.4721  4.0365
#> 2    2 3.1269 0.13629 22.62064  2.8447  3.4091
#> 3    3 3.2128 0.13629 22.62064  2.9306  3.4950
#> 4    4 3.5023 0.13629 22.62064  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.43749  0.42291  0.42802 
#> attr(,"Significant level")
#> [1] 0.05
#> attr(,"Degree of freedom")
#> [1] 27.03
#> 
#> $mean_table
#>    col   Mean      SE       Df LL(95%) UL(95%)
#> 1    1 3.5053 0.15496 27.03111  3.1873  3.8232
#> 2    2 3.4996 0.15649 27.03111  3.1786  3.8207
#> 3    3 3.8509 0.15493 27.03111  3.5330  4.1687
#> 4    4 3.8280 0.15494 27.03111  3.5101  4.1459
#> 5    5 3.5947 0.15490 27.03111  3.2769  3.9125
#> 6    6 3.7829 0.15721 27.03111  3.4604  4.1055
#> 7    7 3.3059 0.15479 27.03111  2.9883  3.6234
#> 8    8 3.5158 0.15699 27.03111  3.1937  3.8379
#> 9    9 3.1776 0.15485 27.03111  2.8599  3.4953
#> 10  10 1.9301 0.15648 27.03111  1.6091  2.2512
#> 
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.70929  0.65696  0.68775 
#> attr(,"Significant level")
#> [1] 0.05
#> attr(,"Degree of freedom")
#> [1] 161.14
#> 
#> $mean_table
#>    prov   Mean      SE       Df LL(95%) UL(95%)
#> 1     1 2.4222 0.25362 161.1392  1.9214  2.9231
#> 2     2 3.1425 0.25446 161.1392  2.6400  3.6450
#> 3     3 2.8646 0.25278 161.1392  2.3654  3.3638
#> 4     4 2.2599 0.25443 161.1392  1.7574  2.7623
#> 5     5 3.7189 0.25423 161.1392  3.2168  4.2209
#> 6     6 3.5445 0.25588 161.1392  3.0392  4.0498
#> 7     7 3.9474 0.25420 161.1392  3.4454  4.4494
#> 8     8 2.7469 0.25338 161.1392  2.2466  3.2473
#> 9    10 2.6459 0.25417 161.1392  2.1440  3.1478
#> 10   11 2.0497 0.25259 161.1392  1.5509  2.5485
#> 11   12 2.7497 0.25306 161.1392  2.2499  3.2494
#> 12   13 2.9833 0.25436 161.1392  2.4809  3.4856
#> 13   14 2.6534 0.25259 161.1392  2.1546  3.1522
#> 14   15 3.4459 0.25232 161.1392  2.9476  3.9442
#> 15   16 3.9209 0.25502 161.1392  3.4173  4.4245
#> 16   17 3.4294 0.25319 161.1392  2.9294  3.9294
#> 17   18 3.3624 0.25340 161.1392  2.8620  3.8628
#> 18   19 3.6540 0.25309 161.1392  3.1542  4.1538
#> 19   20 3.3685 0.25294 161.1392  2.8690  3.8680
#> 20   21 3.2210 0.25204 161.1392  2.7233  3.7188
#> 21   22 3.3315 0.25350 161.1392  2.8309  3.8321
#> 22   23 3.2453 0.25355 161.1392  2.7446  3.7461
#> 23   24 3.6933 0.25358 161.1392  3.1926  4.1941
#> 24   25 3.8460 0.25316 161.1392  3.3460  4.3459
#> 25   26 3.5233 0.25289 161.1392  3.0239  4.0227
#> 26   27 3.5271 0.25294 161.1392  3.0276  4.0266
#> 27   28 3.1001 0.25291 161.1392  2.6007  3.5996
#> 28   29 3.8997 0.25341 161.1392  3.3993  4.4001
#> 29   30 3.5447 0.25254 161.1392  3.0460  4.0434
#> 30   31 4.1294 0.25330 161.1392  3.6292  4.6297
#> 31   32 4.2755 0.25377 161.1392  3.7744  4.7767
#> 32   33 3.7429 0.25387 161.1392  3.2415  4.2442
#> 33   34 3.8890 0.25510 161.1392  3.3853  4.3928
#> 34   35 4.1179 0.25347 161.1392  3.6174  4.6185
#> 35   36 3.8526 0.25351 161.1392  3.3520  4.3533
#> 36   37 3.4538 0.25352 161.1392  2.9531  3.9544
#> 37   38 3.0379 0.25416 161.1392  2.5360  3.5398
#> 38   39 3.5306 0.25393 161.1392  3.0291  4.0320
#> 39   40 3.7412 0.25450 161.1392  3.2386  4.2437
#> 40   41 3.7422 0.25416 161.1392  3.2403  4.2441
#> 41   42 3.8765 0.25303 161.1392  3.3768  4.3762
#> 42   45 4.2126 0.25257 161.1392  3.7138  4.7114
#> 43   46 3.4660 0.25323 161.1392  2.9660  3.9661
#> 44   47 4.4688 0.25316 161.1392  3.9688  4.9687
#> 45   48 3.7503 0.25356 161.1392  3.2496  4.2510
#> 46   50 2.5571 0.25354 161.1392  2.0564  3.0578
#> 47   51 3.4571 0.25372 161.1392  2.9561  3.9582
#> 48   52 3.3227 0.25361 161.1392  2.8219  3.8236
#> 49   53 3.5926 0.25237 161.1392  3.0942  4.0910
#> 50   54 3.5337 0.25321 161.1392  3.0337  4.0338
#> 51   55 2.7157 0.25374 161.1392  2.2146  3.2167
#> 52   56 2.7702 0.25639 161.1392  2.2639  3.2765
#> 53   57 3.9805 0.25173 161.1392  3.4834  4.4776
#> 54   58 2.9968 0.25319 161.1392  2.4968  3.4968
#> 55   59 3.3978 0.25250 161.1392  2.8991  3.8964
#> 56   60 3.0115 0.25303 161.1392  2.5118  3.5112
#> 57   61 3.2238 0.25378 161.1392  2.7227  3.7250
#> 58   62 3.2814 0.25266 161.1392  2.7825  3.7804
#> 59   63 3.5751 0.25500 161.1392  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(RCB, 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>
Example 8.1.1 from Experimental Design and Analysis for Tree Improvement

References

See also

Author

Examples
