0

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!

camille
  • 16,432
  • 18
  • 38
  • 60
M F
  • 11
  • 2
  • The example code you've provided does not seem to run successfully, so we can't see the results of the model or the ANOVA. – jdobres May 18 '21 at 21:52
  • There are bunch of missing pieces here. Your reprexes should include `library()` calls to `lme4`, `car`, `dplyr`, `gtsummary`, and `sjPlot`, otherwise they won't work and provide the right information. Your fundamental question is one about how linear models are parameterized and displayed in R. Since you can't estimate *independent* effects of (e.g.) all three species - if you know the mean and what two species are doing, then the third effect is determined - R only shows you *differences* from a baseline rate ... – Ben Bolker May 18 '21 at 22:14
  • I think this probably deserves to get migrated to [CrossValidated](https://stats.stackexchange.com), but you should try to sort out the problem with your reprexes so we can see the results. – Ben Bolker May 18 '21 at 22:14
  • Thank you for your comments, i've finally made a reprex that works! Buut now i'm starting to understand this is not a problem, rather a misunderstanding on my part, about how lmer and more generally, mixed effect models, are computed, as you mentioned above. So the baseline in this case, would be species A and SedLayer 1, right? It gets very confusing when it comes to interpreting interactions though. – M F May 19 '21 at 10:47
  • I don't want to accidentally delete something that you need here, but the first several lines of code seem to be for creating the reprex, not stuff that we need in order to run your code (loading `datapasta`, calling `dput`, etc) and just creates errors. Aside from that, this is a question for [stats.se] since the question is about interpretation rather than code debugging, but I'd recommend cleaning up that beginning section – camille May 20 '21 at 00:46
  • Yes I was scared to delete anything as well in case it affected the reprex. I'll ask the question in CrossValidated then, I had no idea there was a difference, I thought both were help forums. – M F May 20 '21 at 16:27

0 Answers0