LMEs for Sable Physiology Comparisons

Linear Mixed Effect Models

Michelle Shero, April 3, 2024

# Load and attach files (change to whatever directory your file is in; note I saved the file as a .csv)
setwd("~/Desktop")
phys<- read.csv("HG_Physiology_R_Data.csv")
attach(phys)

The following objects are masked from phys (pos = 3):

    Age, Animal_ID_corrected, Animal_ID_raw, BV, BV_., Handling, Hb, Hct, Mass, MCH_pg, MCHC,
    MCV_fL, PV, PV_., RBC_mL, Sex, Treatment, Treatment.Num, WBC, Year

The following objects are masked from phys (pos = 4):

    Age, Animal_ID_corrected, Animal_ID_raw, BV, BV_., Handling, Hb, Hct, Mass, MCH_pg, MCHC,
    MCV_fL, PV, PV_., RBC_mL, Sex, Treatment, Treatment.Num, WBC, Year

The following objects are masked from phys (pos = 24):

    Age, Animal_ID_corrected, Animal_ID_raw, BV, BV_., Handling, Hb, Hct, Mass, MCH_pg, MCHC,
    MCV_fL, PV, PV_., RBC_mL, Sex, Treatment, WBC, Year

head(phys) #note that for some of your columns, like MCHC it thinks it is a character string. You'l have to do: as.numeric(MCHC) before you run those.

library(lme4)
library(lmerTest)
library(ggplot2)
library(ggfortify)
library(lsmeans)
library(predictmeans)
library(HLMdiag)
library(multcomp)
library(plyr)
library(psych)
library(MuMIn)

Explore

# Dotplots
# Useful for looking for potential outliers. For example, you can see there's one adult that looks like it has a low Hb value that sticks out. Don't remove them yet though --- we'll only remove points if they are still outliers when we look at how the models fit the data.
dotchart(Hb[which(Age == "Adult")], main = "Adult")

dotchart(Hb[which(Age == "Pup")], main = "Pup")

# Test for Normality
# Note that it is good to look at this and often the physiology data will be nice and normally distributed. But this isn't an actual requirement for LME's. The requirement is that the RESIDUALS from the model itself are normally distributed.
shapiro.test(Hb[which(Age == "Adult" & Handling == "A")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Adult" & Handling == "A")]
W = 0.95286, p-value = 0.1287

shapiro.test(Hb[which(Age == "Adult" & Handling == "B")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Adult" & Handling == "B")]
W = 0.95862, p-value = 0.2216

shapiro.test(Hb[which(Age == "Pup" & Handling == "A")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Pup" & Handling == "A")]
W = 0.97715, p-value = 0.6485

shapiro.test(Hb[which(Age == "Pup" & Handling == "B")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Pup" & Handling == "B")]
W = 0.95779, p-value = 0.2095

shapiro.test(Hb[which(Age == "Pup" & Handling == "C")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Pup" & Handling == "C")]
W = 0.96207, p-value = 0.2954

shapiro.test(Hb[which(Age == "Pup" & Handling == "D")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Pup" & Handling == "D")]
W = 0.971, p-value = 0.7129

shapiro.test(Hb[which(Age == "Pup" & Handling == "E")])

    Shapiro-Wilk normality test

data:  Hb[which(Age == "Pup" & Handling == "E")]
W = 0.96424, p-value = 0.8517

# Summary Statistics
describeBy(Hb[which(Age == "Adult")], group = Handling[which(Age == "Adult")], na.rm = TRUE)

 Descriptive statistics by group 
group: A

------------------------------------------------------------------------------------ 
group: B

describeBy(Hb[which(Age == "Pup")], group = Handling[which(Age == "Pup")], na.rm = TRUE)

 Descriptive statistics by group 
group: A

------------------------------------------------------------------------------------ 
group: B

------------------------------------------------------------------------------------ 
group: C

------------------------------------------------------------------------------------ 
group: D

------------------------------------------------------------------------------------ 
group: E

Null model with just individual effect

Hb.pup<- Hb[which(Age == "Pup")]
Handling.pup<- Handling[which(Age == "Pup")]
Mass.pup<- Mass[which(Age == "Pup")]
Sex.pup<- Sex[which(Age == "Pup")]
Year.pup<- as.factor(Year[which(Age == "Pup")]) #note am putting this in as a factor --- so it becomes a category instead of a continuous numeric variable.
Treatment.pup<- Treatment[which(Age == "Pup")]
ID.pup<- Animal_ID_corrected[which(Age == "Pup")]

#Null model just looking at animal ID
g0<- lmer(Hb.pup ~ (1|ID.pup))

boundary (singular) fit: see help('isSingular')

AICc(g0)

[1] 507.315

summary(g0)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ (1 | ID.pup)

REML criterion at convergence: 501.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5649 -0.6140 -0.1430  0.6141  2.6334 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.000    0.000   
 Residual             2.586    1.608   
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    16.11       0.14 131.00   115.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

anova(g0, test = "F")

Type III Analysis of Variance Table with Satterthwaite's method
     Sum Sq Mean Sq NumDF DenDF F value Pr(>F)

#plot(g0)

Effect of Handling

#Add in Handling for LME
g1<- lmer(Hb.pup ~ Handling.pup + (1|ID.pup))
AICc(g1)

[1] 429.8784

summary(g1) #This gives you the results for the overall model. I.e., you can see that handling time has a significant effect. But you can't see the pairwise comparisons here (is it A vs B that are significant? A vs C etc etc)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ Handling.pup + (1 | ID.pup)

REML criterion at convergence: 415

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.86038 -0.47695  0.04363  0.55929  2.46680 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.456    0.6753  
 Residual             1.053    1.0259  
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
              Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)    15.9711     0.2047  96.1169  78.019  < 2e-16 ***
Handling.pupB  -0.8485     0.2464  88.9368  -3.444 0.000877 ***
Handling.pupC  -0.4050     0.2487  89.5331  -1.628 0.107003    
Handling.pupD   1.8750     0.2799  93.8881   6.699 1.54e-09 ***
Handling.pupE   3.3211     0.4781 100.4499   6.947 3.79e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Hndl.B Hndl.C Hndl.D
Handlng.ppB -0.580                     
Handlng.ppC -0.574  0.486              
Handlng.ppD -0.510  0.431  0.431       
Handlng.ppE -0.299  0.253  0.252  0.244

comp.hb <- glht(g1, linfct=mcp(Handling.pup="Tukey"))
print(summary(comp.hb, test = adjusted(type = "bonferroni"))) # Using Bonferroni adjustment, for some reason I've never been able to figure out, you can't do it in first line.

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Hb.pup ~ Handling.pup + (1 | ID.pup))

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
B - A == 0  -0.8485     0.2464  -3.444  0.00574 ** 
C - A == 0  -0.4050     0.2487  -1.628  1.00000    
D - A == 0   1.8750     0.2799   6.699 2.10e-10 ***
E - A == 0   3.3211     0.4781   6.947 3.74e-11 ***
C - B == 0   0.4435     0.2511   1.766  0.77385    
D - B == 0   2.7236     0.2820   9.656  < 2e-16 ***
E - B == 0   4.1696     0.4793   8.699  < 2e-16 ***
D - C == 0   2.2800     0.2833   8.048 8.88e-15 ***
E - C == 0   3.7261     0.4801   7.761 8.44e-14 ***
E - D == 0   1.4461     0.4917   2.941  0.03270 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)

#Ultimately, what you have now is a linear mixed effect models, with Bonferroni posthoc comparisons
r.squaredGLMM(g1)

           R2m       R2c
[1,] 0.4677032 0.6286158

anova(g1, test = "F")

Type III Analysis of Variance Table with Satterthwaite's method
             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
Handling.pup 159.41  39.853     4 92.734  37.863 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(g1) # You typically want to make sure that there's no trend in this plot. But it's kinda silly when all you have is a grouping factor in the model. See later models where we put in continuous variables like mass. More important in this case to see that residuals are normally distributed

#Assess fit
model<- g1

#Residuals versus fitted should show no relationship
par(mfrow = c(1,1))

plot(resid(model)~fitted(model))
summary((lm(resid(model)~fitted(model)))) #confirm that there is no relationship

Call:
lm(formula = resid(model) ~ fitted(model))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8880 -0.5655  0.1008  0.5703  2.4701 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -1.93788    1.09337  -1.772   0.0787 .
fitted(model)  0.12025    0.06767   1.777   0.0779 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9157 on 130 degrees of freedom
Multiple R-squared:  0.02372,   Adjusted R-squared:  0.01621 
F-statistic: 3.158 on 1 and 130 DF,  p-value: 0.07789

#Cooks distance looks for influential values
CookD(model)  # Marked points show influential values from the model. 

plot(model, rstudent(.) ~ hatvalues(.)) # Residuals versus Leverage (the hat values). This is a better way to look for outliers. Wouldn't worry about this until we've finished all the field work though and have everything in hand. Nothing here looks glaring, would be surprised if you end up with much of anything with outliers

#Compare distribution of dataset to a normal distribution
qqnorm(resid(model), pch = 16)
qqline(resid(model)) #Points are a pretty good fit to the line except at the ends

#Want the model residuals to be approximately normally distributed
plot(density(resid(model)))

Effect of Iron Supplement

#Add in Handling for LME
g2<- lmer(Hb.pup ~ Treatment.pup + (1|ID.pup))

boundary (singular) fit: see help('isSingular')

AICc(g2)

[1] 509.1382

summary(g2) #This gives you the results for the overall model. I.e., you can see that handling time has a significant effect. But you can't see the pairwise comparisons here (is it A vs B that are significant? A vs C etc etc)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ Treatment.pup + (1 | ID.pup)

REML criterion at convergence: 498.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6204 -0.6560 -0.1097  0.5324  2.5972 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.000    0.000   
 Residual             2.567    1.602   
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           16.18857    0.30280 129.00000  53.462   <2e-16 ***
Treatment.pup20_iron  -0.55665    0.43639 129.00000  -1.276    0.204    
Treatment.pupsham      0.06092    0.35299 129.00000   0.173    0.863    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Tr.20_
Trtmnt.p20_ -0.694       
Trtmnt.ppsh -0.858  0.595
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

comp.hb <- glht(g2, linfct=mcp(Treatment.pup="Tukey"))
print(summary(comp.hb, test = adjusted(type = "bonferroni"))) # Using Bonferroni adjustment, for some reason I've never been able to figure out, you can't do it in first line.

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Hb.pup ~ Treatment.pup + (1 | ID.pup))

Linear Hypotheses:
                       Estimate Std. Error z value Pr(>|z|)
20_iron - 10_iron == 0 -0.55665    0.43639  -1.276    0.606
sham - 10_iron == 0     0.06092    0.35299   0.173    1.000
sham - 20_iron == 0     0.61756    0.36285   1.702    0.266
(Adjusted p values reported -- bonferroni method)

#Ultimately, what you have now is a linear mixed effect models, with Bonferroni posthoc comparisons
r.squaredGLMM(g2)

            R2m        R2c
[1,] 0.02218296 0.02218296

anova(g2, test = "F")

Type III Analysis of Variance Table with Satterthwaite's method
              Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Treatment.pup 7.6298  3.8149     2   129  1.4859 0.2301

plot(g2) # You can see here how you can't really look for no trend because all you have is treatment groups. More important in this case to see that residuals are normally distributed

#Assess fit
model<- g2

#Residuals versus fitted should show no relationship
par(mfrow = c(1,1))

plot(resid(model)~fitted(model))

summary((lm(resid(model)~fitted(model)))) #confirm that there is no relationship

Call:
lm(formula = resid(model) ~ fitted(model))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1986 -1.0511 -0.1757  0.8530  4.1614 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)    1.854e-12  9.313e+00       0        1
fitted(model) -1.151e-13  5.778e-01       0        1

Residual standard error: 1.596 on 130 degrees of freedom
Multiple R-squared:  3.081e-28, Adjusted R-squared:  -0.007692 
F-statistic: 4.005e-26 on 1 and 130 DF,  p-value: 1

#Cooks distance looks for influential values
#CookD(model)  # Marked points show influential values from the model. 
plot(model, rstudent(.) ~ hatvalues(.)) # Residuals versus Leverage (the hat values). This is a better way to look for outliers. Wouldn't worry about this until we've finished all the field work though and have everything in hand. Nothing here looks glaring, would be surprised if you end up with much of anything with outliers

#Compare distribution of dataset to a normal distribution
qqnorm(resid(model), pch = 16)
qqline(resid(model)) #Points are a pretty good fit to the line except at the ends

#Want the model residuals to be approximately normally distributed
plot(density(resid(model)))

Effect of Handling + Iron Supplement

#Add in Handling for LME
g3<- lmer(Hb.pup ~ Handling.pup + Treatment.pup + (1|ID.pup))
AICc(g3)

[1] 429.3969

summary(g3) #This gives you the results for the overall model. I.e., you can see that handling time has a significant effect. But you can't see the pairwise comparisons here (is it A vs B that are significant? A vs C etc etc)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ Handling.pup + Treatment.pup + (1 | ID.pup)

REML criterion at convergence: 409.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.92488 -0.56737  0.05057  0.59977  2.39436 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.4017   0.6338  
 Residual             1.0437   1.0216  
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           15.7195     0.3477  40.3754  45.210  < 2e-16 ***
Handling.pupB         -0.8476     0.2453  90.7328  -3.455 0.000839 ***
Handling.pupC         -0.4105     0.2476  91.3209  -1.658 0.100760    
Handling.pupD          1.9088     0.2789  95.1609   6.844 7.45e-10 ***
Handling.pupE          3.3132     0.4771 100.7711   6.945 3.78e-10 ***
Treatment.pup20_iron  -0.2919     0.4465  27.5675  -0.654 0.518673    
Treatment.pupsham      0.5047     0.3638  27.3071   1.387 0.176613    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Hndl.B Hndl.C Hndl.D Hndl.E Tr.20_
Handlng.ppB -0.325                                   
Handlng.ppC -0.325  0.485                            
Handlng.ppD -0.327  0.430  0.429                     
Handlng.ppE -0.248  0.250  0.250  0.243              
Trtmnt.p20_ -0.644 -0.019 -0.004 -0.004  0.090       
Trtmnt.ppsh -0.793 -0.016 -0.017  0.045  0.079  0.616

comp.hb <- glht(g3, linfct=mcp(Treatment.pup="Tukey"))
print(summary(comp.hb, test = adjusted(type = "bonferroni"))) # Using Bonferroni adjustment, for some reason I've never been able to figure out, you can't do it in first line.

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Hb.pup ~ Handling.pup + Treatment.pup + (1 | ID.pup))

Linear Hypotheses:
                       Estimate Std. Error z value Pr(>|z|)  
20_iron - 10_iron == 0  -0.2919     0.4465  -0.654   1.0000  
sham - 10_iron == 0      0.5047     0.3638   1.387   0.4962  
sham - 20_iron == 0      0.7965     0.3627   2.196   0.0843 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)

#Ultimately, what you have now is a linear mixed effect models, with Bonferroni posthoc comparisons
r.squaredGLMM(g3)

           R2m       R2c
[1,] 0.4915809 0.6328632

anova(g3, test = "F")

Type III Analysis of Variance Table with Satterthwaite's method
               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)    
Handling.pup  160.855  40.214     4 94.012 38.5286 <2e-16 ***
Treatment.pup   5.835   2.917     2 27.843  2.7952 0.0783 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(g3) # You can see here how you can't really look for no trend because all you have is treatment groups. More important in this case to see that residuals are normally distributed

#Assess fit
model<- g3

#Residuals versus fitted should show no relationship
par(mfrow = c(1,1))

plot(resid(model)~fitted(model))

summary((lm(resid(model)~fitted(model)))) #confirm that there is no relationship

Call:
lm(formula = resid(model) ~ fitted(model))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.99057 -0.61468  0.04065  0.67509  2.38505 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)    -1.6722     1.0762  -1.554    0.123
fitted(model)   0.1038     0.0666   1.558    0.122

Residual standard error: 0.915 on 130 degrees of freedom
Multiple R-squared:  0.01833,   Adjusted R-squared:  0.01078 
F-statistic: 2.427 on 1 and 130 DF,  p-value: 0.1217

#Cooks distance looks for influential values
#CookD(model)  # Marked points show influential values from the model. 
plot(model, rstudent(.) ~ hatvalues(.)) # Residuals versus Leverage (the hat values). This is a better way to look for outliers. Wouldn't worry about this until we've finished all the field work though and have everything in hand. Nothing here looks glaring, would be surprised if you end up with much of anything with outliers

#Compare distribution of dataset to a normal distribution
qqnorm(resid(model), pch = 16)
qqline(resid(model)) #Points are a pretty good fit to the line except at the ends

#Want the model residuals to be approximately normally distributed
plot(density(resid(model)))

Effect of Handling + Iron Supplement + an INTERACTIVE effect between the two

#Add in Handling for LME
g4<- lmer(Hb.pup ~ Handling.pup + Treatment.pup + (Handling.pup * Treatment.pup) +  (1|ID.pup))

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

AICc(g4)

[1] 425.2787

summary(g4) #This gives you the results for the overall model. I.e., you can see that handling time has a significant effect. 

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ Handling.pup + Treatment.pup + (Handling.pup * Treatment.pup) +      (1 | ID.pup)

REML criterion at convergence: 388.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2819 -0.4574  0.0229  0.5465  2.2746 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.3572   0.5976  
 Residual             0.9839   0.9919  
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
                                   Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                        15.64143    0.43770 95.74365  35.736  < 2e-16 ***
Handling.pupB                      -1.20095    0.55747 86.63398  -2.154  0.03400 *  
Handling.pupC                       0.02072    0.55747 86.63398   0.037  0.97044    
Handling.pupD                       1.84405    0.55747 86.63398   3.308  0.00137 ** 
Handling.pupE                       4.08105    0.70642 91.34564   5.777 1.04e-07 ***
Treatment.pup20_iron               -0.90857    0.61900 95.74365  -1.468  0.14543    
Treatment.pupsham                   0.82857    0.50253 95.74365   1.649  0.10247    
Handling.pupB:Treatment.pup20_iron  0.97952    0.76934 84.33245   1.273  0.20645    
Handling.pupC:Treatment.pup20_iron  0.44134    0.78755 85.72895   0.560  0.57667    
Handling.pupD:Treatment.pup20_iron  1.53800    0.78755 85.72895   1.953  0.05410 .  
Handling.pupB:Treatment.pupsham     0.25218    0.63476 85.86436   0.397  0.69214    
Handling.pupC:Treatment.pupsham    -0.82092    0.63476 85.86436  -1.293  0.19939    
Handling.pupD:Treatment.pupsham    -0.55823    0.67299 87.99309  -0.829  0.40908    
Handling.pupE:Treatment.pupsham    -1.77943    0.95656 93.60082  -1.860  0.06599 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 14 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

#It gives a warning. This one doesn't impact the validity of the model --- it's saying that our number of things we are comparing is getting higher while we still have a moderately low sample size. This may just go away after we get more samples next year. And if not, there are some extra simulations you can do (would need to figure out. i've never had to do those)
anova(g4, test = "F")

Missing cells for: Handling.pupE:Treatment.pup20_iron.  
Interpret type III hypotheses with care.

Type III Analysis of Variance Table with Satterthwaite's method
                            Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
Handling.pup               160.509  40.127     4 86.554 40.7843 < 2e-16 ***
Treatment.pup                2.528   1.264     2 30.746  1.2846 0.29121    
Handling.pup:Treatment.pup  15.513   2.216     7 85.716  2.2525 0.03748 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# I don't think you can do a pairwise comparison for the interactive effect
#comp.hb <- glht(g4, linfct=mcp((Handling.pup)="Tukey"))
#print(summary(comp.hb, test = adjusted(type = "bonferroni"))) # Using Bonferroni adjustment, for some reason I've never been able to figure out, you can't do it in first line.
#You can also see that it is showing you some pairwise comparisons in the model summary. But it's doing it in alphabetical order so all comparisons are made to the pup10_iron group --- which kinda sucks. Would much rather see it all compared to the sham group. So I've added a column to your spreadsheet in order to do that.

Treatment.Num.pup<- as.factor(Treatment.Num[which(Age == "Pup")])
g4.2<- lmer(Hb.pup ~ Handling.pup + Treatment.Num.pup + (Handling.pup * Treatment.Num.pup) +  (1|ID.pup))

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

summary(g4.2) #Now you can see all the treatments that are significantly different from Sham! :) 

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hb.pup ~ Handling.pup + Treatment.Num.pup + (Handling.pup * Treatment.Num.pup) +      (1 | ID.pup)

REML criterion at convergence: 388.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2819 -0.4574  0.0229  0.5465  2.2746 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID.pup   (Intercept) 0.3572   0.5976  
 Residual             0.9839   0.9919  
Number of obs: 132, groups:  ID.pup, 36

Fixed effects:
                                  Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                        16.4700     0.2469 95.7436  66.708  < 2e-16 ***
Handling.pupB                      -0.9488     0.3035 83.2786  -3.126 0.002443 ** 
Handling.pupC                      -0.8002     0.3035 83.2786  -2.636 0.010000 *  
Handling.pupD                       1.2858     0.3770 90.9645   3.411 0.000968 ***
Handling.pupE                       2.3016     0.6450 96.2755   3.569 0.000562 ***
Treatment.Num.pup10                -0.8286     0.5025 95.7436  -1.649 0.102467    
Treatment.Num.pup20                -1.7371     0.5025 95.7436  -3.457 0.000817 ***
Handling.pupB:Treatment.Num.pup10  -0.2522     0.6348 85.8644  -0.397 0.692144    
Handling.pupC:Treatment.Num.pup10   0.8209     0.6348 85.8644   1.293 0.199385    
Handling.pupD:Treatment.Num.pup10   0.5582     0.6730 87.9931   0.829 0.409076    
Handling.pupE:Treatment.Num.pup10   1.7794     0.9566 93.6008   1.860 0.065990 .  
Handling.pupB:Treatment.Num.pup20   0.7273     0.6109 82.1689   1.191 0.237271    
Handling.pupC:Treatment.Num.pup20   1.2623     0.6337 84.4670   1.992 0.049623 *  
Handling.pupD:Treatment.Num.pup20   2.0962     0.6720 86.7522   3.119 0.002460 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 14 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

anova(g4.2, test = "F")

Missing cells for: Handling.pupE:Treatment.Num.pup20.  
Interpret type III hypotheses with care.

Type III Analysis of Variance Table with Satterthwaite's method
                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
Handling.pup                   160.509  40.127     4 86.554 40.7843 < 2e-16 ***
Treatment.Num.pup                2.528   1.264     2 30.639  1.2846 0.29126    
Handling.pup:Treatment.Num.pup  15.513   2.216     7 85.722  2.2525 0.03748 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

r.squaredGLMM(g4)

           R2m       R2c
[1,] 0.5267112 0.6527653

anova(g4, test = "F")

Missing cells for: Handling.pupE:Treatment.pup20_iron.  
Interpret type III hypotheses with care.

Type III Analysis of Variance Table with Satterthwaite's method
                            Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
Handling.pup               160.509  40.127     4 86.554 40.7843 < 2e-16 ***
Treatment.pup                2.528   1.264     2 30.746  1.2846 0.29121    
Handling.pup:Treatment.pup  15.513   2.216     7 85.716  2.2525 0.03748 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(g4) # You can see here how you can't really look for no trend because all you have is treatment groups. More important in this case to see that residuals are normally distributed

#Assess fit
model<- g4

#Residuals versus fitted should show no relationship
par(mfrow = c(1,1))

plot(resid(model)~fitted(model))

summary((lm(resid(model)~fitted(model)))) #confirm that there is no relationship

Call:
lm(formula = resid(model) ~ fitted(model))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2807 -0.5068  0.0072  0.5622  2.1880 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   -1.41070    0.97363  -1.449    0.150
fitted(model)  0.08754    0.06024   1.453    0.149

Residual standard error: 0.862 on 130 degrees of freedom
Multiple R-squared:  0.01599,   Adjusted R-squared:  0.008416 
F-statistic: 2.112 on 1 and 130 DF,  p-value: 0.1486

#Cooks distance looks for influential values
#CookD(model)  # Marked points show influential values from the model. 
plot(model, rstudent(.) ~ hatvalues(.)) # Residuals versus Leverage (the hat values). This is a better way to look for outliers. Wouldn't worry about this until we've finished all the field work though and have everything in hand. Nothing here looks glaring, would be surprised if you end up with much of anything with outliers

#Compare distribution of dataset to a normal distribution
qqnorm(resid(model), pch = 16)
qqline(resid(model)) #Points are a pretty good fit to the line except at the ends

#Want the model residuals to be approximately normally distributed
plot(density(resid(model)))

Can start playing around with models and looking for which is best fit

g5<- lmer(Hb.pup ~ Mass.pup +  (1|ID.pup))
g6<- lmer(Hb.pup ~ Year.pup +  (1|ID.pup))

boundary (singular) fit: see help('isSingular')

g7<- lmer(Hb.pup ~ Handling.pup + Treatment.pup + Year.pup +  (1|ID.pup))

##AICc comparisons

# AICc is the Akaike Information Criterion corrected for small sample size (the lowercase 'c'). There's really no reason just not to always use AICc instead of AIC. If the AICc decreases by more than 2, then it significantly improves the model fit to add that variable

AICc(g0)

[1] 507.315

AICc(g1)

[1] 429.8784

AICc(g2)

[1] 509.1382

AICc(g3)

[1] 429.3969

AICc(g4)

[1] 425.2787

AICc(g5)

[1] 474.8472

AICc(g6)

[1] 510.5345

AICc(g7)

[1] 427.1819

# Based on this, the g4 model with AICc 425 is the best fit model. That's the one with Handling + Treatment + Handling x Treatment interactive effect