Chapter 6 Alpha diversity

if ("package:pscl" %in% search()) detach("package:pscl", unload = TRUE)
if ("package:MuMIn" %in% search()) detach("package:MuMIn", unload = TRUE)
if ("package:MASS" %in% search()) detach("package:MASS", unload = TRUE)

load("data/data.Rdata")

6.1 Summary table

library(hilldiv2)
behaviour_colors <- c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c"   
)

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%  # richness
  as.matrix() %>%        # ensure it's a plain matrix for t()
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%  # neutral diversity
  as.matrix() %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%  # phylogenetic diversity
  as.matrix() %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
genome_counts_filt <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  rownames_to_column(., "genome")

genome_gifts <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts <- genome_gifts[, colSums(genome_gifts != 0) > 0]

dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  filter(genome %in% rownames(dist)) %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%  # functional diversity
  as.matrix() %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

alpha	aggr	tame	unsel
richness	38.89±45.23	29.4±39.37	51.25±45.54
neutral	16.43±18.23	14.11±17.85	22.02±19.84
phylogenetic	3.48±2.55	3.14±2.45	4.34±2.48
functional	1.29±0.22	1.29±0.24	1.37±0.21

alpha	fcolon	mcolon	kileum	dileum
richness	97.91±21.81	24.33±25.63	1.73±0.65	25.83±32.59
neutral	42.5±9.35	11.49±12.46	1.15±0.14	11.34±12.01
phylogenetic	6.54±0.75	3.3±2.38	1.11±0.11	3.03±1.91
functional	1.51±0.06	1.28±0.22	1.03±0.03	1.39±0.16

6.2 By location

6.2.1 Plots

richness_mean <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  group_by(gut_location) %>%
  dplyr::summarise_at(.vars = names(.)[2], .funs = c("Richness mean" = "mean", "Richness sd" = "sd"))

neutral_mean <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  group_by(gut_location) %>%
  dplyr::summarise_at(.vars = names(.)[3], .funs = c("Neutral mean" = "mean", "Neutral sd" = "sd"))

phylogenetic_mean <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  group_by(gut_location) %>%
  dplyr::summarise_at(.vars = names(.)[4], .funs = c("Phylogenetic mean" = "mean", "Phylogenetic sd" = "sd"))

cbind(richness_mean, neutral_mean[, 2:3], phylogenetic_mean[, 2:3])%>%
  tt()

gut_location	Richness mean	Richness sd	Neutral mean	Neutral sd	Phylogenetic mean	Phylogenetic sd
D_ileum	25.833333	32.5906661	11.335985	12.0056914	3.032274	1.9060737
F_colon	97.909091	21.8057540	42.504875	9.3501832	6.544581	0.7458711
K_ileum	1.727273	0.6466698	1.153628	0.1406374	1.112373	0.1137430
M_colon	24.333333	25.6314201	11.486763	12.4646180	3.299202	2.3761209

group_n <- alpha_div %>%
  left_join(., sample_metadata, by = join_by(sample == sample))%>%
 dplyr::select(gut_location) %>%
  pull() %>%
  unique() %>%
  length()

#pdf("figures/diversity_location.pdf",width=20, height=9)
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
  ggplot(aes(y = value, x = gut_location, group=gut_location, color=gut_location, fill=gut_location)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE, alpha = 0.5) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(values = c(
    F_colon = "#d6604d",
    M_colon = "#542788",
    D_ileum = "#fdb863",
    K_ileum = "#bb99d8")) +
  scale_fill_manual(values = c(
    F_colon = "#d6604d",
    M_colon = "#542788",
    D_ileum = "#fdb863",
    K_ileum = "#bb99d8", "50"
  )) +
  facet_wrap(. ~ metric, scales = "free", ncol=4) +
  coord_cartesian(xlim = c(1, NA)) +
  theme_classic() +
  theme(
    axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 12, color="black",face="bold"),
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        legend.text=element_text(size=10),
        legend.title = element_text(size=12))+
    guides(fill = guide_legend(override.aes = list(size=3)))

#dev.off()

6.2.2 Generalized linear models

6.2.2.1 Richness

First, we will model the richness only taking into account the gut location, as it seems to hold the greatest variance.

Count data are often modeled with Poisson distribution, but this distribution assumes the mean = variance. In the case of richness, the data is usually overdispersed (variance > mean), so we use a negative binomial model.

Analysis of Deviance Table

Model: Negative Binomial(1.314), link: log

Response: richness

Terms added sequentially (first to last)

             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                            45    123.001              
gut_location  3   76.079        42     46.922 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The p_value (< 2.2e-16) is much lower than 0.05 , so this model strongly suggests that alpha diversity differs greatly between gut locations.

#how well the model fits compared to a null model
pscl::pR2(Modelq0_Loca)

fitting null model for pseudo-r2

Theta(1) = 0.550546, 2(Ls - Lm) = 54.832600

         llh      llhNull           G2     McFadden         r2ML         r2CU 
-182.1918677 -206.5422888   48.7008421    0.1178956    0.6530983    0.6531806 

Even though gut location is highly significant, it only modestly explains richness variation on a McFadden scale

#Estimated marginal means
emmeans(Modelq0_Loca, pairwise ~ gut_location)

$emmeans
 gut_location emmean    SE  df asymp.LCL asymp.UCL
 D_ileum       3.252 0.258 Inf     2.746      3.76
 F_colon       4.584 0.265 Inf     4.065      5.10
 K_ileum       0.547 0.349 Inf    -0.138      1.23
 M_colon       3.192 0.259 Inf     2.685      3.70

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast          estimate    SE  df z.ratio p.value
 D_ileum - F_colon  -1.3324 0.370 Inf  -3.603  0.0018
 D_ileum - K_ileum   2.7051 0.434 Inf   6.231  <.0001
 D_ileum - M_colon   0.0598 0.365 Inf   0.164  0.9984
 F_colon - K_ileum   4.0375 0.438 Inf   9.216  <.0001
 F_colon - M_colon   1.3922 0.370 Inf   3.762  0.0010
 K_ileum - M_colon  -2.6453 0.434 Inf  -6.090  <.0001

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 

K_ileum has the lowest richness, F_colon the highest, and D_ileum / M_colon are intermediate.

summary(Modelq0_Loca)

Call:
MASS::glm.nb(formula = richness ~ gut_location, data = alpha_div_meta, 
    trace = TRUE, init.theta = 1.314009242, link = log)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.25167    0.25816  12.596  < 2e-16 ***
gut_locationF_colon  1.33237    0.36981   3.603 0.000315 ***
gut_locationK_ileum -2.70512    0.43412  -6.231 4.63e-10 ***
gut_locationM_colon -0.05982    0.36536  -0.164 0.869948    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(1.314) family taken to be 1)

    Null deviance: 123.001  on 45  degrees of freedom
Residual deviance:  46.922  on 42  degrees of freedom
AIC: 374.38

Number of Fisher Scoring iterations: 1

              Theta:  1.314 
          Std. Err.:  0.279 

 2 x log-likelihood:  -364.384 

6.2.2.2 Neutral

In this case neutral diversity is measured with the exponential of the Shannon entropy (Hill q =1), this means that the data is likely continuous (unlike richness, which is counts), and often has a normal distribution, therefore, we can use a linear model (lm).

Modelq1_Loca <- lm(formula = neutral ~ gut_location, data = alpha_div_meta) 
anova(Modelq1_Loca)

Analysis of Variance Table

Response: neutral
             Df Sum Sq Mean Sq F value   Pr(>F)    
gut_location  3  10650  3550.2  35.766 1.22e-11 ***
Residuals    42   4169    99.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(Modelq1_Loca)

           R2m       R2c
[1,] 0.7045247 0.7045247

emmeans(Modelq1_Loca, pairwise ~ gut_location)

$emmeans
 gut_location emmean   SE df lower.CL upper.CL
 D_ileum       11.34 2.88 42     5.53    17.14
 F_colon       42.50 3.00 42    36.44    48.57
 K_ileum        1.15 3.00 42    -4.91     7.22
 M_colon       11.49 2.88 42     5.68    17.29

Confidence level used: 0.95 

$contrasts
 contrast          estimate   SE df t.ratio p.value
 D_ileum - F_colon  -31.169 4.16 42  -7.495  <.0001
 D_ileum - K_ileum   10.182 4.16 42   2.448  0.0834
 D_ileum - M_colon   -0.151 4.07 42  -0.037  1.0000
 F_colon - K_ileum   41.351 4.25 42   9.734  <.0001
 F_colon - M_colon   31.018 4.16 42   7.458  <.0001
 K_ileum - M_colon  -10.333 4.16 42  -2.485  0.0771

P value adjustment: tukey method for comparing a family of 4 estimates 

The p-value (1.22e-11) is much lower than 0.05, suggesting that gut_location significantly affects neutral diversity.

R2m= 0.7045, which means that about 70% of the variation in neutral diversity is explained by gut location.

F_colon has the highest neutral diversity (≈42.5). K_ileum has the lowest (≈1.15). D_ileum and M_colon are similar (≈11.3–11.5).

6.2.2.3 Phylogenetic

Modelq1p_Loca <- lm(formula = phylogenetic ~ gut_location, data = alpha_div_meta) 
anova(Modelq1p_Loca)

Analysis of Variance Table

Response: phylogenetic
             Df Sum Sq Mean Sq F value    Pr(>F)    
gut_location  3 167.77  55.923  21.796 1.143e-08 ***
Residuals    42 107.76   2.566                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(Modelq1p_Loca)

           R2m       R2c
[1,] 0.5923431 0.5923431

emmeans(Modelq1p_Loca, pairwise ~ gut_location)

$emmeans
 gut_location emmean    SE df lower.CL upper.CL
 D_ileum        3.03 0.462 42    2.099     3.97
 F_colon        6.54 0.483 42    5.570     7.52
 K_ileum        1.11 0.483 42    0.138     2.09
 M_colon        3.30 0.462 42    2.366     4.23

Confidence level used: 0.95 

$contrasts
 contrast          estimate    SE df t.ratio p.value
 D_ileum - F_colon   -3.512 0.669 42  -5.253  <.0001
 D_ileum - K_ileum    1.920 0.669 42   2.871  0.0312
 D_ileum - M_colon   -0.267 0.654 42  -0.408  0.9767
 F_colon - K_ileum    5.432 0.683 42   7.953  <.0001
 F_colon - M_colon    3.245 0.669 42   4.854  0.0001
 K_ileum - M_colon   -2.187 0.669 42  -3.271  0.0111

P value adjustment: tukey method for comparing a family of 4 estimates 

6.2.2.4 Functional

Modelq1F_Loca <- lm(formula = functional ~ gut_location, data = alpha_div_meta) 
anova(Modelq1F_Loca)

Analysis of Variance Table

Response: functional
             Df  Sum Sq Mean Sq F value    Pr(>F)    
gut_location  3 1.35326 0.45109  22.264 8.729e-09 ***
Residuals    42 0.85094 0.02026                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(Modelq1F_Loca)

           R2m       R2c
[1,] 0.5974711 0.5974711

emmeans(Modelq1F_Loca, pairwise ~ gut_location)

$emmeans
 gut_location emmean     SE df lower.CL upper.CL
 D_ileum        1.39 0.0411 42    1.311     1.48
 F_colon        1.51 0.0429 42    1.418     1.59
 K_ileum        1.03 0.0429 42    0.947     1.12
 M_colon        1.28 0.0411 42    1.195     1.36

Confidence level used: 0.95 

$contrasts
 contrast          estimate     SE df t.ratio p.value
 D_ileum - F_colon   -0.111 0.0594 42  -1.866  0.2582
 D_ileum - K_ileum    0.360 0.0594 42   6.066  <.0001
 D_ileum - M_colon    0.116 0.0581 42   1.992  0.2072
 F_colon - K_ileum    0.471 0.0607 42   7.765  <.0001
 F_colon - M_colon    0.227 0.0594 42   3.814  0.0024
 K_ileum - M_colon   -0.245 0.0594 42  -4.118  0.0010

P value adjustment: tukey method for comparing a family of 4 estimates 

6.3 By behaviour and location

6.3.1 Plots

6.3.1.1 Richness

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(sample_metadata, by = "sample") %>%
  filter(!is.na(gut_location)) %>%
  filter(metric=="richness") %>%
      ggplot(aes(y = value, x = fox_behaviour, group=fox_behaviour, color=fox_behaviour, fill=fox_behaviour)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5) +
      geom_jitter(width = 0.1, alpha=0.5) +
      scale_color_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      scale_fill_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      facet_wrap(. ~ gut_location, scales = "fixed", ncol=6) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )+ labs(title= "Richness by gut location and fox behaviour")

6.3.1.2 Neutral

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(sample_metadata, by = "sample") %>%
  filter(!is.na(gut_location)) %>%
  filter(metric=="neutral") %>%
      ggplot(aes(y = value, x = fox_behaviour, group=fox_behaviour, color=fox_behaviour, fill=fox_behaviour)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5) +
      geom_jitter(width = 0.1, alpha=0.5) +
      scale_color_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      scale_fill_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      facet_wrap(. ~ gut_location, scales = "fixed", ncol=6) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )+ labs(title= "Neutral diversity by gut location and fox behaviour")

6.3.1.3 Phylogenetic

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(sample_metadata, by = "sample") %>%
  filter(!is.na(gut_location)) %>%
  filter(metric=="phylogenetic") %>%
      ggplot(aes(y = value, x = fox_behaviour, group=fox_behaviour, color=fox_behaviour, fill=fox_behaviour)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5) +
      geom_jitter(width = 0.1, alpha=0.5) +
      scale_color_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      scale_fill_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      facet_wrap(. ~ gut_location, scales = "fixed", ncol=6) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )+ labs(title= "Phylogenetic diversity by gut location and fox behaviour")

6.3.1.4 Functional

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(sample_metadata, by = "sample") %>%
  filter(!is.na(gut_location)) %>%
  filter(metric=="functional") %>%
      ggplot(aes(y = value, x = fox_behaviour, group=fox_behaviour, color=fox_behaviour, fill=fox_behaviour)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5) +
      geom_jitter(width = 0.1, alpha=0.5) +
      scale_color_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      scale_fill_manual(values =  c(
  tame = "#1f77b4",   
  aggr = "#d62728",   
  unsel = "#2ca02c")) +
      facet_wrap(. ~ gut_location, scales = "fixed", ncol=6) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )+ labs(title= "Functional diversity by gut location and fox behaviour")

6.3.2 Mixed models

6.3.2.1 Richness

Model_richness <- glmer.nb(richness ~ fox_behaviour+(1|gut_location), data = alpha_div_meta)
summary(Model_richness)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(1.3513)  ( log )
Formula: richness ~ fox_behaviour + (1 | gut_location)
   Data: alpha_div_meta

      AIC       BIC    logLik -2*log(L)  df.resid 
    386.7     395.8    -188.3     376.7        41 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0886 -0.6944 -0.1725  0.3584  3.2176 

Random effects:
 Groups       Name        Variance Std.Dev.
 gut_location (Intercept) 2.059    1.435   
Number of obs: 46, groups:  gut_location, 4

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.0142     0.7503   4.017 5.89e-05 ***
fox_behaviourtame   -0.4865     0.3047  -1.596    0.110    
fox_behaviourunsel   0.3134     0.3846   0.815    0.415    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) fx_bhvrt
fox_behvrtm -0.207         
fox_bhvrnsl -0.161  0.394  

emmeans(Model_richness, pairwise ~ fox_behaviour)

$emmeans
 fox_behaviour emmean    SE  df asymp.LCL asymp.UCL
 aggr            3.01 0.750 Inf      1.54      4.48
 tame            2.53 0.749 Inf      1.06      4.00
 unsel           3.33 0.786 Inf      1.79      4.87

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast     estimate    SE  df z.ratio p.value
 aggr - tame     0.486 0.305 Inf   1.596  0.2472
 aggr - unsel   -0.313 0.385 Inf  -0.815  0.6938
 tame - unsel   -0.800 0.385 Inf  -2.077  0.0946

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

There is no significant effect of fox behaviour on the richness.

6.3.2.2 Neutral

Model_neutral <- lme(fixed = neutral ~ fox_behaviour, data = alpha_div_meta,
               random = ~ 1 | gut_location)#log(seq_depth)+
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  346.8649 355.6709 -168.4324

Random effects:
 Formula: ~1 | gut_location
        (Intercept) Residual
StdDev:     17.6769 9.756911

Fixed effects:  neutral ~ fox_behaviour 
                       Value Std.Error DF    t-value p-value
(Intercept)        16.993796  9.134109 40  1.8604766  0.0702
fox_behaviourtame  -2.880201  3.173896 40 -0.9074655  0.3696
fox_behaviourunsel  5.023322  4.148905 40  1.2107584  0.2331
 Correlation: 
                   (Intr) fx_bhvrt
fox_behaviourtame  -0.183         
fox_behaviourunsel -0.140  0.404  

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.62821779 -0.78911686 -0.04122328  0.46700583  2.23438491 

Number of Observations: 46
Number of Groups: 4 

6.3.2.3 Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ fox_behaviour, data = alpha_div_meta,
               random = ~ 1 | gut_location)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  188.6297 197.4357 -89.31487

Random effects:
 Formula: ~1 | gut_location
        (Intercept) Residual
StdDev:     2.20058 1.577186

Fixed effects:  phylogenetic ~ fox_behaviour 
                       Value Std.Error DF    t-value p-value
(Intercept)         3.518775  1.161669 40  3.0290686  0.0043
fox_behaviourtame  -0.380035  0.513043 40 -0.7407461  0.4632
fox_behaviourunsel  0.822181  0.670654 40  1.2259398  0.2274
 Correlation: 
                   (Intr) fx_bhvrt
fox_behaviourtame  -0.233         
fox_behaviourunsel -0.178  0.404  

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.47646378 -0.74894597  0.02329104  0.31903019  2.12465637 

Number of Observations: 46
Number of Groups: 4 

6.3.2.4 Functional

Model_func <- lme(fixed = functional ~ fox_behaviour, data = alpha_div_meta,
               random = ~ 1 | gut_location)
summary(Model_func)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -18.20571 -9.399711 14.10286

Random effects:
 Formula: ~1 | gut_location
        (Intercept)  Residual
StdDev:    0.197249 0.1424272

Fixed effects:  functional ~ fox_behaviour 
                       Value  Std.Error DF   t-value p-value
(Intercept)        1.2866404 0.10420652 40 12.347024  0.0000
fox_behaviourtame  0.0054055 0.04633013 40  0.116674  0.9077
fox_behaviourunsel 0.0799250 0.06056314 40  1.319697  0.1944
 Correlation: 
                   (Intr) fx_bhvrt
fox_behaviourtame  -0.235         
fox_behaviourunsel -0.179  0.404  

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.3122951 -0.4504445 -0.0874495  0.5457938  1.8860386 

Number of Observations: 46
Number of Groups: 4 

detach("package:pscl", unload = TRUE)
detach("package:MuMIn", unload = TRUE)
detach("package:MASS", unload = TRUE)

Warning: namespace 'MASS' no puede ser descargado:
  namespace 'MASS' is imported by 'lmerTest', 'clusterGeneration', 'vegan', 'ANCOMBC', 'lme4', 'multtest', 'mia', 'geiger', 'DescTools', 'ade4', 'TH.data', 'phytools' so cannot be unloaded

detach("package:hilldiv2", unload = TRUE)