-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathIntroToAgriculturalExperimentsinR.html
2192 lines (2132 loc) · 197 KB
/
IntroToAgriculturalExperimentsinR.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html >
<head>
<meta charset="UTF-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<title>Analysis of Common Agricultural Designs in R</title>
<meta name="description" content="A series of guidance, tutorials and examples for dealing with the data analysis of experiments following common agricultural designs using R and RStudio">
<meta name="generator" content="bookdown <!--bookdown:version--> and GitBook 2.6.7">
<meta property="og:title" content="Analysis of Common Agricultural Designs in R" />
<meta property="og:type" content="book" />
<meta property="og:description" content="A series of guidance, tutorials and examples for dealing with the data analysis of experiments following common agricultural designs using R and RStudio" />
<meta name="twitter:card" content="summary" />
<meta name="twitter:title" content="Analysis of Common Agricultural Designs in R" />
<meta name="twitter:description" content="A series of guidance, tutorials and examples for dealing with the data analysis of experiments following common agricultural designs using R and RStudio" />
<meta name="author" content="Sam Dumble">
<meta name="date" content="2018-11-15">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta name="apple-mobile-web-app-capable" content="yes">
<meta name="apple-mobile-web-app-status-bar-style" content="black">
<!--bookdown:link_prev-->
<!--bookdown:link_next-->
<script src="libs/jquery-2.2.3/jquery.min.js"></script>
<link href="libs/gitbook-2.6.7/css/style.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-bookdown.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-highlight.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-search.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-fontsettings.css" rel="stylesheet" />
<script src="libs/gitbook-2.6.7/js/app.min.js"></script>
<script src="libs/gitbook-2.6.7/js/lunr.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-search.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-sharing.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-fontsettings.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-bookdown.js"></script>
<script src="libs/gitbook-2.6.7/js/jquery.highlight.js"></script>
<script src="libs/htmlwidgets-1.3/htmlwidgets.js"></script>
<link href="libs/datatables-css-0.0.0/datatables-crosstalk.css" rel="stylesheet" />
<script src="libs/datatables-binding-0.5/datatables.js"></script>
<link href="libs/dt-core-1.10.16/css/jquery.dataTables.min.css" rel="stylesheet" />
<link href="libs/dt-core-1.10.16/css/jquery.dataTables.extra.css" rel="stylesheet" />
<script src="libs/dt-core-1.10.16/js/jquery.dataTables.min.js"></script>
<link href="libs/crosstalk-1.0.0/css/crosstalk.css" rel="stylesheet" />
<script src="libs/crosstalk-1.0.0/js/crosstalk.min.js"></script>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>
<link rel="stylesheet" href="style.css" type="text/css" />
</head>
<body>
<!--bookdown:title:start-->
<div id="header">
<h1 class="title">Analysis of Common Agricultural Designs in R</h1>
<h4 class="author"><em>Sam Dumble</em></h4>
<h4 class="date"><em>2018-11-15</em></h4>
</div>
<!--bookdown:title:end-->
<!--bookdown:toc:start-->
<div class="book without-animation with-summary font-size-2 font-family-1" data-basepath=".">
<div class="book-summary">
<nav role="navigation">
<!--bookdown:toc2:start-->
<ul>
<li><a href="#preface"><span class="toc-section-number">1</span> Preface</a></li>
<li><a href="#intro"><span class="toc-section-number">2</span> Introduction</a><ul>
<li><a href="#general-structure-r-syntax"><span class="toc-section-number">2.1</span> General Structure: R Syntax</a><ul>
<li><a href="#step-1-load-libraries"><span class="toc-section-number">2.1.1</span> Step 1: Load Libraries</a></li>
<li><a href="#step-2-import-data"><span class="toc-section-number">2.1.2</span> Step 2: Import Data</a></li>
<li><a href="#step-3-check-and-update-data"><span class="toc-section-number">2.1.3</span> Step 3: Check and update data</a></li>
<li><a href="#step-4.-explore-data"><span class="toc-section-number">2.1.4</span> Step 4. Explore data</a></li>
<li><a href="#step-5.-specify-a-model-for-data"><span class="toc-section-number">2.1.5</span> Step 5. Specify a model for data</a></li>
<li><a href="#step-6.-check-the-model"><span class="toc-section-number">2.1.6</span> Step 6. Check the model</a></li>
<li><a href="#step-7.-interpret-the-model"><span class="toc-section-number">2.1.7</span> Step 7. Interpret the model</a></li>
<li><a href="#step-8.-present-the-results-from-the-model"><span class="toc-section-number">2.1.8</span> Step 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#general-structure-explanation-of-each-step"><span class="toc-section-number">2.2</span> General Structure: Explanation of Each Step</a><ul>
<li><a href="#step-1-load-libraries-1"><span class="toc-section-number">2.2.1</span> Step 1: Load Libraries</a></li>
<li><a href="#step-2-import-data-1"><span class="toc-section-number">2.2.2</span> Step 2: Import Data</a></li>
<li><a href="#step-3-check-and-update-data-1"><span class="toc-section-number">2.2.3</span> Step 3: Check and update data</a></li>
<li><a href="#step-4.-explore-data-1"><span class="toc-section-number">2.2.4</span> Step 4. Explore data</a></li>
<li><a href="#step-5.-specify-a-model-for-data-1"><span class="toc-section-number">2.2.5</span> Step 5. Specify a model for data</a></li>
<li><a href="#step-6.-check-the-model-1"><span class="toc-section-number">2.2.6</span> Step 6. Check the model</a></li>
<li><a href="#step-7.-interpret-the-model-1"><span class="toc-section-number">2.2.7</span> Step 7. Interpret the model</a></li>
<li><a href="#step-8.-present-the-results-from-the-model-1"><span class="toc-section-number">2.2.8</span> Step 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#other-resources"><span class="toc-section-number">2.3</span> Other resources</a></li>
</ul></li>
<li><a href="#randomised-complete-block-design-rcbd"><span class="toc-section-number">3</span> Randomised Complete Block Design (RCBD)</a><ul>
<li><a href="#about-the-data"><span class="toc-section-number">3.1</span> About the data</a></li>
<li><a href="#section-1-steps-in-analysis-using-r"><span class="toc-section-number">3.2</span> Section 1: Steps in analysis using R</a></li>
<li><a href="#section-2-explanation-of-steps"><span class="toc-section-number">3.3</span> Section 2: Explanation of Steps</a><ul>
<li><a href="#install-r-packages-needed"><span class="toc-section-number">3.3.1</span> 1. Install R packages needed</a></li>
<li><a href="#import-data"><span class="toc-section-number">3.3.2</span> 2. Import data</a></li>
<li><a href="#check-and-update-data"><span class="toc-section-number">3.3.3</span> 3. Check and update data</a></li>
<li><a href="#explore-data"><span class="toc-section-number">3.3.4</span> 4. Explore data</a></li>
<li><a href="#specify-a-model-for-data"><span class="toc-section-number">3.3.5</span> 5. Specify a model for data</a></li>
<li><a href="#check-the-model"><span class="toc-section-number">3.3.6</span> 6. Check the model</a></li>
<li><a href="#interpret-model"><span class="toc-section-number">3.3.7</span> 7. Interpret Model</a></li>
<li><a href="#present-the-results-from-the-model"><span class="toc-section-number">3.3.8</span> 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#section-3-methodological-principles"><span class="toc-section-number">3.4</span> Section 3 – Methodological Principles</a></li>
</ul></li>
<li><a href="#split-plot-designs"><span class="toc-section-number">4</span> Split Plot Designs</a><ul>
<li><a href="#about-the-data-1"><span class="toc-section-number">4.1</span> About the data</a></li>
<li><a href="#section-1-steps-in-analysis-using-r-1"><span class="toc-section-number">4.2</span> Section 1: Steps in analysis using R</a></li>
<li><a href="#section-2-explanation-of-steps-1"><span class="toc-section-number">4.3</span> Section 2: Explanation of Steps</a><ul>
<li><a href="#install-r-packages-needed-1"><span class="toc-section-number">4.3.1</span> 1. Install R packages needed</a></li>
<li><a href="#import-data-1"><span class="toc-section-number">4.3.2</span> 2. Import data</a></li>
<li><a href="#check-and-update-data-1"><span class="toc-section-number">4.3.3</span> 3. Check and update data</a></li>
<li><a href="#explore-data-1"><span class="toc-section-number">4.3.4</span> 4. Explore data</a></li>
<li><a href="#plots-1"><span class="toc-section-number">4.3.5</span> Plots</a></li>
<li><a href="#summary-statistics-1"><span class="toc-section-number">4.3.6</span> Summary Statistics</a></li>
<li><a href="#specify-a-model-for-data-1"><span class="toc-section-number">4.3.7</span> 5. Specify a model for data</a></li>
<li><a href="#check-the-model-1"><span class="toc-section-number">4.3.8</span> 6. Check the model</a></li>
<li><a href="#interpret-model-1"><span class="toc-section-number">4.3.9</span> 7. Interpret Model</a></li>
<li><a href="#present-the-results-from-the-model-1"><span class="toc-section-number">4.3.10</span> 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#section-3-methodological-principles-1"><span class="toc-section-number">4.4</span> Section 3 – Methodological Principles</a></li>
</ul></li>
<li><a href="#adjusting-for-covariates"><span class="toc-section-number">5</span> Adjusting for Covariates</a><ul>
<li><a href="#about-the-data-2"><span class="toc-section-number">5.1</span> About the data</a></li>
<li><a href="#section-1-steps-in-analysis-using-r-2"><span class="toc-section-number">5.2</span> Section 1: Steps in analysis using R</a></li>
<li><a href="#section-2-explanation-of-steps-2"><span class="toc-section-number">5.3</span> Section 2: Explanation of Steps</a><ul>
<li><a href="#install-r-packages-needed-2"><span class="toc-section-number">5.3.1</span> 1. Install R packages needed</a></li>
<li><a href="#import-data-2"><span class="toc-section-number">5.3.2</span> 2. Import data</a></li>
<li><a href="#check-and-update-data-2"><span class="toc-section-number">5.3.3</span> 3. Check and update data</a></li>
<li><a href="#explore-data-2"><span class="toc-section-number">5.3.4</span> 4. Explore data</a></li>
<li><a href="#specify-a-model-for-data-2"><span class="toc-section-number">5.3.5</span> 5. Specify a model for data</a></li>
<li><a href="#check-the-model-2"><span class="toc-section-number">5.3.6</span> 6. Check the model</a></li>
<li><a href="#interpret-model-2"><span class="toc-section-number">5.3.7</span> 7. Interpret Model</a></li>
<li><a href="#present-the-results-from-the-model-2"><span class="toc-section-number">5.3.8</span> 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#section-3-methodological-principles-2"><span class="toc-section-number">5.4</span> Section 3 – Methodological Principles</a></li>
</ul></li>
<li><a href="#relay-planting-example"><span class="toc-section-number">6</span> Relay Planting Example</a><ul>
<li><a href="#section-1-steps-in-analysis-using-r-3"><span class="toc-section-number">6.1</span> Section 1: Steps in analysis using R</a></li>
<li><a href="#section-2-explanation-of-steps-3"><span class="toc-section-number">6.2</span> Section 2: Explanation of Steps</a><ul>
<li><a href="#install-r-packages-needed-3"><span class="toc-section-number">6.2.1</span> 1. Install R packages needed</a></li>
<li><a href="#import-data-3"><span class="toc-section-number">6.2.2</span> 2. Import data</a></li>
<li><a href="#check-and-update-data-3"><span class="toc-section-number">6.2.3</span> 3. Check and update data</a></li>
<li><a href="#explore-data-3"><span class="toc-section-number">6.2.4</span> 4. Explore data</a></li>
<li><a href="#specify-a-model-for-data-3"><span class="toc-section-number">6.2.5</span> 5. Specify a model for data</a></li>
<li><a href="#check-the-model-3"><span class="toc-section-number">6.2.6</span> 6. Check the model</a></li>
<li><a href="#interpret-model-3"><span class="toc-section-number">6.2.7</span> 7. Interpret Model</a></li>
<li><a href="#present-the-results-from-the-model-3"><span class="toc-section-number">6.2.8</span> 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#section-3-methodological-principles-3"><span class="toc-section-number">6.3</span> Section 3 – Methodological Principles</a></li>
</ul></li>
<li><a href="#multi-environment-trial-analysis"><span class="toc-section-number">7</span> Multi-Environment Trial Analysis</a><ul>
<li><a href="#section-1-steps-in-analysis-using-r-4"><span class="toc-section-number">7.1</span> Section 1: Steps in analysis using R</a></li>
<li><a href="#section-2-explanation-of-steps-4"><span class="toc-section-number">7.2</span> Section 2: Explanation of Steps</a><ul>
<li><a href="#install-r-packages-needed-4"><span class="toc-section-number">7.2.1</span> 1. Install R packages needed</a></li>
<li><a href="#import-data-4"><span class="toc-section-number">7.2.2</span> 2. Import data</a></li>
<li><a href="#check-and-update-data-4"><span class="toc-section-number">7.2.3</span> 3. Check and update data</a></li>
<li><a href="#explore-data-4"><span class="toc-section-number">7.2.4</span> 4. Explore data</a></li>
<li><a href="#specify-a-model-for-data-4"><span class="toc-section-number">7.2.5</span> 5. Specify a model for data</a></li>
<li><a href="#check-the-model-4"><span class="toc-section-number">7.2.6</span> 6. Check the model</a></li>
<li><a href="#interpret-model-4"><span class="toc-section-number">7.2.7</span> 7. Interpret Model</a></li>
<li><a href="#present-the-results-from-the-model-4"><span class="toc-section-number">7.2.8</span> 8. Present the results from the model</a></li>
</ul></li>
<li><a href="#section-3-methodological-principles-4"><span class="toc-section-number">7.3</span> Section 3 – Methodological Principles</a></li>
</ul></li>
<li><a href="#references">References</a></li>
</ul>
<!--bookdown:toc2:end-->
</nav>
</div>
<div class="book-body">
<div class="body-inner">
<div class="book-header" role="navigation">
<h1>
<i class="fa fa-circle-o-notch fa-spin"></i><a href="./">Analysis of Common Agricultural Designs in R</a>
</h1>
</div>
<div class="page-wrapper" tabindex="-1" role="main">
<div class="page-inner">
<section class="normal" id="section-">
<!--bookdown:toc:end-->
<!--bookdown:body:start-->
<div id="preface" class="section level1">
<h1><span class="header-section-number">1</span> Preface</h1>
<p>All of these tutorials assume that you have already been able to install R and RStudio onto your computer and that you have a reliable internet connection. For help with orientation of R for new users please see [add cross reference to an intro document].<br />
1. RCBDs (Randomised complete block design) [add cross reference]<br />
2. Split Plot Design [add cross reference]<br />
3. Adjusting for Covariates [add cross reference]<br />
4. Factorial designs and interactions [add cross reference]<br />
5. Multi Environment Trials [add cross reference]</p>
<!--chapter:end:index.Rmd-->
</div>
<div id="intro" class="section level1">
<h1><span class="header-section-number">2</span> Introduction</h1>
<p>Different designs require different models. But in R, nearly all other steps are identical before and after model fitting – assuming that, regardless of the design, you are interested in more or less the same question:</p>
<p>Assessing how a numeric response variable (e.g. yield) varies by a treatment factor, or factors.</p>
<p>Being able to learn and understand these steps, will let you analyse any data you have available from on-station trials! In these guides we will use the lmer function within R to fit (nearly) all the models we may want to consider for these agricultural designs. This fits a linear mixed effects regression model. A detailed explanation of these statistical models, and their applicability to agricultural analyses can be found here: <a href="https://www.jic.ac.uk/services/statistics/readingadvice/booklets/topmix.html" class="uri">https://www.jic.ac.uk/services/statistics/readingadvice/booklets/topmix.html</a> . In short, these models enable us to separate out factors that are of interest to us (e.g. treatments, varieties) to factors which are not of interest to us, but that still introduce (e.g. blocks). On-farm trials, less standard designs, and more complex outcome variables (e.g. disease scores, incidence rates, growth patterns) may require more care with analysis and more consideration in how to analyse and interpret results. Many of the general principles are the same, as is a large portion of the R syntax, but in these cases more care is needed to ensure a coherent analysis. There is no “recipe” which will work in the same way every time, each analysis may bring up new or unexpected considerations that need to be addressed rather than forcing the analysis to fit within a standard framework.</p>
<div id="general-structure-r-syntax" class="section level2">
<h2><span class="header-section-number">2.1</span> General Structure: R Syntax</h2>
<div id="step-1-load-libraries" class="section level3">
<h3><span class="header-section-number">2.1.1</span> Step 1: Load Libraries</h3>
<pre><code>library(ggplot2)
library(emmeans)
library(doBy)
library(lmerTest)
library(multcompView)</code></pre>
</div>
<div id="step-2-import-data" class="section level3">
<h3><span class="header-section-number">2.1.2</span> Step 2: Import Data</h3>
<pre><code>mydata <- read.csv("C:/Users/Admin/Desktop/mydata.csv")</code></pre>
</div>
<div id="step-3-check-and-update-data" class="section level3">
<h3><span class="header-section-number">2.1.3</span> Step 3: Check and update data</h3>
<pre><code>summary(mydata)
str(mydata)
mydata$treatment<-factor(mydata$treatment)</code></pre>
</div>
<div id="step-4.-explore-data" class="section level3">
<h3><span class="header-section-number">2.1.4</span> Step 4. Explore data</h3>
<pre><code>ggplot(data= mydata,aes(y=response,x=treatment,col= block))+
geom_point()
summaryBy(response ~ treatment, data= mydata, FUN=c(mean,median,sd))</code></pre>
</div>
<div id="step-5.-specify-a-model-for-data" class="section level3">
<h3><span class="header-section-number">2.1.5</span> Step 5. Specify a model for data</h3>
<pre><code>mymodel <-lmer(response~treatment+(1|block), data=mydata)</code></pre>
</div>
<div id="step-6.-check-the-model" class="section level3">
<h3><span class="header-section-number">2.1.6</span> Step 6. Check the model</h3>
<pre><code>plot(mymodel)
qqnorm(resid(mymodel))
qqline(resid(mymodel))</code></pre>
</div>
<div id="step-7.-interpret-the-model" class="section level3">
<h3><span class="header-section-number">2.1.7</span> Step 7. Interpret the model</h3>
<pre><code>anova(mymodel, ddf="Kenward-Roger")
print(VarCorr(mymodel), comp=("Variance"))</code></pre>
</div>
<div id="step-8.-present-the-results-from-the-model" class="section level3">
<h3><span class="header-section-number">2.1.8</span> Step 8. Present the results from the model</h3>
<pre><code>emmip(mymodel,~treatment,CIs = TRUE)
emmeans(mymodel, ~ treatment)
cld(emmeans(mymodel, ~ treatment))</code></pre>
</div>
</div>
<div id="general-structure-explanation-of-each-step" class="section level2">
<h2><span class="header-section-number">2.2</span> General Structure: Explanation of Each Step</h2>
<div id="step-1-load-libraries-1" class="section level3">
<h3><span class="header-section-number">2.2.1</span> Step 1: Load Libraries</h3>
<p>R is an open-source piece of software. One major benefit of this is that many useful functions for importing, manipulating, analysing and presenting data have been created by other R users, beyond what is available in the “base” R packages. Many of these functions are implemented in packages, or libraries, which need to be downloaded and installed separately from your main R and RStudio installation. The main ones you will need to be able to follow this set of guides are: <code>ggplot2</code>: A powerful graphing package, allowing high quality graphs to be produced <code>doBy</code>: A package for easy and customisable calculations of summary statistics <code>lmerTest</code>: A package for fitting and evaluating linear mixed effects regression models, using REML (restricted maximum likelihood) methods as found in Genstat <code>emmeans</code>: A package to calculate estimated marginal means and confidence intervals from statistical models. Similar to EMMEANS in Genstat or LSMEANS in SAS. <code>multcompView</code>: A package for conducting mean separation analysis from mixed effects regression models To install these packages onto your computer you need an internet connection, and for a clean installation of R and RStudio onto your computer. You only need to install an R package once, using install.packages() or through the menus, but you do need to load the packages every time you come to use them using library(). You can learn more about libraries here: <a href="https://www.datacamp.com/community/tutorials/r-packages-guide" class="uri">https://www.datacamp.com/community/tutorials/r-packages-guide</a></p>
</div>
<div id="step-2-import-data-1" class="section level3">
<h3><span class="header-section-number">2.2.2</span> Step 2: Import Data</h3>
<pre><code>mydata <- read.csv("C:/Users/Admin/Desktop/mydata.csv")</code></pre>
<p>Key things to consider before even attempting to read your data into R: • Is your data in a single sheet, in a continuous rectangle, with no blank rows or columns? • Is there a single row at the top of your data containing the variable names? • Are the variable names concise, but informative, and contain no spaces or punctuation? • Are missing values consistently coded in your dataset? • Are factor levels consistently coded in your dataset (even including case sensitive – R will consider “treatment A” and “Treatment a” as 2 different treatments. • If you have dates in your data then are they always written in the same format? You can learn more about some of the important considerations of preparing your data for importing into R here: <a href="http://www.sthda.com/english/wiki/best-practices-in-preparing-data-files-for-importing-into-r" class="uri">http://www.sthda.com/english/wiki/best-practices-in-preparing-data-files-for-importing-into-r</a></p>
</div>
<div id="step-3-check-and-update-data-1" class="section level3">
<h3><span class="header-section-number">2.2.3</span> Step 3: Check and update data</h3>
<p>There can be many unforeseen issues when importing your dataset if it is not cleaned in the way you would like it to be. Checking the data, both visually, and using functions like summary() and str() can help you see if there have been any issues which may need addressing. Common problems you might see at this point would be: • Variable names changing: if your variable names contained spaces, or punctuation, then R will change them and introduce extra dots into the name. Ideally you want variable names in R to be concise, and contain no punctuation. This will make writing the syntax much easier • Missing value codes: If you have missing values in your dataset, check that R has imported these as missing values. If in Excel you have a blank cell then this will be imported correctly into R. If you are using a code (like -999 for example) R will not automatically recognise this as a missing value. • Factors being treated as numbers</p>
<p>These are largely the same concerns as in step 2; but being checked from within R rather than within Excel.</p>
<p>Why is it important to make sure factor variables are treated as factors?</p>
<p>We are often taught to use codes when entering and collecting data for categorical variables, such as treatment or variety. If we use numeric codes, i.e. 1,2,3,4 for 4 treatments, then we can potentially see problems with our analysis unless we specify explicitly that this is the case. This problem is not an issue is we use non-numeric codes for treatments, e.g. A,B,C,D. The same data is presented below twice; once with estimates of the treatment means from an analysis of a numeric treatment variable and once from a factor treatment variable.</p>
<p>With the numeric variable the model tries to fit the treatment effect as if it is a continuous scale; i.e. that treatment 2 is 1 point higher than treatment 1. With the factor variable the model treats all 4 treatment groups as being independent of each other. In this case if we had not converted the treatment to a factor we would have had a completely useless model, telling us that there was no treatment effect and providing severe over-estimates of treatments 2 and 4 and a sever underestimate of treatment 3. In fact there is a very highly significant treatment effect in this data, which can only be identified from the analysis when the variable is treated as a factor.</p>
</div>
<div id="step-4.-explore-data-1" class="section level3">
<h3><span class="header-section-number">2.2.4</span> Step 4. Explore data</h3>
<p>Exploratory analysis helps us to understand the results we have found in our data. It can show us • if there are clear effects from visual inspection • the magnitude of any effects, • the variability in our results • if our data is distributed in a way that will lead to a standard modelling approach We can also calculate summary statistics, such as means and percentages.</p>
<p><a href="http://r4ds.had.co.nz/exploratory-data-analysis.html" class="uri">http://r4ds.had.co.nz/exploratory-data-analysis.html</a></p>
</div>
<div id="step-5.-specify-a-model-for-data-1" class="section level3">
<h3><span class="header-section-number">2.2.5</span> Step 5. Specify a model for data</h3>
<pre><code>mymodel <-lmer(response~treatment+(1|block),data=mydata)</code></pre>
<p>Cross link to slides of examples for model construction.</p>
</div>
<div id="step-6.-check-the-model-1" class="section level3">
<h3><span class="header-section-number">2.2.6</span> Step 6. Check the model</h3>
<pre><code>plot(mymodel)</code></pre>
<p>There are three main assumptions that are worthwhile considering when assessing if the model being fitted is valid from a statistical perspective. 1. “Independence”: This assumption can be met by including the dependencies within the design of the experiment within the model through the use of random effects. For example - two plots within the same block, may have some level of inter-relatedness. Including a “block” term in the model allows this assumption to be met in this instance. 2. “Homogeneity”: This assumption relates to whether the variability in each treatment group is similar. In order to calculate standard errors and p-values from the model an assumption is made that there is constant variance across all treatments. If this assumption does not hold then these standard errors and p-values will not be accurate. It is common in many situations to have more variability in high yielding treatments than in low yielding treatments. E.g. 4 treatments, each replicated 8 times</p>
</div>
<div id="step-7.-interpret-the-model-1" class="section level3">
<h3><span class="header-section-number">2.2.7</span> Step 7. Interpret the model</h3>
<pre><code>anova(mymodel, ddf="Kenward-Roger")
print(VarCorr(mymodel), comp=("Variance"))</code></pre>
<p>Summary of Kenward-Rogers degree of freedom from mixed models: <a href="https://www.jstatsoft.org/article/view/v082i13/v82i13.pdf" class="uri">https://www.jstatsoft.org/article/view/v082i13/v82i13.pdf</a></p>
</div>
<div id="step-8.-present-the-results-from-the-model-1" class="section level3">
<h3><span class="header-section-number">2.2.8</span> Step 8. Present the results from the model</h3>
<pre><code>emmip(mymodel,~treatment,CIs = TRUE)
emmeans(mymodel, ~ treatment)
cld(emmeans(mymodel, ~ treatment))
</code></pre>
<p><a href="https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html" class="uri">https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html</a> <a href="https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html" class="uri">https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html</a></p>
</div>
</div>
<div id="other-resources" class="section level2">
<h2><span class="header-section-number">2.3</span> Other resources</h2>
<p>For other agricultural trials, particularly if you have slightly different hypotheses to this standard framework, this provides a useful resource and overview of using R for agricultural analyses: <a href="http://rstats4ag.org" class="uri">http://rstats4ag.org</a> There are also specific examples of agricultural experiments with more complex designs, particularly in dealing with repeated measurements over time, in an R package called agriTutorial, which provides 5 specific case-studies of analysing field trial data in R. <a href="https://cran.r-project.org/web/packages/agriTutorial/agriTutorial.pdf" class="uri">https://cran.r-project.org/web/packages/agriTutorial/agriTutorial.pdf</a></p>
<!--chapter:end:01-intro.Rmd-->
</div>
</div>
<div id="randomised-complete-block-design-rcbd" class="section level1">
<h1><span class="header-section-number">3</span> Randomised Complete Block Design (RCBD)</h1>
<p>Aim: make it easy to do standard analysis of standard experimental designs used in field trials Assumptions: you know some basic R, have R and RStudio already installed on your compuiter and you are familiar with the standard analyses of field trials.<br />
This document will focus initially on the simple analysis of an RCBD trial using R. Section 1 provides the steps used to produce the analysis; Section 2 provides some commentary on how these commands work, what output is created, and why these commands were chosen; Section 3 deals with aspects of the statistical methodology.</p>
<div id="about-the-data" class="section level2">
<h2><span class="header-section-number">3.1</span> About the data</h2>
<p>The data used in this example is from a study was conducted in Eastern Zambia and the main aim was to improve on the efficiency of the natural fallows by using appropriate trees that may have relevance in soil fertility regeneration within permissible fallow periods.</p>
<p>The design was a randomized complete block design experiment with 4 blocks and 9 treatments was conducted. The primary outcome variable was crop yield (yield).</p>
<p>The objective for this analysis is to study the impact of different fallow types on crop yields.</p>
<p>The following steps were followed to generate the output in this document. The data was organized in excel rectangle columns with the different variables appearing in excel columns. All data checks were done in excel, meaningful data was selected and a copy of this data file was stored as a CSV file to make data import easy in R. The data file used in this analysis can be downloaded here: <a href="https://bit.ly/2rfLBEt" class="uri">https://bit.ly/2rfLBEt</a></p>
</div>
<div id="section-1-steps-in-analysis-using-r" class="section level2">
<h2><span class="header-section-number">3.2</span> Section 1: Steps in analysis using R</h2>
<ol style="list-style-type: decimal">
<li>Install R packages needed</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(ggplot2)
<span class="kw">library</span>(emmeans)
<span class="kw">library</span>(doBy)
<span class="kw">library</span>(lmerTest)
<span class="kw">library</span>(multcompView)</code></pre></div>
<ol start="2" style="list-style-type: decimal">
<li>Import data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/Fallow N2.csv"</span>)</code></pre></div>
<ol start="3" style="list-style-type: decimal">
<li>Check and update data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fallow)
<span class="kw">str</span>(fallow)
fallow<span class="op">$</span>rep<-<span class="kw">factor</span>(fallow<span class="op">$</span>rep)
fallow<span class="op">$</span>plot<-<span class="kw">factor</span>(fallow<span class="op">$</span>plot)</code></pre></div>
<ol start="4" style="list-style-type: decimal">
<li>Explore data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>treat,<span class="dt">col=</span>rep))<span class="op">+</span><span class="kw">geom_point</span>()
<span class="kw">summaryBy</span>(yield<span class="op">~</span>treat, <span class="dt">data=</span>fallow, <span class="dt">FUN=</span><span class="kw">c</span>(min,max,mean,median,sd))</code></pre></div>
<ol start="5" style="list-style-type: decimal">
<li>Specify a model for data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">rcbdmodel1<-<span class="kw">lmer</span>(yield<span class="op">~</span>treat<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>rep),<span class="dt">data=</span>fallow)</code></pre></div>
<ol start="6" style="list-style-type: decimal">
<li>Check the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(rcbdmodel1)
<span class="kw">qqnorm</span>(<span class="kw">resid</span>(rcbdmodel1))
<span class="kw">qqline</span>(<span class="kw">resid</span>(rcbdmodel1))</code></pre></div>
<ol start="7" style="list-style-type: decimal">
<li>Interpret the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">anova</span>(rcbdmodel1,<span class="dt">ddf=</span><span class="st">"Kenward-Roger"</span>)
<span class="kw">print</span>(<span class="kw">VarCorr</span>(rcbdmodel1), <span class="dt">comp=</span>(<span class="st">"Variance"</span>))</code></pre></div>
<ol start="8" style="list-style-type: decimal">
<li>Present the results from the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmip</span>(rcbdmodel1,<span class="op">~</span>treat,<span class="dt">CIs =</span> <span class="ot">TRUE</span>)
<span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat)
<span class="kw">cld</span>(<span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat))</code></pre></div>
</div>
<div id="section-2-explanation-of-steps" class="section level2">
<h2><span class="header-section-number">3.3</span> Section 2: Explanation of Steps</h2>
<div id="install-r-packages-needed" class="section level3">
<h3><span class="header-section-number">3.3.1</span> 1. Install R packages needed</h3>
<p>A number of packages following packages were used during data exploration and analysis. For a general introduction explaining what R packages are and how they work, this is a really useful guide <a href="https://www.datacamp.com/community/tutorials/r-packages-guide" class="uri">https://www.datacamp.com/community/tutorials/r-packages-guide</a>. For each of these packages to be installed, using install.packages(), this requires a reliable internet connection and a correctly installed version of R and RStudio. If you are having difficulties installing these packages please ask for help.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"ggplot2"</span>)
<span class="kw">library</span>(ggplot2)</code></pre></div>
<p><code>ggplot2</code> This package provides a powerful graphics language for creating elegant and complex graphs in R.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"emmeans"</span>)
<span class="kw">library</span>(emmeans)</code></pre></div>
<p><code>emmeans</code> Estimated marginal means (also known as least squares means) helps provide expected mean values and confidence intervals from statistical models.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"doBy"</span>)
<span class="kw">library</span>(doBy)</code></pre></div>
<p><code>doBy</code>Allows easy production of summary statistic tables</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"lmerTest"</span>)
<span class="kw">library</span>(lmerTest)</code></pre></div>
<p><code>lmerTest</code> Allows produce of flexible mixed effects regression models, similar to REML in Genstat.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"multcompView"</span>)
<span class="kw">library</span>(multcompView)</code></pre></div>
<p><code>multcompView</code> allows for mean seperation methods on analyses</p>
</div>
<div id="import-data" class="section level3">
<h3><span class="header-section-number">3.3.2</span> 2. Import data</h3>
<p>Our data set saved as a CSV file, so we can use the read.csv commmand to import the data. We are going to assign the name of the data with R to be <code>fallow2</code>. Remember in R Studio you could also use the “Import Dataset” menu to import a dataset.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/Fallow N2.csv"</span>)</code></pre></div>
</div>
<div id="check-and-update-data" class="section level3">
<h3><span class="header-section-number">3.3.3</span> 3. Check and update data</h3>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">DT<span class="op">::</span><span class="kw">datatable</span>(fallow)</code></pre></div>
<div id="htmlwidget-845525c2c62065a3b2aa" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-845525c2c62065a3b2aa">{"x":{"filter":"none","data":[["1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36"],[1,4,4,1,1,3,3,1,3,2,2,4,2,1,3,2,3,4,3,3,3,2,1,4,4,2,2,1,2,4,4,1,2,3,1,4],[2,3,6,9,7,3,8,6,9,9,5,5,3,8,5,6,4,8,6,1,2,7,4,7,1,8,2,3,4,9,4,1,1,7,5,2],["8 maize","5 C.siamea","8 maize","7 groundnut","5 C.siamea","8 maize","5 C.siamea","9 nat fallow","6 C.calo","5 C.siamea","6 C.calo","6 C.calo","8 maize","6 C.calo","9 nat fallow","9 nat fallow","7 groundnut","9 nat fallow","4 F.congesta","2 G.sepium","3 L.leuco","4 F.congesta","3 L.leuco","4 F.congesta","7 groundnut","3 L.leuco","7 groundnut","4 F.congesta","2 G.sepium","2 G.sepium","3 L.leuco","2 G.sepium","1 S.sesban","1 S.sesban","1 S.sesban","1 S.sesban"],[1.14,1.74,1.95,2.06,2.09,2.15,2.21,2.22,2.34,2.38,2.62,2.65,2.66,2.74,2.9,3,3.03,3.12,3.16,3.19,3.39,3.45,3.48,3.48,3.52,3.62,3.72,3.75,3.85,3.95,4.21,4.39,4.72,5.29,5.64,6.54],[2798,0,1787,129,1,1144,0,228,0,0,0,105,400,71,0,203,120,89,1,1371,0,49,0,1741,270,0,1,1,721,755,41,1,0,1,0,0]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th> <\/th>\n <th>rep<\/th>\n <th>plot<\/th>\n <th>treat<\/th>\n <th>yield<\/th>\n <th>striga<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-right","targets":[1,2,4,5]},{"orderable":false,"targets":0}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<p>When reading data into R it is always useful to check that data is in the format expected. How many variables are there? How many rows? How have the columns been read in? The summary command can help to show if the data is being treated correctly.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fallow)</code></pre></div>
<pre><code>## rep plot treat yield
## Min. :1.00 Min. :1 1 S.sesban : 4 Min. :1.140
## 1st Qu.:1.75 1st Qu.:3 2 G.sepium : 4 1st Qu.:2.370
## Median :2.50 Median :5 3 L.leuco : 4 Median :3.140
## Mean :2.50 Mean :5 4 F.congesta: 4 Mean :3.232
## 3rd Qu.:3.25 3rd Qu.:7 5 C.siamea : 4 3rd Qu.:3.728
## Max. :4.00 Max. :9 6 C.calo : 4 Max. :6.540
## (Other) :12
## striga
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 21.0
## Mean : 334.1
## 3rd Qu.: 238.5
## Max. :2798.0
## </code></pre>
<p>Where data is being treated as a numeric variable (i.e. a number) <code>summary</code> provides statistics like the mean, min and max. Where data is being treated like a categorical variable (i.e. a group) then summary provides frequency tables.</p>
<p>From the results we can see that the variables rep and plot are being considered as numeric variables. However these are grouping variables, not number variables, the numbers used are simply codes. If we do not rectify this then our analysis later will be incorrect and meaningless.<br />
This can also be seen more explicitly using the str() function.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">str</span>(fallow)</code></pre></div>
<pre><code>## 'data.frame': 36 obs. of 5 variables:
## $ rep : int 1 4 4 1 1 3 3 1 3 2 ...
## $ plot : int 2 3 6 9 7 3 8 6 9 9 ...
## $ treat : Factor w/ 9 levels "1 S.sesban","2 G.sepium",..: 8 5 8 7 5 8 5 9 6 5 ...
## $ yield : num 1.14 1.74 1.95 2.06 2.09 2.15 2.21 2.22 2.34 2.38 ...
## $ striga: int 2798 0 1787 129 1 1144 0 228 0 0 ...</code></pre>
<p>So we need to convert these variables into factors.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow<span class="op">$</span>rep<-<span class="kw">factor</span>(fallow<span class="op">$</span>rep)
fallow<span class="op">$</span>plot<-<span class="kw">factor</span>(fallow<span class="op">$</span>plot)</code></pre></div>
<p>These commands take the column rep within the data frame fallow, converts into a factor and saves the result in a column called rep within fallow.</p>
</div>
<div id="explore-data" class="section level3">
<h3><span class="header-section-number">3.3.4</span> 4. Explore data</h3>
<div id="plots" class="section level4">
<h4><span class="header-section-number">3.3.4.1</span> Plots</h4>
<p>With this code we want to summarize data fallow by yield as the response and treatment as a factor using points.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>treat))<span class="op">+</span><span class="kw">geom_point</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-22-1.png" width="672" /> We could also extend this to identify which points came from which reps.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>treat,<span class="dt">col=</span>rep))<span class="op">+</span><span class="kw">geom_point</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-23-1.png" width="672" /> Using ggplot2 we can easily change between different types of graph with small changes to the code. Boxplots are very useful if we have lots of data in each group, but in this example we only have 4 points so it is easy to visualise all of our data using a scatter plot. But the only change we would need to make to our original code is to change geom_point() to geom_boxplot().</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>treat))<span class="op">+</span><span class="kw">geom_boxplot</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-24-1.png" width="672" /></p>
<p>From the figures produced we can see that treatment 1 has consistently high yields. The lowest yield recorded for treatment 1 is higher than the highest yield recorded for any of the other treatments. Treatments 5 and 8 had consistently low yields.</p>
</div>
<div id="summary-statistics" class="section level4">
<h4><span class="header-section-number">3.3.4.2</span> Summary Statistics</h4>
<p>To produce summary statistics, by group, there are many options within R. One option is to use the summaryBy function, from the doBy library. The code used for this is quite similar to the code we will use to produce models in a later step.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summaryBy</span>(yield<span class="op">~</span>treat, <span class="dt">data=</span>fallow, <span class="dt">FUN=</span>mean)</code></pre></div>
<pre><code>## treat yield.mean
## 1 1 S.sesban 5.5475
## 2 2 G.sepium 3.8450
## 3 3 L.leuco 3.6750
## 4 4 F.congesta 3.4600
## 5 5 C.siamea 2.1050
## 6 6 C.calo 2.5875
## 7 7 groundnut 3.0825
## 8 8 maize 1.9750
## 9 9 nat fallow 2.8100</code></pre>
<p>We can also calculate multiple statistics in the same line of code</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summaryBy</span>(yield<span class="op">~</span>treat, <span class="dt">data=</span>fallow, <span class="dt">FUN=</span><span class="kw">c</span>(min,max,mean,median,sd))</code></pre></div>
<pre><code>## treat yield.min yield.max yield.mean yield.median yield.sd
## 1 1 S.sesban 4.72 6.54 5.5475 5.465 0.7625997
## 2 2 G.sepium 3.19 4.39 3.8450 3.900 0.4956813
## 3 3 L.leuco 3.39 4.21 3.6750 3.550 0.3690077
## 4 4 F.congesta 3.16 3.75 3.4600 3.465 0.2412468
## 5 5 C.siamea 1.74 2.38 2.1050 2.150 0.2708628
## 6 6 C.calo 2.34 2.74 2.5875 2.635 0.1726992
## 7 7 groundnut 2.06 3.72 3.0825 3.275 0.7407372
## 8 8 maize 1.14 2.66 1.9750 2.050 0.6318491
## 9 9 nat fallow 2.22 3.12 2.8100 2.950 0.4034848</code></pre>
</div>
</div>
<div id="specify-a-model-for-data" class="section level3">
<h3><span class="header-section-number">3.3.5</span> 5. Specify a model for data</h3>
<p>In this design, an RCBD, we have one treatment factor, “treat”, and one layout factor “rep”. More information about model fitting can be found in section 2.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">rcbdmodel1<-<span class="kw">lmer</span>(yield<span class="op">~</span>treat<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>rep),<span class="dt">data=</span>fallow)</code></pre></div>
<p>R is unlike many other software packages in how it fits models. The best way of handling models in R is to assign the model to a name (in this case rcbdmodel1) and then ask R to provide different sorts of output for this model. When you run the above line you will get now output from the data - this is what we expected to see!</p>
</div>
<div id="check-the-model" class="section level3">
<h3><span class="header-section-number">3.3.6</span> 6. Check the model</h3>
<p>Before interpretting the model any further we should investigate the model validity, to ensure any conclusions we draw are valid. There are 3 assumptions that we can check for using standard model checking plots. 1. Homogeneity (equal variance) 2. Values with high leverage 3. Normality of residuals</p>
<p>The function plot() when used with a model will plot the fitted values from the model against the expected values.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(rcbdmodel1)</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-28-1.png" width="672" /> The residual Vs fitted plot is a scatter plot of the Residuals on the y-axis and the fitted on the x-axis and the aim for this plot is to test the assumption of equal variance of the residuals across the range of fitted values. Since the residuals do not funnel out (to form triangular/diamond shape) the assumption of equal variance is met.</p>
<p>We can also see that there are no extreme values in the residuals which might be potentially causing problems with the validity of our conclusions (leverage)</p>
<p>To assess the assumption of normality we can produce a qqplot. This shows us how closely the residuals follow a normal distribution - if there are severe and syste,matic deviations from the line then we may want to consider an alternative distribution.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">qqnorm</span>(<span class="kw">resid</span>(rcbdmodel1))
<span class="kw">qqline</span>(<span class="kw">resid</span>(rcbdmodel1))</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-29-1.png" width="672" /> In this case the residuals seem to fit the assumption required for normality.</p>
</div>
<div id="interpret-model" class="section level3">
<h3><span class="header-section-number">3.3.7</span> 7. Interpret Model</h3>
<p>The anova() function only prints the rows of analysis of variance table for treatment effects when looking at a mixed model fitted using lmer().</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">anova</span>(rcbdmodel1,<span class="dt">ddf=</span><span class="st">"Kenward-Roger"</span>)</code></pre></div>
<pre><code>## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treat 37.806 4.7258 8 24 20.146 6.981e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>ddf=Kenward-Roger tells R which method to use for determining the calculations of the table; this option matches the defaults found within SAS or Genstat. The ANOVA table suggests a highly significant effect of the treatment on the yield.</p>
<p>To obtain the residual variance, and the variance attributed to the blocks we need an additional command. From these number it is possible to reconstruct a more classic ANOVA table, if so desired.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">print</span>(<span class="kw">VarCorr</span>(rcbdmodel1), <span class="dt">comp=</span>(<span class="st">"Variance"</span>))</code></pre></div>
<pre><code>## Groups Name Variance
## rep (Intercept) 0.013817
## Residual 0.234577</code></pre>
</div>
<div id="present-the-results-from-the-model" class="section level3">
<h3><span class="header-section-number">3.3.8</span> 8. Present the results from the model</h3>
<p>To help understand what the significant result from the ANOVA table means we can produce several plots and tables to help us. First we can use the function emmip() to produce plots of the modelled results, including 95% confidence intervals.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmip</span>(rcbdmodel1,<span class="op">~</span>treat,<span class="dt">CIs =</span> <span class="ot">TRUE</span>)</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-32-1.png" width="672" /></p>
<p>To obtain the numbers used in creating this graph we can use the function emmeans.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat)</code></pre></div>
<pre><code>## treat emmean SE df lower.CL upper.CL
## 1 S.sesban 5.5475 0.2491955 26.35 5.0356 6.0594
## 2 G.sepium 3.8450 0.2491955 26.35 3.3331 4.3569
## 3 L.leuco 3.6750 0.2491955 26.35 3.1631 4.1869
## 4 F.congesta 3.4600 0.2491955 26.35 2.9481 3.9719
## 5 C.siamea 2.1050 0.2491955 26.35 1.5931 2.6169
## 6 C.calo 2.5875 0.2491955 26.35 2.0756 3.0994
## 7 groundnut 3.0825 0.2491955 26.35 2.5706 3.5944
## 8 maize 1.9750 0.2491955 26.35 1.4631 2.4869
## 9 nat fallow 2.8100 0.2491955 26.35 2.2981 3.3219
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95</code></pre>
<p>And one method for conducting mean separation analysis we can use the function cld().</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">cld</span>(<span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat))</code></pre></div>
<pre><code>## treat emmean SE df lower.CL upper.CL .group
## 8 maize 1.9750 0.2491955 26.35 1.4631 2.4869 1
## 5 C.siamea 2.1050 0.2491955 26.35 1.5931 2.6169 1
## 6 C.calo 2.5875 0.2491955 26.35 2.0756 3.0994 12
## 9 nat fallow 2.8100 0.2491955 26.35 2.2981 3.3219 123
## 7 groundnut 3.0825 0.2491955 26.35 2.5706 3.5944 123
## 4 F.congesta 3.4600 0.2491955 26.35 2.9481 3.9719 23
## 3 L.leuco 3.6750 0.2491955 26.35 3.1631 4.1869 23
## 2 G.sepium 3.8450 0.2491955 26.35 3.3331 4.3569 3
## 1 S.sesban 5.5475 0.2491955 26.35 5.0356 6.0594 4
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 9 estimates
## significance level used: alpha = 0.05</code></pre>
<p>In the output, groups sharing a letter in the .group are not statistically different from each other.</p>
</div>
</div>
<div id="section-3-methodological-principles" class="section level2">
<h2><span class="header-section-number">3.4</span> Section 3 – Methodological Principles</h2>
<p>There are always many different ways of doing all that we have done here in R. The less complex the method/code is, the better it is for you so that you can easily grasp the method.</p>
<p>For instance, we have fitted our model as a linear mixed effect model rather than traditional ANOVA because lmer model has the following advantages:</p>
<ol style="list-style-type: decimal">
<li>They are very flexible especially where we have repeated measures, for instance you don’t need to have the same number of observations per subject/treatment.</li>
<li>Ability to account for a series of random effects. Not only are farms/farmers/plots…. different from each other, but things with in farms/plots….. also differ . Not taking these sources of variation into account will lead to underestimations of accuracy.</li>
<li>Allows for generalization of non-normal data.<br />
</li>
<li>Handling missing data: If the percentage of missing data is small and that data missing is a random sample of the data set,data from the observations with missing data can be analysed with lmer (unlike other packages that would do listwise deletion.</li>
<li>Takes into account variation that is explained by the predictor variables of interest ie fixed effects and variation that is not explained by these predictors ie random effects.</li>
</ol>
<p>Not forgetting that selecting variables to include in our model generally depends on theory, statistics and practical knowledge the following (general) rules will be considered while fitting our models:</p>
<ol style="list-style-type: lower-roman">
<li>Consider the Treatments (A, B,….) as fixed effects and hence presented as A*B in our model.</li>
<li>Consider the layout factors as random effects and hence presented as (1|block/plot…) in our model. Generally, our model is in the form of Model<-lmer(Response~ (1|Block/Plot)+Treatment A + Treatment B…, data=Dataframe)</li>
</ol>
<p>In this example using the fallow data, note that if we had a “completely randomised” design rather than a “blocked randomised design”, where each treatment was replicated 4 times but there were not blocks, this is a rare example of a design which cannot be handled by lmer. In this case there would be no random effects, so the function needed would be lm() rather than lmer().</p>
<p>Food for thought: Your best model will certainly be as good as the data you collected!!!</p>
<!--chapter:end:02-RCBD.Rmd-->
</div>
</div>
<div id="split-plot-designs" class="section level1">
<h1><span class="header-section-number">4</span> Split Plot Designs</h1>
<p>Aim: make it easy to do standard analysis of standard experimental designs used in field trials Assumptions: you know some basic R, have R and RStudio already installed on your compuiter and you are familiar with the standard analyses of field trials.<br />
This document will focus initially on the simple analysis of a split plot design trial using R. Section 1 provides the steps used to produce the analysis; Section 2 provides some commentary on how these commands work, what output is created, and why these commadns were chosen; Section 3 deals with aspects of the statistical methodology.</p>
<p>It would be beneficial to also read through Part 1 in this series, analysis of RCBD single factor experiments. You may notice many similarities in the R syntax used in these guides.</p>
<div id="about-the-data-1" class="section level2">
<h2><span class="header-section-number">4.1</span> About the data</h2>
<p>The data for this example involves a split plot designed experiment. Treatments are 4 cropping patterns, and two nitrogen levels. The design is a split Both N and P could limit maize growth in the –N subplots, whereas N will not limit maize growth in the +N subplots. The comparison of +N and –N subplots within a mainplot will assess whether the fallows have eliminated N deficiency for maize.</p>
<p>Differences in maize yield among treatments for the +N subplot will result from differences in P plus “fallow benefits” to maize. Differences in maize yield among treatments for the -N subplot will result from differences in N plus P plus “fallow benefits” to maize.</p>
<p>The following steps were followed to generate the output in this document. The data was organized in excel rectangle columns with the different variables appearing in excel columns. All data checks were done in excel, meaningful data was selected and a copy of this data file was stored as a CSV file to make data import easy in R. The data file used in this analysis can be downloaded here: <a href="https://bit.ly/2rfLBEt" class="uri">https://bit.ly/2rfLBEt</a></p>
</div>
<div id="section-1-steps-in-analysis-using-r-1" class="section level2">
<h2><span class="header-section-number">4.2</span> Section 1: Steps in analysis using R</h2>
<ol style="list-style-type: decimal">
<li>Install R packages needed</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(ggplot2)
<span class="kw">library</span>(emmeans)
<span class="kw">library</span>(doBy)
<span class="kw">library</span>(lmerTest)
<span class="kw">library</span>(multcompView)</code></pre></div>
<ol start="2" style="list-style-type: decimal">
<li>Import data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fphosphorus <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/FPhosphorus.csv"</span>)</code></pre></div>
<ol start="3" style="list-style-type: decimal">
<li>Check and update data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fphosphorus)
<span class="kw">str</span>(fphosphorus)
fphosphorus<span class="op">$</span>mainplot<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>mainplot)
fphosphorus<span class="op">$</span>subplot<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>subplot)
fphosphorus<span class="op">$</span>block<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>block)</code></pre></div>
<ol start="4" style="list-style-type: decimal">
<li>Explore data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fphosphorus,<span class="kw">aes</span>(<span class="dt">y=</span>grain,<span class="dt">x=</span>fallow))<span class="op">+</span><span class="kw">geom_boxplot</span>(<span class="kw">aes</span>(<span class="dt">colour=</span>nitrogen))
<span class="kw">summaryBy</span>(grain<span class="op">~</span>fallow<span class="op">+</span>nitrogen, <span class="dt">data=</span>fphosphorus, <span class="dt">FUN=</span><span class="kw">c</span>(mean,sd))</code></pre></div>
<ol start="5" style="list-style-type: decimal">
<li>Specify a model for data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">splitplotmodel1<-<span class="kw">lmer</span>(grain<span class="op">~</span>fallow<span class="op">*</span>nitrogen<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>block<span class="op">/</span>mainplot), <span class="dt">data=</span>fphosphorus)</code></pre></div>
<ol start="6" style="list-style-type: decimal">
<li>Check the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(splitplotmodel1)
<span class="kw">qqnorm</span>(<span class="kw">resid</span>(splitplotmodel1))
<span class="kw">qqline</span>(<span class="kw">resid</span>(splitplotmodel1))
splitplotmodel2<-<span class="kw">lmer</span>(<span class="kw">sqrt</span>(grain)<span class="op">~</span>fallow<span class="op">*</span>nitrogen<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>block<span class="op">/</span>mainplot), <span class="dt">data=</span>fphosphorus)
<span class="kw">plot</span>(splitplotmodel2)
<span class="kw">qqnorm</span>(<span class="kw">resid</span>(splitplotmodel2))
<span class="kw">qqline</span>(<span class="kw">resid</span>(splitplotmodel2))</code></pre></div>
<ol start="7" style="list-style-type: decimal">
<li>Interpret the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">anova</span>(splitplotmodel2, <span class="dt">ddf=</span><span class="st">"Kenward-Roger"</span>)
<span class="kw">print</span>(<span class="kw">VarCorr</span>(splitplotmodel2), <span class="dt">comp=</span>(<span class="st">"Variance"</span>))</code></pre></div>
<ol start="8" style="list-style-type: decimal">
<li>Present the results from the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmip</span>(splitplotmodel2,nitrogen<span class="op">~</span>fallow,<span class="dt">CIs =</span> <span class="ot">TRUE</span>,<span class="dt">type=</span><span class="st">"response"</span>)
<span class="kw">emmeans</span>(splitplotmodel2,<span class="op">~</span>fallow,<span class="dt">type=</span><span class="st">"response"</span>)
<span class="kw">cld</span>(<span class="kw">emmeans</span>(splitplotmodel2,<span class="op">~</span>fallow,<span class="dt">type=</span><span class="st">"response"</span>))</code></pre></div>
</div>
<div id="section-2-explanation-of-steps-1" class="section level2">
<h2><span class="header-section-number">4.3</span> Section 2: Explanation of Steps</h2>
<div id="install-r-packages-needed-1" class="section level3">
<h3><span class="header-section-number">4.3.1</span> 1. Install R packages needed</h3>
<p>A number of packages following packages were used during data exploration and analysis. For a general introduction explaining what R packages are and how they work, this is a really useful guide <a href="https://www.datacamp.com/community/tutorials/r-packages-guide" class="uri">https://www.datacamp.com/community/tutorials/r-packages-guide</a>. For each of these packages to be installed, using install.packages(), this requires a reliable internet connection and a correctly installed version of R and RStudio. If you are having difficulties installing these packages please ask for help.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"ggplot2"</span>)
<span class="kw">library</span>(ggplot2)</code></pre></div>
<p><code>ggplot2</code> This package provides a powerful graphics language for creating elegant and complex graphs in R.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"emmeans"</span>)
<span class="kw">library</span>(emmeans)</code></pre></div>
<p><code>emmeans</code> Estimated marginal means (also known as least squares means) helps provide expected mean values and confidence intervals from statistical models.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"doBy"</span>)
<span class="kw">library</span>(doBy)</code></pre></div>
<p><code>doBy</code>Allows easy production of summary statistic tables</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"lmerTest"</span>)
<span class="kw">library</span>(lmerTest)</code></pre></div>
<p><code>lmerTest</code> Allows produce of flexible mixed effects regression models, similar to REML in Genstat.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"multcompView"</span>)
<span class="kw">library</span>(multcompView)</code></pre></div>
<p><code>multcompView</code> allows for mean seperation methods on analyses</p>
</div>
<div id="import-data-1" class="section level3">
<h3><span class="header-section-number">4.3.2</span> 2. Import data</h3>
<p>Our data set saved as a CSV file, so we can use the read.csv commmand to import the data. We are going to assign the name of the data with R to be <code>fphosphorus</code>. Remember in R Studio you could also use the “Import Dataset” menu to import a dataset.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fphosphorus <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/FPhosphorus.csv"</span>)</code></pre></div>
</div>
<div id="check-and-update-data-1" class="section level3">
<h3><span class="header-section-number">4.3.3</span> 3. Check and update data</h3>
<p>When reading data into R it is always useful to check that data is in the format expected. How many variables are there? How many rows? How have the columns been read in? The summary command can help to show if the data is being treated correctly.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fphosphorus)</code></pre></div>
<pre><code>## farmer block mainplot subplot
## BTS-10 : 8 Min. :1 Min. :1.00 Min. :1.0
## BTS-10B: 8 1st Qu.:3 1st Qu.:1.75 1st Qu.:1.0
## BTS-18 : 8 Median :5 Median :2.50 Median :1.5
## BTS-25 : 8 Mean :5 Mean :2.50 Mean :1.5
## BTS-27 : 8 3rd Qu.:7 3rd Qu.:3.25 3rd Qu.:2.0
## BTS-32 : 8 Max. :9 Max. :4.00 Max. :2.0
## (Other):24
## fallow nitrogen grain striga1
## Continous Maize:18 No :36 Min. :0.00 Min. : 7.00
## Crotolaria :18 Yes:36 1st Qu.:0.60 1st Qu.: 86.75
## Tephrosia :18 Median :1.25 Median : 323.50
## Tithonia :18 Mean :1.70 Mean : 582.04
## 3rd Qu.:2.40 3rd Qu.: 852.25
## Max. :5.30 Max. :2999.00
##
## striga2 striga3 striga4
## Min. : 2.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 13.75 1st Qu.: 10.0 1st Qu.: 17.0
## Median : 137.50 Median : 29.5 Median : 172.0
## Mean : 566.11 Mean : 84.4 Mean : 337.7
## 3rd Qu.: 926.00 3rd Qu.: 78.5 3rd Qu.: 448.5
## Max. :2645.00 Max. :1208.0 Max. :1406.0
## NA's :32</code></pre>
<p>Where data is being treated as a numeric variable (i.e. a number) <code>summary</code> provides statistics like the mean, min and max. Where data is being treated like a categorical variable (i.e. a group) then summary provides frequency tables.</p>
<p>From the results we can see that the variables block, mainplor and subplot are being considered as numeric variables. However these are grouping variables, not number variables, the numbers used are simply codes. If we do not rectify this then our analysis later will be incorrect and meaningless.<br />
This can also be seen more explicitly using the str() function.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">str</span>(fphosphorus)</code></pre></div>
<pre><code>## 'data.frame': 72 obs. of 11 variables:
## $ farmer : Factor w/ 9 levels "BTS-10","BTS-10B",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ block : int 7 7 7 7 7 7 7 7 8 8 ...
## $ mainplot: int 1 1 2 2 3 3 4 4 1 1 ...
## $ subplot : int 1 2 1 2 1 2 1 2 1 2 ...
## $ fallow : Factor w/ 4 levels "Continous Maize",..: 1 1 4 4 2 2 3 3 3 3 ...
## $ nitrogen: Factor w/ 2 levels "No","Yes": 1 2 2 1 2 1 2 1 1 2 ...
## $ grain : num 0.8 3 2.2 2.4 1.2 3 0.9 4.1 0.4 0.5 ...
## $ striga1 : int 1438 1340 482 340 98 90 232 120 2854 1715 ...
## $ striga2 : int 1736 960 2092 660 680 921 2645 1033 1709 941 ...
## $ striga3 : int 37 16 63 32 15 57 57 95 12 24 ...
## $ striga4 : int NA NA NA NA NA NA NA NA NA NA ...</code></pre>
<p>So we need to convert these variables into factors.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fphosphorus<span class="op">$</span>block<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>block)
fphosphorus<span class="op">$</span>mainplot<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>mainplot)
fphosphorus<span class="op">$</span>subplot<-<span class="kw">factor</span>(fphosphorus<span class="op">$</span>subplot)</code></pre></div>
<p>These commands take the column block within the data frame fphosphorus, converts into a factor and saves the result in a column called block within fphosphorus</p>
</div>
<div id="explore-data-1" class="section level3">
<h3><span class="header-section-number">4.3.4</span> 4. Explore data</h3>
</div>
<div id="plots-1" class="section level3">
<h3><span class="header-section-number">4.3.5</span> Plots</h3>
<p>In Tutorial 1 we produced plots showing all of the data plotted as points, like this:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fphosphorus,<span class="kw">aes</span>(<span class="dt">y=</span>grain,<span class="dt">x=</span>fallow,<span class="dt">colour=</span>nitrogen))<span class="op">+</span><span class="kw">geom_point</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-54-1.png" width="672" /> But in this instance there are too many points to be able to fully understand how the results are distributed. In this case we would get better information through looking at some boxplots.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fphosphorus,<span class="kw">aes</span>(<span class="dt">y=</span>grain,<span class="dt">x=</span>fallow,<span class="dt">colour=</span>nitrogen))<span class="op">+</span><span class="kw">geom_boxplot</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-55-1.png" width="672" /></p>
</div>
<div id="summary-statistics-1" class="section level3">
<h3><span class="header-section-number">4.3.6</span> Summary Statistics</h3>
<p>Using the summaryBy() function makes it easy to split summary statsitics into groups based on more than one factor. So the combination of fallow treatment and nitrogen treatment can be obtained using a + sign between the two variables.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summaryBy</span>(grain<span class="op">~</span>fallow<span class="op">+</span>nitrogen, <span class="dt">data=</span>fphosphorus, <span class="dt">FUN=</span><span class="kw">c</span>(mean,sd))</code></pre></div>
<pre><code>## fallow nitrogen grain.mean grain.sd
## 1 Continous Maize No 1.277778 1.446356
## 2 Continous Maize Yes 2.100000 1.948718
## 3 Crotolaria No 2.322222 1.475447
## 4 Crotolaria Yes 1.988889 1.334583
## 5 Tephrosia No 1.255556 1.181219
## 6 Tephrosia Yes 1.733333 1.024695
## 7 Tithonia No 1.255556 1.125956
## 8 Tithonia Yes 1.666667 1.736376</code></pre>
</div>
<div id="specify-a-model-for-data-1" class="section level3">
<h3><span class="header-section-number">4.3.7</span> 5. Specify a model for data</h3>
<p>In this design, a split plot desin, we have two treatment factors, “fallow” and “nitrogen”, and two layout factors “block” and “mainplot”.</p>
<p>In order to test the “main effects” of the treatmetns as well as the interaction between the two factors, then we need to make sure the formula is specified as factor1*factor2. Using factor1+factor2 will only include the main effects and not include the interaction.</p>
<p>When dealing with the split plot design, across multiple blocks, then the random effects need to be nested hierachically, from biggest down to smallest. This is done with a random effect that includes a / and looks like (1|biggestlayoutunit/nextbiggestlayout unit).</p>
<p>So the model we want to fit therefore looks like:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">splitplotmodel1<-<span class="kw">lmer</span>(grain<span class="op">~</span>fallow<span class="op">*</span>nitrogen<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>block<span class="op">/</span>mainplot), <span class="dt">data=</span>fphosphorus)</code></pre></div>
<p>R is unlike many other software packages in how it fits models. The best way of handling models in R is to assign the model to a name (in this case rcbdmodel1) and then ask R to provide different sorts of output for this model. When you run the above line you will get now output from the data - this is what we expected to see!</p>
</div>
<div id="check-the-model-1" class="section level3">
<h3><span class="header-section-number">4.3.8</span> 6. Check the model</h3>
<p>Before interpretting the model any further we should investigate the model validity, to ensure any conclusions we draw are valid. There are 3 assumptions that we can check for using standard model checking plots. 1. Homogeneity (equal variance) 2. Values with high leverage 3. Normality of residuals</p>
<p>The function plot() when used with a model will plot the fitted values from the model against the expected values.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(splitplotmodel1)</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-58-1.png" width="672" /> The residual Vs fitted plot is a scatter plot of the Residuals on the y-axis and the fitted on the x-axis and the aim for this plot is to test the assumption of equal variance of the residuals across the range of fitted values. There is someevidence of non-constant variance in our plot - residual values are less varaable around lower fitted values, and more variable around higher fitted values. This issue can often be solved by using a logaithmic or square root transformation. In this case, because there are some zero values within our data, it may be better to use a square root transformation.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">splitplotmodel2<-<span class="kw">lmer</span>(<span class="kw">sqrt</span>(grain)<span class="op">~</span>fallow<span class="op">*</span>nitrogen<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>block<span class="op">/</span>mainplot), <span class="dt">data=</span>fphosphorus)</code></pre></div>
<p>Refitting the plot shows a better approximation of heterogeneity, that is more acceptable to the assumptions required.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(splitplotmodel2)</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-60-1.png" width="672" /> We can also see that there are no extreme values in the residuals which might be potentially causing problems with the validity of our conclusions (leverage)</p>
<p>To assess the assumption of normality we can produce a qqplot. This shows us how closely the residuals follow a normal distribution - if there are severe and syste,matic deviations from the line then we may want to consider an alternative distribution.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">qqnorm</span>(<span class="kw">resid</span>(splitplotmodel2))
<span class="kw">qqline</span>(<span class="kw">resid</span>(splitplotmodel2))</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-61-1.png" width="672" /> In this case the residuals seem to fit the assumption required for normality.</p>
</div>
<div id="interpret-model-1" class="section level3">
<h3><span class="header-section-number">4.3.9</span> 7. Interpret Model</h3>
<p>The anova() function only prints the rows of analysis of variance table for treatment effects when looking at a mixed model fitted using lmer().</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">anova</span>(splitplotmodel2,<span class="dt">ddf=</span><span class="st">"Kenward-Roger"</span>)</code></pre></div>
<pre><code>## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## fallow 0.38006 0.12669 3 24 1.0112 0.4050
## nitrogen 0.27129 0.27129 1 32 2.1653 0.1509
## fallow:nitrogen 0.37536 0.12512 3 32 0.9987 0.4061</code></pre>
<p>ddf=Kenward-Roger tells R which method to use for determining the calculations of the table; this option matches the defaults found within SAS or Genstat. The ANOVA table suggests a highly significant effect of the treatment on the yield.</p>
<p>To obtain the residual variance, and the variance attributed to the blocks we need an additional command. From these number it is possible to reconstruct a more classic ANOVA table, if so desired.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">print</span>(<span class="kw">VarCorr</span>(splitplotmodel1), <span class="dt">comp=</span>(<span class="st">"Variance"</span>))</code></pre></div>
<pre><code>## Groups Name Variance
## mainplot:block (Intercept) 0.38238
## block (Intercept) 0.64581
## Residual 1.04375</code></pre>
</div>
<div id="present-the-results-from-the-model-1" class="section level3">
<h3><span class="header-section-number">4.3.10</span> 8. Present the results from the model</h3>
<p>To help understand what the significant result from the ANOVA table means we can produce several plots and tables to help us. First we can use the function emmip() to produce plots of the modelled results, including 95% confidence intervals.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmip</span>(splitplotmodel2,nitrogen<span class="op">~</span>fallow,<span class="dt">CIs =</span> <span class="ot">TRUE</span>,<span class="dt">type =</span> <span class="st">"response"</span>)</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-64-1.png" width="672" /></p>
<p>To obtain the numbers used in creating this graph we can use the function emmeans.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmeans</span>(splitplotmodel1, <span class="op">~</span>fallow)</code></pre></div>
<pre><code>## NOTE: Results may be misleading due to involvement in interactions</code></pre>
<pre><code>## fallow emmean SE df lower.CL upper.CL
## Continous Maize 1.688889 0.415006 21.04 0.8259423 2.551835
## Crotolaria 2.155556 0.415006 21.04 1.2926090 3.018502
## Tephrosia 1.494444 0.415006 21.04 0.6314978 2.357391
## Tithonia 1.461111 0.415006 21.04 0.5981645 2.324058
##
## Results are averaged over the levels of: nitrogen
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95</code></pre>
<p>And one method for conducting mean separation analysis we can use the function cld().</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">cld</span>(<span class="kw">emmeans</span>(splitplotmodel1, <span class="op">~</span>fallow))</code></pre></div>
<pre><code>## NOTE: Results may be misleading due to involvement in interactions</code></pre>
<pre><code>## fallow emmean SE df lower.CL upper.CL .group
## Tithonia 1.461111 0.415006 21.04 0.5981645 2.324058 1
## Tephrosia 1.494444 0.415006 21.04 0.6314978 2.357391 1
## Continous Maize 1.688889 0.415006 21.04 0.8259423 2.551835 1
## Crotolaria 2.155556 0.415006 21.04 1.2926090 3.018502 1
##
## Results are averaged over the levels of: nitrogen
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05</code></pre>
<p>In the output, groups sharing a letter in the .group are not statistically different from each other.</p>
</div>
</div>
<div id="section-3-methodological-principles-1" class="section level2">
<h2><span class="header-section-number">4.4</span> Section 3 – Methodological Principles</h2>
<p>There are always many different ways of doing all that we have done here in R. The less complex the method/code is, the better it is for you so that you can easily grasp the method.</p>
<p>In this example using the phosphurus data, we have a split plot design. This means that a single plot where the fallow treatment has been applied is split into 2, and each half receives a different nitrogen treatment. It is useful to have seperate columns denoting treamtent factors and layout factors - even if these may be somewhat replicating the same information. The split plots are nested within the plots, which are nested within the blocks. So the random effect needs to incorporate this nesting. Remember that the lowest level design factor, the split plot, does not get included in the model. This is similar to the RCBD analysis, where the lowest level factor - plot, does not get included in the model.</p>
<p>Note that the difference in the specification of random effects in the model is effectively the only difference needed in the R syntax used to produce this analysis, as compared to Tutorial 1, the RCBD. All other syntax has been modified to reflect differences in the data collected, but the same functions (ggplot, summaryBy, emmeans) are being applied in the same way.</p>
<p>Food for thought: Your best model will certainly be as good as the data you collected!!!</p>
<!--chapter:end:03-splitplot.Rmd-->
</div>
</div>
<div id="adjusting-for-covariates" class="section level1">
<h1><span class="header-section-number">5</span> Adjusting for Covariates</h1>
<div id="about-the-data-2" class="section level2">
<h2><span class="header-section-number">5.1</span> About the data</h2>
<p>The data used in this example is from a study was conducted in Eastern Zambia and the main aim was to improve on the efficiency of the natural fallows by using appropriate trees that may have relevance in soil fertility regeneration within permissible fallow periods. This is the same data used in the first part of this series.</p>
<p>The design was a randomized complete block design experiment with 4 blocks and 9 treatments was conducted. The primary outcome variable was crop yield (yield). We also have data collected on striga infestation.</p>
<p>The objective for this analysis is to investigate the relationship between striga infestation and yield across the different treatments.</p>
<p>The following steps were followed to generate the output in this document. The data was organized in excel rectangle columns with the different variables appearing in excel columns. All data checks were done in excel, meaningful data was selected and a copy of this data file was stored as a CSV file to make data import easy in R. The data file used in this analysis can be downloaded here: <a href="https://bit.ly/2rfLBEt" class="uri">https://bit.ly/2rfLBEt</a></p>
</div>
<div id="section-1-steps-in-analysis-using-r-2" class="section level2">
<h2><span class="header-section-number">5.2</span> Section 1: Steps in analysis using R</h2>
<ol style="list-style-type: decimal">
<li>Install R packages needed</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(ggplot2)
<span class="kw">library</span>(emmeans)
<span class="kw">library</span>(doBy)
<span class="kw">library</span>(lmerTest)
<span class="kw">library</span>(multcompView)</code></pre></div>
<ol start="2" style="list-style-type: decimal">
<li>Import data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/Fallow N2.csv"</span>)</code></pre></div>
<ol start="3" style="list-style-type: decimal">
<li>Check and update data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fallow)
<span class="kw">str</span>(fallow)
fallow<span class="op">$</span>rep<-<span class="kw">factor</span>(fallow<span class="op">$</span>rep)
fallow<span class="op">$</span>plot<-<span class="kw">factor</span>(fallow<span class="op">$</span>plot)</code></pre></div>
<ol start="4" style="list-style-type: decimal">
<li>Explore data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>treat,<span class="dt">col=</span>rep))<span class="op">+</span><span class="kw">geom_point</span>()
<span class="kw">summaryBy</span>(yield<span class="op">~</span>treat, <span class="dt">data=</span>fallow, <span class="dt">FUN=</span><span class="kw">c</span>(min,max,mean,median,sd))</code></pre></div>
<ol start="5" style="list-style-type: decimal">
<li>Specify a model for data</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">rcbdmodel1<-<span class="kw">lmer</span>(yield<span class="op">~</span>treat<span class="op">+</span>(<span class="dv">1</span><span class="op">|</span>rep),<span class="dt">data=</span>fallow)</code></pre></div>
<ol start="6" style="list-style-type: decimal">
<li>Check the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(rcbdmodel1)
<span class="kw">qqnorm</span>(<span class="kw">resid</span>(rcbdmodel1))
<span class="kw">qqline</span>(<span class="kw">resid</span>(rcbdmodel1))</code></pre></div>
<ol start="7" style="list-style-type: decimal">
<li>Interpret the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">anova</span>(rcbdmodel1,<span class="dt">ddf=</span><span class="st">"Kenward-Roger"</span>)
<span class="kw">print</span>(<span class="kw">VarCorr</span>(rcbdmodel1), <span class="dt">comp=</span>(<span class="st">"Variance"</span>))</code></pre></div>
<ol start="8" style="list-style-type: decimal">
<li>Present the results from the model</li>
</ol>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">emmip</span>(rcbdmodel1,<span class="op">~</span>treat,<span class="dt">CIs =</span> <span class="ot">TRUE</span>)
<span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat)
<span class="kw">cld</span>(<span class="kw">emmeans</span>(rcbdmodel1, <span class="op">~</span>treat))</code></pre></div>
</div>
<div id="section-2-explanation-of-steps-2" class="section level2">
<h2><span class="header-section-number">5.3</span> Section 2: Explanation of Steps</h2>
<div id="install-r-packages-needed-2" class="section level3">
<h3><span class="header-section-number">5.3.1</span> 1. Install R packages needed</h3>
<p>A number of packages following packages were used during data exploration and analysis. For a general introduction explaining what R packages are and how they work, this is a really useful guide <a href="https://www.datacamp.com/community/tutorials/r-packages-guide" class="uri">https://www.datacamp.com/community/tutorials/r-packages-guide</a>. For each of these packages to be installed, using install.packages(), this requires a reliable internet connection and a correctly installed version of R and RStudio. If you are having difficulties installing these packages please ask for help.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"ggplot2"</span>)
<span class="kw">library</span>(ggplot2)</code></pre></div>
<p><code>ggplot2</code> This package provides a powerful graphics language for creating elegant and complex graphs in R.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"emmeans"</span>)
<span class="kw">library</span>(emmeans)</code></pre></div>
<p><code>emmeans</code> Estimated marginal means (also known as least squares means) helps provide expected mean values and confidence intervals from statistical models.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"doBy"</span>)
<span class="kw">library</span>(doBy)</code></pre></div>
<p><code>doBy</code>Allows easy production of summary statistic tables</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"lmerTest"</span>)
<span class="kw">library</span>(lmerTest)</code></pre></div>
<p><code>lmerTest</code> Allows produce of flexible mixed effects regression models, similar to REML in Genstat.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"multcompView"</span>)
<span class="kw">library</span>(multcompView)</code></pre></div>
<p><code>multcompView</code> allows for mean seperation methods on analyses</p>
</div>
<div id="import-data-2" class="section level3">
<h3><span class="header-section-number">5.3.2</span> 2. Import data</h3>
<p>Our data set saved as a CSV file, so we can use the read.csv commmand to import the data. We are going to assign the name of the data with R to be <code>fallow2</code>. Remember in R Studio you could also use the “Import Dataset” menu to import a dataset.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow <-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">"C:/Users/Admin/Desktop/Fallow N2.csv"</span>)</code></pre></div>
</div>
<div id="check-and-update-data-2" class="section level3">
<h3><span class="header-section-number">5.3.3</span> 3. Check and update data</h3>
<p>When reading data into R it is always useful to check that data is in the format expected. How many variables are there? How many rows? How have the columns been read in? The summary command can help to show if the data is being treated correctly.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(fallow)</code></pre></div>
<pre><code>## rep plot treat yield
## Min. :1.00 Min. :1 1 S.sesban : 4 Min. :1.140
## 1st Qu.:1.75 1st Qu.:3 2 G.sepium : 4 1st Qu.:2.370
## Median :2.50 Median :5 3 L.leuco : 4 Median :3.140
## Mean :2.50 Mean :5 4 F.congesta: 4 Mean :3.232
## 3rd Qu.:3.25 3rd Qu.:7 5 C.siamea : 4 3rd Qu.:3.728
## Max. :4.00 Max. :9 6 C.calo : 4 Max. :6.540
## (Other) :12
## striga
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 21.0
## Mean : 334.1
## 3rd Qu.: 238.5
## Max. :2798.0
## </code></pre>
<p>Where data is being treated as a numeric variable (i.e. a number) <code>summary</code> provides statistics like the mean, min and max. Where data is being treated like a categorical variable (i.e. a group) then summary provides frequency tables.</p>
<p>From the results we can see that the variables rep and plot are being considered as numeric variables. However these are grouping variables, not number variables, the numbers used are simply codes. If we do not rectify this then our analysis later will be incorrect and meaningless.<br />
This can also be seen more explicitly using the str() function.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">str</span>(fallow)</code></pre></div>
<pre><code>## 'data.frame': 36 obs. of 5 variables:
## $ rep : int 1 4 4 1 1 3 3 1 3 2 ...
## $ plot : int 2 3 6 9 7 3 8 6 9 9 ...
## $ treat : Factor w/ 9 levels "1 S.sesban","2 G.sepium",..: 8 5 8 7 5 8 5 9 6 5 ...
## $ yield : num 1.14 1.74 1.95 2.06 2.09 2.15 2.21 2.22 2.34 2.38 ...
## $ striga: int 2798 0 1787 129 1 1144 0 228 0 0 ...</code></pre>
<p>So we need to convert these variables into factors.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fallow<span class="op">$</span>rep<-<span class="kw">factor</span>(fallow<span class="op">$</span>rep)
fallow<span class="op">$</span>plot<-<span class="kw">factor</span>(fallow<span class="op">$</span>plot)</code></pre></div>
<p>These commands take the column rep within the data frame fallow, converts into a factor and saves the result in a column called rep within fallow.</p>
</div>
<div id="explore-data-2" class="section level3">
<h3><span class="header-section-number">5.3.4</span> 4. Explore data</h3>
<div id="plots-2" class="section level4">
<h4><span class="header-section-number">5.3.4.1</span> Plots</h4>
<p>We are now interesting in assessing the relationship between yield and striga - so we want to produce a plot of striga against yield, with different coloured points denoting each treatment.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span>striga,<span class="dt">col=</span>treat))<span class="op">+</span><span class="kw">geom_point</span>()</code></pre></div>
<p><img src="IntroToAgriculturalExperimentsinR_files/figure-html/unnamed-chunk-86-1.png" width="672" /> We can see from the distribution of striga that there are some farms with very high levels of striga, and some farms with no striga. The big range of values makes it hard to make interpretations from this plot, so taking a square root transformation may help to visualise the relationship. A log transformation will not help here because of the large number of 0 values of striga.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ggplot</span>(<span class="dt">data=</span>fallow,<span class="kw">aes</span>(<span class="dt">y=</span>yield,<span class="dt">x=</span><span class="kw">sqrt</span>(striga),<span class="dt">col=</span>treat))<span class="op">+</span><span class="kw">geom_point</span>()</code></pre></div>