I have a study design in which I planted three different plants (S1 - S3) with two different treatments (T1, T2) in five blocks (B1-B5). I sampled the soil microbes under 'em in two timepoints (Y1-Y2).
I want to find out plantstreatmentyear effect on microbial community composition (all three-way and two-way interactions are important to me). The way I understand it, I need to include both blocks and plots (subjects of repeated measures; block x plant x treatment combos) as random effects (first as blocking and latter for repeated measure).
After much sleuthing, I found permute package and related online questions which almost got me to where I need to be. I think I am structuring my model right, but Adonis2 is not cooperating. Here is my code:
library(vegan)
library(tidyr)
library(dplyr)
library(tibble) #i think these would do.. if not, I apologize!
meta.raw<- as.factor(c(1:60))
meta.raw %>%
as.data.frame() %>%
`colnames<-`("sample") %>%
mutate(block = rep(paste0("B", sample(1:5)), 2, each = 6),
crop = rep(c("S1","S2","S3"), 10, each = 2),
trt = rep (c("T1", "T2"), 30),
year = rep(c("Y1", "Y2"), 2, each = 15),
plot = paste0(block,crop,trt)) -> meta
matrix(round(runif(n=6000, min=0, max=2000), 0), nrow=60) %>%
as.data.frame() %>%
mutate(sample = as.factor(c(1:60))) %>%
column_to_rownames("sample")-> df
# how1: Syntax as it appears on "permute" vignette, (as far as I can tell)
h1 <- how(within = Within(type = "series"),
plots = Plots(strata = meta$plot),
blocks = Blocks(strata = meta$block),
nperm = 499)
adonis2(df ~ crop*trt*year,
data = meta,
dist = "bray",
perm = h1,
by = "margin")
# this doesn't work:
# Error in check(sn, control = control, quietly = quietly) :
# Number of observations and length of Block 'strata' do not match.
# how 2: passing within and plots inside how() and passing blocking to strata = argument in adonis2()
h2 <- how(within = Within(type = "series"),
plots = Plots(strata = meta$plot),
nperm = 499) # no blocks here
adonis2(df ~ crop*trt*year,
data = meta,
dist = "bray",
perm = h2,
strata = meta$block, # block here instead
by = "margin")
# this works... but is it "correct"?
# Permutation test for adonis under reduced model
# Marginal effects of terms
# Blocks: strata
# Plots: meta$plot, plot permutation: none
# Permutation: series
# Number of permutations: 499
#
# adonis2(formula = df ~ crop * trt * year, data = meta, permutations = h2, by = "margin", strata = meta$block, dist = "bray")
# Df SumOfSqs R2 F Pr(>F)
# crop:trt:year 2 0.1027 0.03085 0.9117 1
# Residual 48 2.7043 0.81219
# Total 59 3.3296 1.00000
#how 3: passing blocking vairable within how() but not in the same syntax as it appears in the "permute" vingette
h3 <- how(within = Within(type = "series"),
plots = Plots(strata = meta$plot),
blocks = meta$block,
nperm = 499)
adonis2(df ~ crop*trt*year,
data = meta,
dist = "bray",
perm = h3,
by = "margin")
# this also works???
# Permutation test for adonis under reduced model
# Marginal effects of terms
# Blocks: meta$block
# Plots: meta$plot, plot permutation: none
# Permutation: series
# Number of permutations: 499
#
# adonis2(formula = df ~ crop * trt * year, data = meta, permutations = h3, by = "margin", dist = "bray")
# Df SumOfSqs R2 F Pr(>F)
# crop:trt:year 2 0.1007 0.02936 0.8709 1
# Residual 48 2.7760 0.80897
# Total 59 3.4316 1.00000
What am I doing wrong (or right)? I feel like the first option should be the right one, but it is obviously not working... The same error messages have shown up for folks but it seems to be just for unbalanced data (maybe), whereas mine is (at least the made up one).
Thank you all in advance :)
The first option, using a Blocks() function, is wrong. I know there is a Blocks() function in permute but it isn't used and I can't recall now ever using this. Perhaps I thought it might be neater to set the blocks in a design using a function like you set the within plot and plot level restrictions. But I have not made this function work and I will most likely remove this function, or make it work as it needs to.
You approach 2 works, because adonis2() is setting the blocks element on the thing you passed to perm internally. It is updating the how(), but now you have two hows and I thin that is not desirable.
Approach three is the preferred way of doing this. Actually, I personally prefer to do it with with(), to avoid having to use meta$ repeatedly to refer to columns:
h3 <- with(
how(within = Within(type = "series"),
plots = Plots(strata = plot),
blocks = block,
nperm = 499)
)
If you saw me using Blocks() anywhere in the documentation etc, please let me know where so I can correct it.