Plant Height Analysis: Species as Fixed Effect Nested within Nativity
Your Name
2026-02-05
Introduction
This analysis examines whether native vs non-native plants differ in height change using a model where:
- Species is treated as a FIXED EFFECT nested within Nativity
- Plot remains a random effect for spatial clustering
Important Conceptual Difference:
- Previous model: Species as random effect (generalizes to other species)
- This model: Species as fixed effect (specific to these species only)
Step 1: Load Required Packages
library(readxl) # Read Excel files
library(lme4) # Fit Mixed-Effects Models
library(lmerTest) # Get p-values for mixed models
library(ggplot2) # Create plots
library(dplyr) # Data wrangling
library(car) # ANOVA and assumption tests
library(emmeans) # Estimated marginal means
library(multcomp) # Multiple comparisonsStep 2: Import and Clean Data
getwd()## [1] "/Users/ivyvelasco/Downloads"setwd("/Users/ivyvelasco/Desktop/Greenhouse 2026")
getwd()## [1] "/Users/ivyvelasco/Desktop/Greenhouse 2026"list.files()## [1] "arthropod_barplot.png"
## [2] "arthropod_boxplot.png"
## [3] "Arthropods_R.xlsx"
## [4] "Book5.xlsx"
## [5] "CarbonFluxLMM.R"
## [6] "CFlux_Anova#1.Rmd"
## [7] "CFlux_R.xlsx"
## [8] "emmeans_Nativity_by_Light_Edited.xlsx"
## [9] "emmeans_Nativity_by_Light.xlsx"
## [10] "emmeans_Nativity_Light.xlsx"
## [11] "emmeans&contrasts_Light_by_Nativity_Edited.xlsx"
## [12] "emmeans&contrasts_PairwiseComparison_Edited.xlsx"
## [13] "height_change_barplot.png"
## [14] "height_change_boxplot.png"
## [15] "height_change_by_species.png"
## [16] "HeightDataCleaning.xlsx"
## [17] "Nativity_"
## [18] "Nativity_WaterContent 18-19"
## [19] "Plant_Height.html"
## [20] "Plant_Height.Rmd"
## [21] "Plant_Height.xlsx"
## [22] "Rplot.pdf"# Read the Excel file
data <- read_excel("/Users/ivyvelasco/Desktop/Greenhouse 2026/Plant_Height.xlsx")
# Display the raw data
head(data, 10)## # A tibble: 10 × 7
## Nativity Plot Species BioRep HeightT1 HeightT3 HeightChange
## <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Native 1 CEATHY 1 83 92.5 9.5
## 2 Native 1 CEATHY 3 122 135 13
## 3 Native 1 CEATHY 4 59 90 31
## 4 Native 1 CEATHY 5 95 57 -38
## 5 Native 3 CEATHY 6 78 76 -2
## 6 Native 3 CEATHY 7 66 76 10
## 7 Native 3 CEATHY 8 110 64 -46
## 8 Native 3 CEATHY 9 100 100 0
## 9 Native 3 CEATHY 10 61 55 -6
## 10 Native 1 ENCCAL 1 23 77 54cat("Total rows:", nrow(data), "\n")## Total rows: 198# Clean and prepare variables
data$Nativity <- trimws(data$Nativity)
data$Species <- trimws(data$Species)
data$Nativity <- factor(data$Nativity, levels = c("Non-Native", "Native"))
data$Species <- factor(data$Species)
data$Plot <- factor(data$Plot)
# Check the cleaned data
str(data)## tibble [198 × 7] (S3: tbl_df/tbl/data.frame)
## $ Nativity : Factor w/ 2 levels "Non-Native","Native": 2 2 2 2 2 2 2 2 2 2 ...
## $ Plot : Factor w/ 3 levels "1","2","3": 1 1 1 1 3 3 3 3 3 1 ...
## $ Species : Factor w/ 22 levels "CALSPP","CEATHY",..: 2 2 2 2 2 2 2 2 2 6 ...
## $ BioRep : num [1:198] 1 3 4 5 6 7 8 9 10 1 ...
## $ HeightT1 : num [1:198] 83 122 59 95 78 66 110 100 61 23 ...
## $ HeightT3 : num [1:198] 92.5 135 90 57 76 76 64 100 55 77 ...
## $ HeightChange: num [1:198] 9.5 13 31 -38 -2 10 -46 0 -6 54 ...Step 3: Verify Nesting Structure
# Check if species are truly nested within nativity
nesting_check <- data %>%
group_by(Species) %>%
summarise(
Nativities = n_distinct(Nativity),
Nativity_Type = paste(unique(Nativity), collapse = ", ")
)
print("Species Nesting Check:")## [1] "Species Nesting Check:"print(nesting_check)## # A tibble: 22 × 3
## Species Nativities Nativity_Type
## <fct> <int> <chr>
## 1 CALSPP 1 Non-Native
## 2 CEATHY 1 Native
## 3 CISxPUR 1 Non-Native
## 4 ELAPUN 1 Non-Native
## 5 ELYCON 1 Native
## 6 ENCCAL 1 Native
## 7 EPICAN 1 Native
## 8 ERIFAS 1 Native
## 9 FRACAL 1 Native
## 10 HETARB 1 Native
## # ℹ 12 more rowsspecies_in_multiple <- nesting_check %>% filter(Nativities > 1)
if(nrow(species_in_multiple) > 0) {
cat("\nWARNING: These species appear in multiple nativity levels:\n")
print(species_in_multiple)
cat("Species are NOT fully nested - the nesting structure may not be appropriate!\n")
} else {
cat("\n✓ CONFIRMED: All species are nested within a single nativity level\n")
cat(" Each species is EITHER Native OR Non-Native (never both)\n")
}##
## ✓ CONFIRMED: All species are nested within a single nativity level
## Each species is EITHER Native OR Non-Native (never both)# Show distribution
species_counts <- data %>%
group_by(Nativity) %>%
summarise(n_species = n_distinct(Species))
cat("\nSpecies per Nativity:\n")##
## Species per Nativity:print(species_counts)## # A tibble: 2 × 2
## Nativity n_species
## <fct> <int>
## 1 Non-Native 10
## 2 Native 12Step 4: Descriptive Statistics
# Overall by Nativity
summary_nativity <- data %>%
group_by(Nativity) %>%
summarise(
n = n(),
n_species = n_distinct(Species),
Mean = mean(HeightChange, na.rm = TRUE),
SD = sd(HeightChange, na.rm = TRUE),
SE = sd(HeightChange, na.rm = TRUE) / sqrt(n()),
Median = median(HeightChange, na.rm = TRUE),
Min = min(HeightChange, na.rm = TRUE),
Max = max(HeightChange, na.rm = TRUE)
)
cat("By Nativity:\n")## By Nativity:print(summary_nativity)## # A tibble: 2 × 9
## Nativity n n_species Mean SD SE Median Min Max
## <fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Non-Native 88 10 3.09 15.8 1.69 4.5 -58 30
## 2 Native 110 12 12.6 28.1 2.68 5 -46 131# By Species (nested within Nativity)
summary_species <- data %>%
group_by(Nativity, Species) %>%
summarise(
n = n(),
Mean = mean(HeightChange, na.rm = TRUE),
SD = sd(HeightChange, na.rm = TRUE),
.groups = 'drop'
) %>%
arrange(Nativity, desc(Mean))
cat("\nBy Species (showing species-level variation within each nativity):\n")##
## By Species (showing species-level variation within each nativity):print(summary_species)## # A tibble: 22 × 5
## Nativity Species n Mean SD
## <fct> <fct> <int> <dbl> <dbl>
## 1 Non-Native SAL GRExMIC 3 14 9.54
## 2 Non-Native LAVGIN 10 11.1 6.76
## 3 Non-Native ROSOFF 10 10.9 6.30
## 4 Non-Native SAL GRExMIC 6 10.1 13.3
## 5 Non-Native CISxPUR 10 7.58 8.01
## 6 Non-Native PERATR 10 5.2 7.93
## 7 Non-Native ELAPUN 9 0 4.33
## 8 Non-Native RHAIND 10 -0.68 9.69
## 9 Non-Native CALSPP 10 -3.6 4.67
## 10 Non-Native PENSPA 10 -13.5 36.5
## # ℹ 12 more rowsStep 5: Visualizations
Overall Boxplot by Nativity
p1 <- ggplot(data, aes(x = Nativity, y = HeightChange, fill = Nativity)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.4, size = 2) +
scale_fill_manual(values = c("Non-Native" = "#E74C3C", "Native" = "#27AE60")) +
labs(
title = "Plant Height Change by Nativity",
x = "Plant Nativity",
y = "Height Change (mm)",
fill = "Nativity"
) +
theme_minimal(base_size = 14) +
theme(legend.position = "none") +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50")
print(p1)
Species Means within Each Nativity
p2 <- ggplot(summary_species, aes(x = reorder(Species, Mean), y = Mean, fill = Nativity)) +
geom_bar(stat = "identity", alpha = 0.7) +
scale_fill_manual(values = c("Non-Native" = "#E74C3C", "Native" = "#27AE60")) +
coord_flip() +
facet_wrap(~ Nativity, scales = "free_y", ncol = 1) +
labs(
title = "Mean Height Change by Species (Nested within Nativity)",
subtitle = "Species are grouped by their nativity status",
x = "Species",
y = "Mean Height Change (mm)",
fill = "Nativity"
) +
theme_minimal(base_size = 11) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50")
print(p2)
Species-Level Responses
p3 <- ggplot(data, aes(x = Nativity, y = HeightChange, group = Species, color = Species)) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun = mean, geom = "line", aes(group = Species), alpha = 0.5) +
scale_color_viridis_d() +
labs(
title = "Species-Level Responses within Each Nativity",
subtitle = "Each line represents a species; shows variation within nativity groups",
x = "Plant Nativity",
y = "Mean Height Change (mm)"
) +
theme_minimal(base_size = 14) +
theme(legend.position = "right") +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50")
print(p3)
Step 6: Fit Mixed Model with Species as Fixed Effect
Model Formula:
HeightChange ~ Nativity/Species + (1|Plot)
This is equivalent to:
HeightChange ~ Nativity + Nativity:Species + (1|Plot)
What this means:
- Nativity: Main effect of being Native vs Non-Native
- Nativity:Species: Species-specific effects WITHIN each nativity level
- (1|Plot): Random effect accounting for spatial clustering
CRITICAL: Species is now FIXED, meaning we’re only making inferences about THESE specific species, not generalizing to other species.
model_nested <- lmer(HeightChange ~ Nativity/Species + (1|Plot),
data = data,
REML = TRUE)
summary(model_nested)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: HeightChange ~ Nativity/Species + (1 | Plot)
## Data: data
##
## REML criterion at convergence: 1580.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3827 -0.3854 -0.0036 0.3540 5.8232
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 29.99 5.477
## Residual 348.82 18.677
## Number of obs: 198, groups: Plot, 3
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -4.5462 6.7160 24.1022 -0.677
## NativityNative 4.2919 8.3740 174.4315 0.513
## NativityNative:SpeciesCEATHY -5.1495 8.7667 175.9993 -0.587
## NativityNon-Native:SpeciesCISxPUR 11.1281 8.3740 174.4315 1.329
## NativityNon-Native:SpeciesELAPUN 3.9725 8.6115 174.5359 0.461
## NativityNative:SpeciesELYCON -14.6685 8.8590 172.0202 -1.656
## NativityNative:SpeciesENCCAL 46.9009 8.8919 174.5159 5.275
## NativityNative:SpeciesEPICAN 7.7783 8.5544 175.9986 0.909
## NativityNative:SpeciesERIFAS 18.7988 8.6386 174.9062 2.176
## NativityNative:SpeciesFRACAL -1.3919 8.3740 174.4315 -0.166
## NativityNative:SpeciesHETARB 17.6792 8.5186 175.9964 2.075
## NativityNon-Native:SpeciesLAVGIN 14.7000 8.3525 173.9961 1.760
## NativityNative:SpeciesMYRCAL 0.2947 8.4842 175.9067 0.035
## NativityNon-Native:SpeciesPENSPA -9.1267 8.6789 175.3008 -1.052
## NativityNon-Native:SpeciesPERATR 8.8000 8.3525 173.9961 1.054
## NativityNon-Native:SpeciesRHAIND 2.8681 8.3740 174.4315 0.343
## NativityNative:SpeciesRHUINT 4.8480 11.1266 175.1616 0.436
## NativityNon-Native:SpeciesROSOFF 14.4081 8.3740 174.4315 1.721
## NativityNon-Native:SpeciesSAL GRExMIC 11.3603 9.6921 174.8520 1.172
## NativityNon-Native:SpeciesSAL GRExMIC 15.2770 12.3318 174.5521 1.239
## NativityNative:SpeciesSALAPI 44.6408 8.3595 174.1407 5.340
## NativityNative:SpeciesSALCLE 29.4370 8.4358 175.2739 3.490
## Pr(>|t|)
## (Intercept) 0.504900
## NativityNative 0.608932
## NativityNative:SpeciesCEATHY 0.557694
## NativityNon-Native:SpeciesCISxPUR 0.185617
## NativityNon-Native:SpeciesELAPUN 0.645154
## NativityNative:SpeciesELYCON 0.099591 .
## NativityNative:SpeciesENCCAL 3.90e-07 ***
## NativityNative:SpeciesEPICAN 0.364449
## NativityNative:SpeciesERIFAS 0.030886 *
## NativityNative:SpeciesFRACAL 0.868182
## NativityNative:SpeciesHETARB 0.039408 *
## NativityNon-Native:SpeciesLAVGIN 0.080172 .
## NativityNative:SpeciesMYRCAL 0.972328
## NativityNon-Native:SpeciesPENSPA 0.294433
## NativityNon-Native:SpeciesPERATR 0.293535
## NativityNon-Native:SpeciesRHAIND 0.732382
## NativityNative:SpeciesRHUINT 0.663579
## NativityNon-Native:SpeciesROSOFF 0.087099 .
## NativityNon-Native:SpeciesSAL GRExMIC 0.242743
## NativityNon-Native:SpeciesSAL GRExMIC 0.217072
## NativityNative:SpeciesSALAPI 2.87e-07 ***
## NativityNative:SpeciesSALCLE 0.000612 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 22 columns / coefficientsStep 7: ANOVA Table for Fixed Effects
anova_results <- anova(model_nested, type = "III")
print(anova_results)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Nativity 92 91.63 1 174.43 0.2627 0.6089
## Nativity:Species 41166 2058.31 20 172.87 5.9008 1.252e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\nInterpretation of ANOVA:\n")##
## Interpretation of ANOVA:cat("- Nativity: Tests if there's an overall difference between Native and Non-Native\n")## - Nativity: Tests if there's an overall difference between Native and Non-Nativecat("- Nativity:Species: Tests if species within each nativity differ from each other\n")## - Nativity:Species: Tests if species within each nativity differ from each othercat(" (i.e., is there significant species-level variation within nativity groups?)\n")## (i.e., is there significant species-level variation within nativity groups?)Step 8: Main Fixed Effects
# Extract fixed effects coefficients
fixed_effects <- summary(model_nested)$coefficients
cat("Fixed Effects Coefficients (first 10):\n")## Fixed Effects Coefficients (first 10):print(head(fixed_effects, 10))## Estimate Std. Error df t value
## (Intercept) -4.546198 6.715965 24.10221 -0.6769239
## NativityNative 4.291863 8.373961 174.43149 0.5125248
## NativityNative:SpeciesCEATHY -5.149453 8.766651 175.99927 -0.5873911
## NativityNon-Native:SpeciesCISxPUR 11.128137 8.373961 174.43149 1.3288977
## NativityNon-Native:SpeciesELAPUN 3.972511 8.611465 174.53593 0.4613049
## NativityNative:SpeciesELYCON -14.668542 8.859024 172.02023 -1.6557740
## NativityNative:SpeciesENCCAL 46.900944 8.891877 174.51590 5.2745832
## NativityNative:SpeciesEPICAN 7.778320 8.554434 175.99860 0.9092734
## NativityNative:SpeciesERIFAS 18.798779 8.638598 174.90615 2.1761377
## NativityNative:SpeciesFRACAL -1.391863 8.373961 174.43149 -0.1662132
## Pr(>|t|)
## (Intercept) 5.048999e-01
## NativityNative 6.089321e-01
## NativityNative:SpeciesCEATHY 5.576939e-01
## NativityNon-Native:SpeciesCISxPUR 1.856174e-01
## NativityNon-Native:SpeciesELAPUN 6.451542e-01
## NativityNative:SpeciesELYCON 9.959071e-02
## NativityNative:SpeciesENCCAL 3.903736e-07
## NativityNative:SpeciesEPICAN 3.644493e-01
## NativityNative:SpeciesERIFAS 3.088591e-02
## NativityNative:SpeciesFRACAL 8.681817e-01# Focus on the main Nativity effect
nativity_coef <- fixed_effects["NativityNative", ]
cat("\nMain Nativity Effect:\n")##
## Main Nativity Effect:cat("Estimate:", round(nativity_coef[1], 2), "mm\n")## Estimate: 4.29 mmcat("Std. Error:", round(nativity_coef[2], 2), "\n")## Std. Error: 8.37cat("t-value:", round(nativity_coef[4], 2), "\n")## t-value: 0.51cat("p-value:", round(nativity_coef[5], 4), "\n")## p-value: 0.6089if(nativity_coef[5] < 0.001) {
cat("→ HIGHLY SIGNIFICANT main effect of Nativity (p < 0.001)\n")
} else if(nativity_coef[5] < 0.01) {
cat("→ VERY SIGNIFICANT main effect of Nativity (p < 0.01)\n")
} else if(nativity_coef[5] < 0.05) {
cat("→ SIGNIFICANT main effect of Nativity (p < 0.05)\n")
} else {
cat("→ NO SIGNIFICANT main effect of Nativity (p ≥ 0.05)\n")
}## → NO SIGNIFICANT main effect of Nativity (p ≥ 0.05)# Random effects
cat("\nRandom Effects (Variance Components):\n")##
## Random Effects (Variance Components):random_effects <- as.data.frame(VarCorr(model_nested))
print(random_effects)## grp var1 var2 vcov sdcor
## 1 Plot (Intercept) <NA> 29.99269 5.476559
## 2 Residual <NA> <NA> 348.81774 18.676663plot_var <- random_effects$vcov[random_effects$grp == "Plot"]
residual_var <- random_effects$vcov[random_effects$grp == "Residual"]
cat("\nPlot variance:", round(plot_var, 2), "\n")##
## Plot variance: 29.99cat("Residual variance:", round(residual_var, 2), "\n")## Residual variance: 348.82Step 9: Estimated Marginal Means
# Get marginal means for Nativity
emm_nativity <- emmeans(model_nested, ~ Nativity)
cat("Estimated Marginal Means for Nativity:\n")## Estimated Marginal Means for Nativity:print(emm_nativity)## Nativity emmean SE df lower.CL upper.CL
## Non-Native 2.79 3.87 2.81 -10.01 15.6
## Native 12.18 3.67 2.33 -1.66 26.0
##
## Results are averaged over the levels of: Species
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Pairwise comparison
contrast_nativity <- pairs(emm_nativity)
cat("\nPairwise Comparison (Native vs Non-Native):\n")##
## Pairwise Comparison (Native vs Non-Native):print(contrast_nativity)## contrast estimate SE df t.ratio p.value
## (Non-Native) - Native -9.38 2.84 175 -3.300 0.0012
##
## Results are averaged over the levels of: Species
## Degrees-of-freedom method: kenward-roger# Get species-level means within each nativity
emm_species <- emmeans(model_nested, ~ Species | Nativity)
cat("\nEstimated Means by Species within Nativity:\n")##
## Estimated Means by Species within Nativity:print(summary(emm_species))## Nativity = Non-Native:
## Species emmean SE df lower.CL upper.CL
## CALSPP -4.5462 6.72 24.2 -18.42 9.324
## CISxPUR 6.5819 6.72 24.1 -7.28 20.448
## ELAPUN -0.5737 7.00 27.9 -14.92 13.770
## LAVGIN 10.1538 6.72 24.2 -3.72 24.024
## PENSPA -13.6729 7.14 24.4 -28.39 1.042
## PERATR 4.2538 6.72 24.2 -9.62 18.124
## RHAIND -1.6781 6.72 24.1 -15.54 12.188
## ROSOFF 9.8619 6.72 24.1 -4.00 23.728
## SAL GRExMIC 6.8141 8.37 48.2 -10.02 23.645
## SAL GRExMIC 10.7308 11.30 103.5 -11.73 33.187
##
## Nativity = Native:
## Species emmean SE df lower.CL upper.CL
## CEATHY -5.4038 7.12 27.8 -19.99 9.180
## ELYCON -14.9229 7.14 24.4 -29.64 -0.208
## ENCCAL 46.6466 7.33 33.0 31.74 61.556
## EPICAN 7.5240 7.07 24.1 -7.06 22.107
## ERIFAS 18.5444 6.98 28.0 4.24 32.848
## FRACAL -1.6462 6.72 24.2 -15.52 12.224
## HETARB 17.4248 6.83 24.1 3.33 31.520
## MYRCAL 0.0404 6.83 24.2 -14.05 14.134
## RHUINT 4.5937 9.98 78.3 -15.27 24.455
## SALAPI 44.3864 6.81 24.2 30.34 58.437
## SALCLE 29.1827 6.70 24.1 15.35 43.014
## SIMCHI -0.2543 6.81 24.3 -14.29 13.783
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Step 10: Model Diagnostics
# Extract residuals and fitted values
data$residuals <- residuals(model_nested)
data$fitted <- fitted(model_nested)
# Residuals vs Fitted
plot(data$fitted, data$residuals,
xlab = "Fitted Values", ylab = "Residuals",
main = "Residuals vs Fitted Values")
abline(h = 0, col = "red", lty = 2)
# Q-Q plot
qqnorm(data$residuals, main = "Q-Q Plot of Residuals")
qqline(data$residuals, col = "red")
# Residuals by Nativity
boxplot(residuals ~ Nativity, data = data,
main = "Residuals by Nativity",
ylab = "Residuals", xlab = "Nativity")
abline(h = 0, col = "red", lty = 2)
# Shapiro-Wilk test
if(nrow(data) > 5000) {
shapiro_sample <- sample(data$residuals, 5000)
shapiro_test <- shapiro.test(shapiro_sample)
} else {
shapiro_test <- shapiro.test(data$residuals)
}
cat("\nShapiro-Wilk Test for Normality:\n")##
## Shapiro-Wilk Test for Normality:print(shapiro_test)##
## Shapiro-Wilk normality test
##
## data: data$residuals
## W = 0.87102, p-value = 5.996e-12Step 11: Model Comparison
# Compare with the previous model (Species as random effect)
model_random <- lmer(HeightChange ~ Nativity + (1|Species) + (1|Plot),
data = data,
REML = FALSE)
model_nested_ML <- lmer(HeightChange ~ Nativity/Species + (1|Plot),
data = data,
REML = FALSE)
cat("Comparing Models:\n")## Comparing Models:cat("Model 1: Species as random effect\n")## Model 1: Species as random effectcat("Model 2: Species as fixed effect (nested)\n\n")## Model 2: Species as fixed effect (nested)comparison <- data.frame(
Model = c("Species Random", "Species Fixed (Nested)"),
AIC = c(AIC(model_random), AIC(model_nested_ML)),
BIC = c(BIC(model_random), BIC(model_nested_ML)),
LogLik = c(logLik(model_random), logLik(model_nested_ML))
)
print(comparison)## Model AIC BIC LogLik
## 1 Species Random 1772.043 1788.484 -881.0213
## 2 Species Fixed (Nested) 1751.246 1830.165 -851.6232cat("\nInterpretation:\n")##
## Interpretation:cat("- Lower AIC/BIC indicates better model fit\n")## - Lower AIC/BIC indicates better model fitcat("- The nested (fixed) model will typically have lower AIC/BIC\n")## - The nested (fixed) model will typically have lower AIC/BICcat(" because it uses more parameters (one per species)\n")## because it uses more parameters (one per species)cat("- But this comes at a cost: less generalizability\n")## - But this comes at a cost: less generalizabilityStep 12: Post-hoc Comparisons
cat("Comparing species within Native group:\n")## Comparing species within Native group:emm_native_species <- emmeans(model_nested, ~ Species | Nativity,
at = list(Nativity = "Native"))
pairs_native <- pairs(emm_native_species)
print(summary(pairs_native))## Nativity = Native:
## contrast estimate SE df t.ratio p.value
## CALSPP - CEATHY nonEst NA NA NA NA
## CALSPP - CISxPUR nonEst NA NA NA NA
## CALSPP - ELAPUN nonEst NA NA NA NA
## CALSPP - ELYCON nonEst NA NA NA NA
## CALSPP - ENCCAL nonEst NA NA NA NA
## CALSPP - EPICAN nonEst NA NA NA NA
## CALSPP - ERIFAS nonEst NA NA NA NA
## CALSPP - FRACAL nonEst NA NA NA NA
## CALSPP - HETARB nonEst NA NA NA NA
## CALSPP - LAVGIN nonEst NA NA NA NA
## CALSPP - MYRCAL nonEst NA NA NA NA
## CALSPP - PENSPA nonEst NA NA NA NA
## CALSPP - PERATR nonEst NA NA NA NA
## CALSPP - RHAIND nonEst NA NA NA NA
## CALSPP - RHUINT nonEst NA NA NA NA
## CALSPP - ROSOFF nonEst NA NA NA NA
## CALSPP - SAL GRExMIC nonEst NA NA NA NA
## CALSPP - SAL GRExMIC nonEst NA NA NA NA
## CALSPP - SALAPI nonEst NA NA NA NA
## CALSPP - SALCLE nonEst NA NA NA NA
## CALSPP - SIMCHI nonEst NA NA NA NA
## CEATHY - CISxPUR nonEst NA NA NA NA
## CEATHY - ELAPUN nonEst NA NA NA NA
## CEATHY - ELYCON 9.519 8.74 175 1.089 0.9948
## CEATHY - ENCCAL -52.050 9.22 176 -5.646 <0.0001
## CEATHY - EPICAN -12.928 9.31 171 -1.389 0.9642
## CEATHY - ERIFAS -23.948 8.91 175 -2.687 0.2415
## CEATHY - FRACAL -3.758 8.71 175 -0.431 1.0000
## CEATHY - HETARB -22.829 8.58 174 -2.659 0.2558
## CEATHY - LAVGIN nonEst NA NA NA NA
## CEATHY - MYRCAL -5.444 8.60 174 -0.633 1.0000
## CEATHY - PENSPA nonEst NA NA NA NA
## CEATHY - PERATR nonEst NA NA NA NA
## CEATHY - RHAIND nonEst NA NA NA NA
## CEATHY - RHUINT -9.998 11.30 174 -0.886 0.9992
## CEATHY - ROSOFF nonEst NA NA NA NA
## CEATHY - SAL GRExMIC nonEst NA NA NA NA
## CEATHY - SAL GRExMIC nonEst NA NA NA NA
## CEATHY - SALAPI -49.790 8.92 176 -5.581 <0.0001
## CEATHY - SALCLE -34.587 8.67 175 -3.991 0.0053
## CEATHY - SIMCHI -5.149 8.86 176 -0.581 1.0000
## CISxPUR - ELAPUN nonEst NA NA NA NA
## CISxPUR - ELYCON nonEst NA NA NA NA
## CISxPUR - ENCCAL nonEst NA NA NA NA
## CISxPUR - EPICAN nonEst NA NA NA NA
## CISxPUR - ERIFAS nonEst NA NA NA NA
## CISxPUR - FRACAL nonEst NA NA NA NA
## CISxPUR - HETARB nonEst NA NA NA NA
## CISxPUR - LAVGIN nonEst NA NA NA NA
## CISxPUR - MYRCAL nonEst NA NA NA NA
## CISxPUR - PENSPA nonEst NA NA NA NA
## CISxPUR - PERATR nonEst NA NA NA NA
## CISxPUR - RHAIND nonEst NA NA NA NA
## CISxPUR - RHUINT nonEst NA NA NA NA
## CISxPUR - ROSOFF nonEst NA NA NA NA
## CISxPUR - SAL GRExMIC nonEst NA NA NA NA
## CISxPUR - SAL GRExMIC nonEst NA NA NA NA
## CISxPUR - SALAPI nonEst NA NA NA NA
## CISxPUR - SALCLE nonEst NA NA NA NA
## CISxPUR - SIMCHI nonEst NA NA NA NA
## ELAPUN - ELYCON nonEst NA NA NA NA
## ELAPUN - ENCCAL nonEst NA NA NA NA
## ELAPUN - EPICAN nonEst NA NA NA NA
## ELAPUN - ERIFAS nonEst NA NA NA NA
## ELAPUN - FRACAL nonEst NA NA NA NA
## ELAPUN - HETARB nonEst NA NA NA NA
## ELAPUN - LAVGIN nonEst NA NA NA NA
## ELAPUN - MYRCAL nonEst NA NA NA NA
## ELAPUN - PENSPA nonEst NA NA NA NA
## ELAPUN - PERATR nonEst NA NA NA NA
## ELAPUN - RHAIND nonEst NA NA NA NA
## ELAPUN - RHUINT nonEst NA NA NA NA
## ELAPUN - ROSOFF nonEst NA NA NA NA
## ELAPUN - SAL GRExMIC nonEst NA NA NA NA
## ELAPUN - SAL GRExMIC nonEst NA NA NA NA
## ELAPUN - SALAPI nonEst NA NA NA NA
## ELAPUN - SALCLE nonEst NA NA NA NA
## ELAPUN - SIMCHI nonEst NA NA NA NA
## ELYCON - ENCCAL -61.569 9.28 176 -6.634 <0.0001
## ELYCON - EPICAN -22.447 9.34 164 -2.403 0.4102
## ELYCON - ERIFAS -33.467 8.93 176 -3.749 0.0125
## ELYCON - FRACAL -13.277 8.83 175 -1.503 0.9381
## ELYCON - HETARB -32.348 8.56 176 -3.777 0.0114
## ELYCON - LAVGIN nonEst NA NA NA NA
## ELYCON - MYRCAL -14.963 8.66 176 -1.729 0.8525
## ELYCON - PENSPA nonEst NA NA NA NA
## ELYCON - PERATR nonEst NA NA NA NA
## ELYCON - RHAIND nonEst NA NA NA NA
## ELYCON - RHUINT -19.517 11.40 176 -1.711 0.8610
## ELYCON - ROSOFF nonEst NA NA NA NA
## ELYCON - SAL GRExMIC nonEst NA NA NA NA
## ELYCON - SAL GRExMIC nonEst NA NA NA NA
## ELYCON - SALAPI -59.309 9.13 171 -6.497 <0.0001
## ELYCON - SALCLE -44.106 8.64 176 -5.105 <0.0001
## ELYCON - SIMCHI -14.669 9.10 172 -1.612 0.9025
## ENCCAL - EPICAN 39.123 9.10 176 4.298 0.0017
## ENCCAL - ERIFAS 28.102 9.08 174 3.095 0.0919
## ENCCAL - FRACAL 48.293 8.87 174 5.445 <0.0001
## ENCCAL - HETARB 29.222 8.99 175 3.249 0.0601
## ENCCAL - LAVGIN nonEst NA NA NA NA
## ENCCAL - MYRCAL 46.606 8.98 175 5.188 <0.0001
## ENCCAL - PENSPA nonEst NA NA NA NA
## ENCCAL - PERATR nonEst NA NA NA NA
## ENCCAL - RHAIND nonEst NA NA NA NA
## ENCCAL - RHUINT 42.053 11.50 175 3.643 0.0179
## ENCCAL - ROSOFF nonEst NA NA NA NA
## ENCCAL - SAL GRExMIC nonEst NA NA NA NA
## ENCCAL - SAL GRExMIC nonEst NA NA NA NA
## ENCCAL - SALAPI 2.260 8.91 175 0.254 1.0000
## ENCCAL - SALCLE 17.464 8.88 174 1.967 0.7145
## ENCCAL - SIMCHI 46.901 8.91 175 5.266 <0.0001
## EPICAN - ERIFAS -11.020 8.87 176 -1.242 0.9847
## EPICAN - FRACAL 9.170 8.70 176 1.055 0.9961
## EPICAN - HETARB -9.901 9.09 171 -1.090 0.9948
## EPICAN - LAVGIN nonEst NA NA NA NA
## EPICAN - MYRCAL 7.484 9.08 171 0.824 0.9996
## EPICAN - PENSPA nonEst NA NA NA NA
## EPICAN - PERATR nonEst NA NA NA NA
## EPICAN - RHAIND nonEst NA NA NA NA
## EPICAN - RHUINT 2.930 11.60 175 0.252 1.0000
## EPICAN - ROSOFF nonEst NA NA NA NA
## EPICAN - SAL GRExMIC nonEst NA NA NA NA
## EPICAN - SAL GRExMIC nonEst NA NA NA NA
## EPICAN - SALAPI -36.862 8.55 176 -4.311 0.0016
## EPICAN - SALCLE -21.659 8.69 176 -2.492 0.3518
## EPICAN - SIMCHI 7.778 8.65 176 0.899 0.9991
## ERIFAS - FRACAL 20.191 8.60 174 2.348 0.4473
## ERIFAS - HETARB 1.120 8.68 175 0.129 1.0000
## ERIFAS - LAVGIN nonEst NA NA NA NA
## ERIFAS - MYRCAL 18.504 8.68 175 2.131 0.6017
## ERIFAS - PENSPA nonEst NA NA NA NA
## ERIFAS - PERATR nonEst NA NA NA NA
## ERIFAS - RHAIND nonEst NA NA NA NA
## ERIFAS - RHUINT 13.951 11.30 175 1.232 0.9857
## ERIFAS - ROSOFF nonEst NA NA NA NA
## ERIFAS - SAL GRExMIC nonEst NA NA NA NA
## ERIFAS - SAL GRExMIC nonEst NA NA NA NA
## ERIFAS - SALAPI -25.842 8.67 175 -2.981 0.1232
## ERIFAS - SALCLE -10.638 8.58 174 -1.239 0.9850
## ERIFAS - SIMCHI 18.799 8.66 175 2.170 0.5738
## FRACAL - HETARB -19.071 8.47 175 -2.252 0.5146
## FRACAL - LAVGIN nonEst NA NA NA NA
## FRACAL - MYRCAL -1.687 8.44 175 -0.200 1.0000
## FRACAL - PENSPA nonEst NA NA NA NA
## FRACAL - PERATR nonEst NA NA NA NA
## FRACAL - RHAIND nonEst NA NA NA NA
## FRACAL - RHUINT -6.240 11.10 175 -0.561 1.0000
## FRACAL - ROSOFF nonEst NA NA NA NA
## FRACAL - SAL GRExMIC nonEst NA NA NA NA
## FRACAL - SAL GRExMIC nonEst NA NA NA NA
## FRACAL - SALAPI -46.033 8.40 175 -5.477 <0.0001
## FRACAL - SALCLE -30.829 8.39 174 -3.676 0.0160
## FRACAL - SIMCHI -1.392 8.38 174 -0.166 1.0000
## HETARB - LAVGIN nonEst NA NA NA NA
## HETARB - MYRCAL 17.384 8.36 174 2.079 0.6381
## HETARB - PENSPA nonEst NA NA NA NA
## HETARB - PERATR nonEst NA NA NA NA
## HETARB - RHAIND nonEst NA NA NA NA
## HETARB - RHUINT 12.831 11.10 174 1.157 0.9914
## HETARB - ROSOFF nonEst NA NA NA NA
## HETARB - SAL GRExMIC nonEst NA NA NA NA
## HETARB - SAL GRExMIC nonEst NA NA NA NA
## HETARB - SALAPI -26.962 8.67 176 -3.110 0.0882
## HETARB - SALCLE -11.758 8.44 175 -1.394 0.9634
## HETARB - SIMCHI 17.679 8.61 176 2.054 0.6559
## LAVGIN - MYRCAL nonEst NA NA NA NA
## LAVGIN - PENSPA nonEst NA NA NA NA
## LAVGIN - PERATR nonEst NA NA NA NA
## LAVGIN - RHAIND nonEst NA NA NA NA
## LAVGIN - RHUINT nonEst NA NA NA NA
## LAVGIN - ROSOFF nonEst NA NA NA NA
## LAVGIN - SAL GRExMIC nonEst NA NA NA NA
## LAVGIN - SAL GRExMIC nonEst NA NA NA NA
## LAVGIN - SALAPI nonEst NA NA NA NA
## LAVGIN - SALCLE nonEst NA NA NA NA
## LAVGIN - SIMCHI nonEst NA NA NA NA
## MYRCAL - PENSPA nonEst NA NA NA NA
## MYRCAL - PERATR nonEst NA NA NA NA
## MYRCAL - RHAIND nonEst NA NA NA NA
## MYRCAL - RHUINT -4.553 11.10 174 -0.412 1.0000
## MYRCAL - ROSOFF nonEst NA NA NA NA
## MYRCAL - SAL GRExMIC nonEst NA NA NA NA
## MYRCAL - SAL GRExMIC nonEst NA NA NA NA
## MYRCAL - SALAPI -44.346 8.63 176 -5.142 <0.0001
## MYRCAL - SALCLE -29.142 8.45 175 -3.450 0.0331
## MYRCAL - SIMCHI 0.295 8.56 176 0.034 1.0000
## PENSPA - PERATR nonEst NA NA NA NA
## PENSPA - RHAIND nonEst NA NA NA NA
## PENSPA - RHUINT nonEst NA NA NA NA
## PENSPA - ROSOFF nonEst NA NA NA NA
## PENSPA - SAL GRExMIC nonEst NA NA NA NA
## PENSPA - SAL GRExMIC nonEst NA NA NA NA
## PENSPA - SALAPI nonEst NA NA NA NA
## PENSPA - SALCLE nonEst NA NA NA NA
## PENSPA - SIMCHI nonEst NA NA NA NA
## PERATR - RHAIND nonEst NA NA NA NA
## PERATR - RHUINT nonEst NA NA NA NA
## PERATR - ROSOFF nonEst NA NA NA NA
## PERATR - SAL GRExMIC nonEst NA NA NA NA
## PERATR - SAL GRExMIC nonEst NA NA NA NA
## PERATR - SALAPI nonEst NA NA NA NA
## PERATR - SALCLE nonEst NA NA NA NA
## PERATR - SIMCHI nonEst NA NA NA NA
## RHAIND - RHUINT nonEst NA NA NA NA
## RHAIND - ROSOFF nonEst NA NA NA NA
## RHAIND - SAL GRExMIC nonEst NA NA NA NA
## RHAIND - SAL GRExMIC nonEst NA NA NA NA
## RHAIND - SALAPI nonEst NA NA NA NA
## RHAIND - SALCLE nonEst NA NA NA NA
## RHAIND - SIMCHI nonEst NA NA NA NA
## RHUINT - ROSOFF nonEst NA NA NA NA
## RHUINT - SAL GRExMIC nonEst NA NA NA NA
## RHUINT - SAL GRExMIC nonEst NA NA NA NA
## RHUINT - SALAPI -39.793 11.20 176 -3.544 0.0246
## RHUINT - SALCLE -24.589 11.20 175 -2.204 0.5490
## RHUINT - SIMCHI 4.848 11.20 175 0.434 1.0000
## ROSOFF - SAL GRExMIC nonEst NA NA NA NA
## ROSOFF - SAL GRExMIC nonEst NA NA NA NA
## ROSOFF - SALAPI nonEst NA NA NA NA
## ROSOFF - SALCLE nonEst NA NA NA NA
## ROSOFF - SIMCHI nonEst NA NA NA NA
## SAL GRExMIC - SAL GRExMIC nonEst NA NA NA NA
## SAL GRExMIC - SALAPI nonEst NA NA NA NA
## SAL GRExMIC - SALCLE nonEst NA NA NA NA
## SAL GRExMIC - SIMCHI nonEst NA NA NA NA
## SAL GRExMIC - SALAPI nonEst NA NA NA NA
## SAL GRExMIC - SALCLE nonEst NA NA NA NA
## SAL GRExMIC - SIMCHI nonEst NA NA NA NA
## SALAPI - SALCLE 15.204 8.48 175 1.793 0.8199
## SALAPI - SIMCHI 44.641 8.36 174 5.338 <0.0001
## SALCLE - SIMCHI 29.437 8.47 175 3.474 0.0308
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 12 estimatescat("\n\nComparing species within Non-Native group:\n")##
##
## Comparing species within Non-Native group:emm_nonnative_species <- emmeans(model_nested, ~ Species | Nativity,
at = list(Nativity = "Non-Native"))
pairs_nonnative <- pairs(emm_nonnative_species)
print(summary(pairs_nonnative))## Nativity = Non-Native:
## contrast estimate SE df t.ratio p.value
## CALSPP - CEATHY nonEst NA NA NA NA
## CALSPP - CISxPUR -11.128 8.38 174 -1.327 0.9461
## CALSPP - ELAPUN -3.973 8.63 175 -0.461 1.0000
## CALSPP - ELYCON nonEst NA NA NA NA
## CALSPP - ENCCAL nonEst NA NA NA NA
## CALSPP - EPICAN nonEst NA NA NA NA
## CALSPP - ERIFAS nonEst NA NA NA NA
## CALSPP - FRACAL nonEst NA NA NA NA
## CALSPP - HETARB nonEst NA NA NA NA
## CALSPP - LAVGIN -14.700 8.35 174 -1.760 0.7591
## CALSPP - MYRCAL nonEst NA NA NA NA
## CALSPP - PENSPA 9.127 8.83 175 1.033 0.9898
## CALSPP - PERATR -8.800 8.35 174 -1.054 0.9883
## CALSPP - RHAIND -2.868 8.38 174 -0.342 1.0000
## CALSPP - RHUINT nonEst NA NA NA NA
## CALSPP - ROSOFF -14.408 8.38 174 -1.718 0.7837
## CALSPP - SAL GRExMIC -11.360 9.72 175 -1.169 0.9761
## CALSPP - SAL GRExMIC -15.277 12.40 175 -1.237 0.9653
## CALSPP - SALAPI nonEst NA NA NA NA
## CALSPP - SALCLE nonEst NA NA NA NA
## CALSPP - SIMCHI nonEst NA NA NA NA
## CEATHY - CISxPUR nonEst NA NA NA NA
## CEATHY - ELAPUN nonEst NA NA NA NA
## CEATHY - ELYCON nonEst NA NA NA NA
## CEATHY - ENCCAL nonEst NA NA NA NA
## CEATHY - EPICAN nonEst NA NA NA NA
## CEATHY - ERIFAS nonEst NA NA NA NA
## CEATHY - FRACAL nonEst NA NA NA NA
## CEATHY - HETARB nonEst NA NA NA NA
## CEATHY - LAVGIN nonEst NA NA NA NA
## CEATHY - MYRCAL nonEst NA NA NA NA
## CEATHY - PENSPA nonEst NA NA NA NA
## CEATHY - PERATR nonEst NA NA NA NA
## CEATHY - RHAIND nonEst NA NA NA NA
## CEATHY - RHUINT nonEst NA NA NA NA
## CEATHY - ROSOFF nonEst NA NA NA NA
## CEATHY - SAL GRExMIC nonEst NA NA NA NA
## CEATHY - SAL GRExMIC nonEst NA NA NA NA
## CEATHY - SALAPI nonEst NA NA NA NA
## CEATHY - SALCLE nonEst NA NA NA NA
## CEATHY - SIMCHI nonEst NA NA NA NA
## CISxPUR - ELAPUN 7.156 8.58 174 0.834 0.9979
## CISxPUR - ELYCON nonEst NA NA NA NA
## CISxPUR - ENCCAL nonEst NA NA NA NA
## CISxPUR - EPICAN nonEst NA NA NA NA
## CISxPUR - ERIFAS nonEst NA NA NA NA
## CISxPUR - FRACAL nonEst NA NA NA NA
## CISxPUR - HETARB nonEst NA NA NA NA
## CISxPUR - LAVGIN -3.572 8.38 174 -0.426 1.0000
## CISxPUR - MYRCAL nonEst NA NA NA NA
## CISxPUR - PENSPA 20.255 8.63 176 2.348 0.3637
## CISxPUR - PERATR 2.328 8.38 174 0.278 1.0000
## CISxPUR - RHAIND 8.260 8.35 174 0.989 0.9926
## CISxPUR - RHUINT nonEst NA NA NA NA
## CISxPUR - ROSOFF -3.280 8.35 174 -0.393 1.0000
## CISxPUR - SAL GRExMIC -0.232 9.69 174 -0.024 1.0000
## CISxPUR - SAL GRExMIC -4.149 12.30 174 -0.336 1.0000
## CISxPUR - SALAPI nonEst NA NA NA NA
## CISxPUR - SALCLE nonEst NA NA NA NA
## CISxPUR - SIMCHI nonEst NA NA NA NA
## ELAPUN - ELYCON nonEst NA NA NA NA
## ELAPUN - ENCCAL nonEst NA NA NA NA
## ELAPUN - EPICAN nonEst NA NA NA NA
## ELAPUN - ERIFAS nonEst NA NA NA NA
## ELAPUN - FRACAL nonEst NA NA NA NA
## ELAPUN - HETARB nonEst NA NA NA NA
## ELAPUN - LAVGIN -10.727 8.63 175 -1.244 0.9641
## ELAPUN - MYRCAL nonEst NA NA NA NA
## ELAPUN - PENSPA 13.099 8.81 176 1.486 0.8959
## ELAPUN - PERATR -4.827 8.63 175 -0.560 0.9999
## ELAPUN - RHAIND 1.104 8.58 174 0.129 1.0000
## ELAPUN - RHUINT nonEst NA NA NA NA
## ELAPUN - ROSOFF -10.436 8.58 174 -1.216 0.9690
## ELAPUN - SAL GRExMIC -7.388 9.91 175 -0.746 0.9991
## ELAPUN - SAL GRExMIC -11.304 12.50 174 -0.904 0.9962
## ELAPUN - SALAPI nonEst NA NA NA NA
## ELAPUN - SALCLE nonEst NA NA NA NA
## ELAPUN - SIMCHI nonEst NA NA NA NA
## ELYCON - ENCCAL nonEst NA NA NA NA
## ELYCON - EPICAN nonEst NA NA NA NA
## ELYCON - ERIFAS nonEst NA NA NA NA
## ELYCON - FRACAL nonEst NA NA NA NA
## ELYCON - HETARB nonEst NA NA NA NA
## ELYCON - LAVGIN nonEst NA NA NA NA
## ELYCON - MYRCAL nonEst NA NA NA NA
## ELYCON - PENSPA nonEst NA NA NA NA
## ELYCON - PERATR nonEst NA NA NA NA
## ELYCON - RHAIND nonEst NA NA NA NA
## ELYCON - RHUINT nonEst NA NA NA NA
## ELYCON - ROSOFF nonEst NA NA NA NA
## ELYCON - SAL GRExMIC nonEst NA NA NA NA
## ELYCON - SAL GRExMIC nonEst NA NA NA NA
## ELYCON - SALAPI nonEst NA NA NA NA
## ELYCON - SALCLE nonEst NA NA NA NA
## ELYCON - SIMCHI nonEst NA NA NA NA
## ENCCAL - EPICAN nonEst NA NA NA NA
## ENCCAL - ERIFAS nonEst NA NA NA NA
## ENCCAL - FRACAL nonEst NA NA NA NA
## ENCCAL - HETARB nonEst NA NA NA NA
## ENCCAL - LAVGIN nonEst NA NA NA NA
## ENCCAL - MYRCAL nonEst NA NA NA NA
## ENCCAL - PENSPA nonEst NA NA NA NA
## ENCCAL - PERATR nonEst NA NA NA NA
## ENCCAL - RHAIND nonEst NA NA NA NA
## ENCCAL - RHUINT nonEst NA NA NA NA
## ENCCAL - ROSOFF nonEst NA NA NA NA
## ENCCAL - SAL GRExMIC nonEst NA NA NA NA
## ENCCAL - SAL GRExMIC nonEst NA NA NA NA
## ENCCAL - SALAPI nonEst NA NA NA NA
## ENCCAL - SALCLE nonEst NA NA NA NA
## ENCCAL - SIMCHI nonEst NA NA NA NA
## EPICAN - ERIFAS nonEst NA NA NA NA
## EPICAN - FRACAL nonEst NA NA NA NA
## EPICAN - HETARB nonEst NA NA NA NA
## EPICAN - LAVGIN nonEst NA NA NA NA
## EPICAN - MYRCAL nonEst NA NA NA NA
## EPICAN - PENSPA nonEst NA NA NA NA
## EPICAN - PERATR nonEst NA NA NA NA
## EPICAN - RHAIND nonEst NA NA NA NA
## EPICAN - RHUINT nonEst NA NA NA NA
## EPICAN - ROSOFF nonEst NA NA NA NA
## EPICAN - SAL GRExMIC nonEst NA NA NA NA
## EPICAN - SAL GRExMIC nonEst NA NA NA NA
## EPICAN - SALAPI nonEst NA NA NA NA
## EPICAN - SALCLE nonEst NA NA NA NA
## EPICAN - SIMCHI nonEst NA NA NA NA
## ERIFAS - FRACAL nonEst NA NA NA NA
## ERIFAS - HETARB nonEst NA NA NA NA
## ERIFAS - LAVGIN nonEst NA NA NA NA
## ERIFAS - MYRCAL nonEst NA NA NA NA
## ERIFAS - PENSPA nonEst NA NA NA NA
## ERIFAS - PERATR nonEst NA NA NA NA
## ERIFAS - RHAIND nonEst NA NA NA NA
## ERIFAS - RHUINT nonEst NA NA NA NA
## ERIFAS - ROSOFF nonEst NA NA NA NA
## ERIFAS - SAL GRExMIC nonEst NA NA NA NA
## ERIFAS - SAL GRExMIC nonEst NA NA NA NA
## ERIFAS - SALAPI nonEst NA NA NA NA
## ERIFAS - SALCLE nonEst NA NA NA NA
## ERIFAS - SIMCHI nonEst NA NA NA NA
## FRACAL - HETARB nonEst NA NA NA NA
## FRACAL - LAVGIN nonEst NA NA NA NA
## FRACAL - MYRCAL nonEst NA NA NA NA
## FRACAL - PENSPA nonEst NA NA NA NA
## FRACAL - PERATR nonEst NA NA NA NA
## FRACAL - RHAIND nonEst NA NA NA NA
## FRACAL - RHUINT nonEst NA NA NA NA
## FRACAL - ROSOFF nonEst NA NA NA NA
## FRACAL - SAL GRExMIC nonEst NA NA NA NA
## FRACAL - SAL GRExMIC nonEst NA NA NA NA
## FRACAL - SALAPI nonEst NA NA NA NA
## FRACAL - SALCLE nonEst NA NA NA NA
## FRACAL - SIMCHI nonEst NA NA NA NA
## HETARB - LAVGIN nonEst NA NA NA NA
## HETARB - MYRCAL nonEst NA NA NA NA
## HETARB - PENSPA nonEst NA NA NA NA
## HETARB - PERATR nonEst NA NA NA NA
## HETARB - RHAIND nonEst NA NA NA NA
## HETARB - RHUINT nonEst NA NA NA NA
## HETARB - ROSOFF nonEst NA NA NA NA
## HETARB - SAL GRExMIC nonEst NA NA NA NA
## HETARB - SAL GRExMIC nonEst NA NA NA NA
## HETARB - SALAPI nonEst NA NA NA NA
## HETARB - SALCLE nonEst NA NA NA NA
## HETARB - SIMCHI nonEst NA NA NA NA
## LAVGIN - MYRCAL nonEst NA NA NA NA
## LAVGIN - PENSPA 23.827 8.83 175 2.697 0.1831
## LAVGIN - PERATR 5.900 8.35 174 0.706 0.9994
## LAVGIN - RHAIND 11.832 8.38 174 1.411 0.9224
## LAVGIN - RHUINT nonEst NA NA NA NA
## LAVGIN - ROSOFF 0.292 8.38 174 0.035 1.0000
## LAVGIN - SAL GRExMIC 3.340 9.72 175 0.344 1.0000
## LAVGIN - SAL GRExMIC -0.577 12.40 175 -0.047 1.0000
## LAVGIN - SALAPI nonEst NA NA NA NA
## LAVGIN - SALCLE nonEst NA NA NA NA
## LAVGIN - SIMCHI nonEst NA NA NA NA
## MYRCAL - PENSPA nonEst NA NA NA NA
## MYRCAL - PERATR nonEst NA NA NA NA
## MYRCAL - RHAIND nonEst NA NA NA NA
## MYRCAL - RHUINT nonEst NA NA NA NA
## MYRCAL - ROSOFF nonEst NA NA NA NA
## MYRCAL - SAL GRExMIC nonEst NA NA NA NA
## MYRCAL - SAL GRExMIC nonEst NA NA NA NA
## MYRCAL - SALAPI nonEst NA NA NA NA
## MYRCAL - SALCLE nonEst NA NA NA NA
## MYRCAL - SIMCHI nonEst NA NA NA NA
## PENSPA - PERATR -17.927 8.83 175 -2.029 0.5798
## PENSPA - RHAIND -11.995 8.63 176 -1.391 0.9287
## PENSPA - RHUINT nonEst NA NA NA NA
## PENSPA - ROSOFF -23.535 8.63 176 -2.729 0.1705
## PENSPA - SAL GRExMIC -20.487 9.97 176 -2.055 0.5616
## PENSPA - SAL GRExMIC -24.404 12.60 176 -1.944 0.6388
## PENSPA - SALAPI nonEst NA NA NA NA
## PENSPA - SALCLE nonEst NA NA NA NA
## PENSPA - SIMCHI nonEst NA NA NA NA
## PERATR - RHAIND 5.932 8.38 174 0.707 0.9994
## PERATR - RHUINT nonEst NA NA NA NA
## PERATR - ROSOFF -5.608 8.38 174 -0.669 0.9996
## PERATR - SAL GRExMIC -2.560 9.72 175 -0.263 1.0000
## PERATR - SAL GRExMIC -6.477 12.40 175 -0.524 1.0000
## PERATR - SALAPI nonEst NA NA NA NA
## PERATR - SALCLE nonEst NA NA NA NA
## PERATR - SIMCHI nonEst NA NA NA NA
## RHAIND - RHUINT nonEst NA NA NA NA
## RHAIND - ROSOFF -11.540 8.35 174 -1.382 0.9313
## RHAIND - SAL GRExMIC -8.492 9.69 174 -0.876 0.9970
## RHAIND - SAL GRExMIC -12.409 12.30 174 -1.006 0.9916
## RHAIND - SALAPI nonEst NA NA NA NA
## RHAIND - SALCLE nonEst NA NA NA NA
## RHAIND - SIMCHI nonEst NA NA NA NA
## RHUINT - ROSOFF nonEst NA NA NA NA
## RHUINT - SAL GRExMIC nonEst NA NA NA NA
## RHUINT - SAL GRExMIC nonEst NA NA NA NA
## RHUINT - SALAPI nonEst NA NA NA NA
## RHUINT - SALCLE nonEst NA NA NA NA
## RHUINT - SIMCHI nonEst NA NA NA NA
## ROSOFF - SAL GRExMIC 3.048 9.69 174 0.314 1.0000
## ROSOFF - SAL GRExMIC -0.869 12.30 174 -0.070 1.0000
## ROSOFF - SALAPI nonEst NA NA NA NA
## ROSOFF - SALCLE nonEst NA NA NA NA
## ROSOFF - SIMCHI nonEst NA NA NA NA
## SAL GRExMIC - SAL GRExMIC -3.917 13.20 174 -0.297 1.0000
## SAL GRExMIC - SALAPI nonEst NA NA NA NA
## SAL GRExMIC - SALCLE nonEst NA NA NA NA
## SAL GRExMIC - SIMCHI nonEst NA NA NA NA
## SAL GRExMIC - SALAPI nonEst NA NA NA NA
## SAL GRExMIC - SALCLE nonEst NA NA NA NA
## SAL GRExMIC - SIMCHI nonEst NA NA NA NA
## SALAPI - SALCLE nonEst NA NA NA NA
## SALAPI - SIMCHI nonEst NA NA NA NA
## SALCLE - SIMCHI nonEst NA NA NA NA
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 10 estimatesStep 13: Effect Sizes
# Overall effect size (Cohen's d) for Nativity
pooled_sd <- sqrt(residual_var)
nativity_effect_size <- nativity_coef[1] / pooled_sd
cat("Cohen's d for Nativity effect:", round(nativity_effect_size, 3), "\n")## Cohen's d for Nativity effect: 0.23if(abs(nativity_effect_size) < 0.2) {
cat("→ SMALL effect size\n")
} else if(abs(nativity_effect_size) < 0.5) {
cat("→ MEDIUM effect size\n")
} else {
cat("→ LARGE effect size\n")
}## → MEDIUM effect size# Confidence interval
confint_nativity <- confint(model_nested, parm = "NativityNative", method = "Wald")
cat("\n95% Confidence Interval for Nativity Effect:\n")##
## 95% Confidence Interval for Nativity Effect:print(confint_nativity)## 2.5 % 97.5 %
## NativityNative -12.1208 20.70452Summary and Interpretation
Model Specification
HeightChange ~ Nativity/Species + (1|Plot)
What this model does:
- Tests the MAIN EFFECT of Nativity (Native vs Non-Native)
- Accounts for SPECIES-SPECIFIC responses within each nativity level
- Treats each species as a unique fixed effect (not generalizable)
- Accounts for spatial clustering via Plot random effect
Key Differences from Random Effects Model
Species as RANDOM Effect (Previous)
- Generalizes to other species beyond these 22
- Conservative test (more difficult to find significance)
- Appropriate when species are random sample from population
- Question: “Do native plants generally grow more?”
Species as FIXED Effect (This Analysis)
- Specific to these exact 22 species only
- Less conservative (easier to find significance)
- Appropriate when interested in these specific species
- Question: “Do THESE native species grow more than THESE non-native species?”
Which Model Should You Use?
Use RANDOM effects (previous model) if:
✓ You want to generalize beyond these species
✓ Species are a random sample from a larger population
✓ You’re asking: “Do native species in general grow differently?”
✓ More conservative, more generalizable
Use FIXED effects (this model) if:
✓ You only care about these specific species
✓ These species were deliberately chosen (not random)
✓ You’re asking: “Do these specific native species differ from these
specific non-native species?”
✓ Less conservative, less generalizable
Recommendation
For most ecological questions about nativity effects, the RANDOM effects model is more appropriate because it allows generalization. However, if your study specifically selected these species for comparison, the fixed effects model may be justified.
The choice between random and fixed effects fundamentally changes your inference scope and should be based on your research question and study design, not just which gives a smaller p-value!
Session Info
sessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.5
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Atlantic/Reykjavik
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] multcomp_1.4-29 TH.data_1.1-5 MASS_7.3-65 survival_3.8-3
## [5] mvtnorm_1.3-3 emmeans_2.0.1 car_3.1-3 carData_3.0-6
## [9] dplyr_1.1.4 ggplot2_4.0.1 lmerTest_3.2-0 lme4_1.1-38
## [13] Matrix_1.7-4 readxl_1.4.5
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.55 bslib_0.9.0
## [4] lattice_0.22-7 numDeriv_2016.8-1.1 vctrs_0.7.1
## [7] tools_4.5.2 Rdpack_2.6.4 generics_0.1.4
## [10] parallel_4.5.2 pbkrtest_0.5.5 sandwich_3.1-1
## [13] tibble_3.3.1 pkgconfig_2.0.3 RColorBrewer_1.1-3
## [16] S7_0.2.1 lifecycle_1.0.5 compiler_4.5.2
## [19] farver_2.1.2 codetools_0.2-20 htmltools_0.5.9
## [22] sass_0.4.10 yaml_2.3.12 Formula_1.2-5
## [25] pillar_1.11.1 nloptr_2.2.1 jquerylib_0.1.4
## [28] tidyr_1.3.2 cachem_1.1.0 reformulas_0.4.3.1
## [31] boot_1.3-32 abind_1.4-8 nlme_3.1-168
## [34] tidyselect_1.2.1 digest_0.6.39 purrr_1.2.1
## [37] labeling_0.4.3 splines_4.5.2 fastmap_1.2.0
## [40] grid_4.5.2 cli_3.6.5 magrittr_2.0.4
## [43] utf8_1.2.6 broom_1.0.11 withr_3.0.2
## [46] scales_1.4.0 backports_1.5.0 estimability_1.5.1
## [49] rmarkdown_2.30 cellranger_1.1.0 zoo_1.8-15
## [52] coda_0.19-4.1 evaluate_1.0.5 knitr_1.51
## [55] rbibutils_2.4 viridisLite_0.4.2 rlang_1.1.7
## [58] Rcpp_1.1.1 glue_1.8.0 rstudioapi_0.18.0
## [61] minqa_1.2.8 jsonlite_2.0.0 R6_2.6.1