Importing from Rmarkdown

Mixed model ANOVA

Very straightforward.

#load packages
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.2     ✓ dplyr   1.0.6
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

flum <- read_csv("http://rstats4ag.org/data/FlumiBeans.csv")

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   year = col_double(),
##   treatment = col_character(),
##   block = col_double(),
##   population.4wk = col_double()
## )

head(flum)

## # A tibble: 6 x 4
##    year treatment                 block population.4wk
##   <dbl> <chr>                     <dbl>          <dbl>
## 1  2009 Nontreated                    1          55757
## 2  2009 Nontreated                    2          45302
## 3  2009 Nontreated                    3          38333
## 4  2009 flumioxazin + trifluralin     1          13939
## 5  2009 flumioxazin + trifluralin     2          27878
## 6  2009 flumioxazin + trifluralin     3          31363

str(flum)

## spec_tbl_df [77 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ year          : num [1:77] 2009 2009 2009 2009 2009 ...
##  $ treatment     : chr [1:77] "Nontreated" "Nontreated" "Nontreated" "flumioxazin + trifluralin" ...
##  $ block         : num [1:77] 1 2 3 1 2 3 1 2 3 1 ...
##  $ population.4wk: num [1:77] 55757 45302 38333 13939 27878 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   year = col_double(),
##   ..   treatment = col_character(),
##   ..   block = col_double(),
##   ..   population.4wk = col_double()
##   .. )

flum <- flum %>%
  mutate(yrblock = with(flum, factor(year):factor(block)))
head(flum)

## # A tibble: 6 x 5
##    year treatment                 block population.4wk yrblock
##   <dbl> <chr>                     <dbl>          <dbl> <fct>  
## 1  2009 Nontreated                    1          55757 2009:1 
## 2  2009 Nontreated                    2          45302 2009:2 
## 3  2009 Nontreated                    3          38333 2009:3 
## 4  2009 flumioxazin + trifluralin     1          13939 2009:1 
## 5  2009 flumioxazin + trifluralin     2          27878 2009:2 
## 6  2009 flumioxazin + trifluralin     3          31363 2009:3