I know that pretty much the same question has been asked here 'list' object cannot be coerced to type 'integer' -- WRS2 package
Unfortunately, I cannot comment yet and ask more questions regarding the solution of the problem.
I am using the pbad2way function of WRS2-package. R version 4.2.2 (2022-10-31 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19044
I want to calculate robust ANOVA as explained on https://rcompanion.org/rcompanion/d_08a.html . I am getting the error:
"Error in sequence.default(piece, ...) : 'list' object cannot be coerced to type 'integer'"
When using the following code:
pbad2way(avg_rad ~ antifungal + dose + antifungal:dose, data = y2, est = "mom", nboot = 5000)
My data.frame looks like this:
structure(list(avg_rad = c(9.91e+08, 8.17e+08, 461200000, 15330000,
175100000, 50320000, 13590000, 22970000, 2778000, 3453000, 12890000,
375900000, 44590000, 1.611e+09, 1e+09, 889900000, 373200000,
5010000, 6549000, 23160000, 32520000, 7707000, 556900000, 634600000,
820900000, 391400000, 498300000, 147900000, 646900000, 22060000,
1e+07, 306800000, 319400000, 41290000, 94100000, 127200000, 117200000,
618300000, 570700000, 617100000, 284900000, 449600000, 3866000,
6918000, 4177000, 14870000, 29380000, 2815000, 1619000, 3126000,
1710000, 2191000), antifungal = structure(c(1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L,
3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels = c("untreated",
"VRZ", "AMB"), class = "factor"), dose = structure(c(1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L,
3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L,
3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L,
3L, 3L), levels = c("0", "5", "20"), class = "factor")), row.names = c(NA,
-52L), class = c("tbl_df", "tbl", "data.frame"))
Is my assumption correct that it does not work because I have cero cells for some combinations of my factor dose?
with(y2, table(antifungal, dose))
| antifungal | dose 0 | dose 5 | dose 20 |
|---|---|---|---|
| untreated | 17 | 0 | 0 |
| AMB | 0 | 10 | 10 |
| VRZ | 0 | 5 | 10 |
Is there a way to solve this problem or to improve my data.frame in its structure from the beginning?
I.
I suspect you're correct, that due to the unbalanced design, the pbad2way() function won't accept it.
II.
There's a question of how to analyze the data given the nature of the dependent variable. There's some heteroscedasticity and skewness to consider.
In this case Gamma regression (with glm() and the Gamma() family appears to be appropriate.
I think simply taking the log of the dependent variable works fine in this case, and leads to relatively simple interpretation. Essentially, you would be doing an analysis of variance on the geometric means. I'll use this approach in the example below. But you could use Gamma regression with almost the same codes. (But you might then modify the emmeans call to marginal = emmeans(model, ~ Treat, type="response").
III.
There's a question about to handle this design, which is essentially a 2 x 2 factorial plus a control treatment. Traditionally, this can be approached by using a one-way design of all treatments, and then comparing whatever contrasts are of interest. Below, you could also use pairs(marginal) to see all the pairwise comparisons. I have some additional examples here: rcompanion.org/rcompanion/h_01.html .
y2 = structure(list(avg_rad = c(9.91e+08, 8.17e+08, 461200000, 15330000,
175100000, 50320000, 13590000, 22970000, 2778000, 3453000, 12890000,
375900000, 44590000, 1.611e+09, 1e+09, 889900000, 373200000,
5010000, 6549000, 23160000, 32520000, 7707000, 556900000, 634600000,
820900000, 391400000, 498300000, 147900000, 646900000, 22060000,
1e+07, 306800000, 319400000, 41290000, 94100000, 127200000, 117200000,
618300000, 570700000, 617100000, 284900000, 449600000, 3866000,
6918000, 4177000, 14870000, 29380000, 2815000, 1619000, 3126000,
1710000, 2191000), antifungal = structure(c(1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L,
3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels = c("untreated",
"VRZ", "AMB"), class = "factor"), dose = structure(c(1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L,
3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L,
3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L,
3L, 3L), levels = c("0", "5", "20"), class = "factor")), row.names = c(NA,
-52L), class = c("tbl_df", "tbl", "data.frame"))
y2$Treat = interaction(y2$antifungal, y2$dose)
y2$Treat = droplevels(y2$Treat)
levels(y2$Treat)[1]="Untreated"
y2$lograd = log(y2$avg_rad)
plot(lograd ~ Treat, data=y2)
model = lm(lograd ~ Treat, data = y2)
library(car)
Anova(model)
hist(residuals(model))
plot(residuals(model) ~ predict(model))
library(emmeans)
marginal = emmeans(model, ~ Treat)
marginal
pairs(marginal)
levels(y2$Treat)
### "Untreated" "VRZ.5" "AMB.5" "VRZ.20" "AMB.20"
Contrasts = list(UntreatedvsRest = c(-4, 1, 1, 1, 1),
VRZvsAMB = c( 0, 1, -1, 1, -1),
c5vsc20 = c( 0, 1, 1, -1, -1),
InteractionAD = c( 0, 1, -1, -1, 1))
contrast(marginal, Contrasts, adjust="sidak")
### contrast estimate SE df t.ratio p.value
### UntreatedvsRest -14.128 1.379 47 -10.245 <.0001
### VRZvsAMB 4.061 0.812 47 5.000 <.0001
### c5vsc20 0.742 0.812 47 0.913 0.8381
### InteractionAD -0.349 0.812 47 -0.430 0.9880
###
### P value adjustment: sidak method for 4 tests