Reputation: 11
For a bit of context: I am trying to see the effect of 3 different mangrove species (A, L and R) on sediment characteristics, in this case the dependent variable is pH. I have 10 sites, which have been coded as random effect.
SedLayer indicates the top layer (1) and the bottom layer (9) of sediment collected over the sites, near each of the three species. Now when I run my lmer code, I get a weird result where only two species show up, and similarly for sediment layer, only layer 9 is shown.
When I run an anova on my object (respH), I get "understandable" results, showing Species, SedLayer and Species:SedLayer. However, when it comes to building my table, the anova results are not accepted. I have tried with sjPlots and gtsummary, both packages refusing to use the anova results, and working only with my lmer result (respH).
library("car")
library("lme4")
library("reprex")
library("dplyr")
library("sjPlot")
library("datapasta")
dat <- dput(Dataset_1_1 %>% select(2:5,))
#> Error in select(., 2:5, ): object 'Dataset_1_1' not found
dat$Species <- as.factor(dat$Species)
#> Error in is.factor(x): object 'dat' not found
dat$`Sed. Layer` <- as.factor(dat$`Sed. Layer`)
#> Error in is.factor(x): object 'dat' not found
dat$Site <- as.factor(dat$Site)
#> Error in is.factor(x): object 'dat' not found
str(dat)
#> Error in str(dat): object 'dat' not found
colnames(dat) <- c("Species", "Site", "SedLayer", "pH")
#> Error in colnames(dat) <- c("Species", "Site", "SedLayer", "pH"): object 'dat' not found
dat <- tibble::tribble(
~Species, ~Site, ~SedLayer, ~pH,
"A", "1", "1", 6.361,
"A", "1", "9", 6.602,
"A", "2", "1", 5.83,
"A", "2", "9", 6.053,
"A", "3", "1", 7.03,
"A", "3", "9", 6.65,
"A", "4", "1", 6.579,
"A", "4", "9", 6.436,
"A", "5", "1", 6.411,
"A", "5", "9", 5.962,
"A", "6", "1", 6.914,
"A", "6", "9", 6.693,
"A", "7", "1", 5.987,
"A", "7", "9", 5.717,
"A", "8", "1", 5.606,
"A", "8", "9", 5.28,
"A", "9", "1", 6.783,
"A", "9", "9", 6.922,
"A", "10", "1", 7.703,
"A", "10", "9", 6.251,
"L", "1", "1", 6.55,
"L", "1", "9", 6.705,
"L", "3", "1", 7.184,
"L", "3", "9", 6.814,
"L", "4", "1", 6.552,
"L", "4", "9", 6.555,
"L", "5", "1", 6.638,
"L", "5", "9", 5.943,
"L", "6", "1", 6.32,
"L", "6", "9", 6.738,
"L", "7", "1", 6.313,
"L", "7", "9", 5.971,
"L", "8", "1", 4.963,
"L", "8", "9", 5.274,
"L", "9", "1", 6.443,
"L", "9", "9", 6.293,
"L", "10", "1", 7.247,
"L", "10", "9", 6.25,
"R", "1", "1", 5.443,
"R", "1", "9", 4.89,
"R", "2", "1", 6.055,
"R", "2", "9", 4.74,
"R", "3", "1", 6.731,
"R", "3", "9", 6.563,
"R", "4", "1", 6.303,
"R", "4", "9", 6.264,
"R", "5", "1", 5.71,
"R", "5", "9", 6.087,
"R", "6", "1", 4.045,
"R", "6", "9", 5.026,
"R", "7", "1", 5.835,
"R", "7", "9", 5.375,
"R", "8", "1", 5.332,
"R", "8", "9", 5.579,
"R", "9", "1", 6.141,
"R", "9", "9", 6.508,
"R", "10", "1", 7.354,
"R", "10", "9", 6.783
)
#lme4
respH <- lmer(pH ~ SedLayer * Species + (+1|Site), data = dat)
summary(respH)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: pH ~ SedLayer * Species + (+1 | Site)
#> Data: dat
#>
#> REML criterion at convergence: 103.8
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.3479 -0.5001 -0.0156 0.5340 1.6956
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> Site (Intercept) 0.2054 0.4532
#> Residual 0.2448 0.4947
#> Number of obs: 58, groups: Site, 10
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 6.52040 0.21217 30.731
#> SedLayer9 -0.26380 0.22126 -1.192
#> SpeciesL -0.09062 0.22847 -0.397
#> SpeciesR -0.62550 0.22126 -2.827
#> SedLayer9:SpeciesL 0.07858 0.32148 0.244
#> SedLayer9:SpeciesR 0.15040 0.31291 0.481
#>
#> Correlation of Fixed Effects:
#> (Intr) SdLyr9 SpecsL SpecsR SL9:SL
#> SedLayer9 -0.521
#> SpeciesL -0.505 0.484
#> SpeciesR -0.521 0.500 0.484
#> SdLyr9:SpcL 0.359 -0.688 -0.704 -0.344
#> SdLyr9:SpcR 0.369 -0.707 -0.342 -0.707 0.487
AnorespH <- Anova(respH)
AnorespH
#> Analysis of Deviance Table (Type II Wald chisquare tests)
#>
#> Response: pH
#> Chisq Df Pr(>Chisq)
#> SedLayer 2.0837 1 0.1488746
#> Species 14.8467 2 0.0005972 ***
#> SedLayer:Species 0.2312 2 0.8908424
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sjPlot
tab_model(respH, show.df = TRUE)
I'm sorry if this is not clear, I'm new to complex stats, new to R and new to this forum.
Here is a screenshot of the table as it's not showing up properly with the reprex sjPlot table
Thanks in advance for any input you can give me!
Upvotes: 0
Views: 404