Mixed model
Kinde
2025-01-15
#mixed effect 2023
library(lme4)## Loading required package: Matrix#install.packages("lmerTest")
#install.packages("lmtest")
library(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric#install.packages("TMB", type = "source", repos = "https://cran.r-project.org/")
#install.packages("glmmTMB", type = "source")
#install.packages("robustlmm")
library(glmmTMB)
library(ordinal)
library(lme4)
library(robustlmm)#Data
TFMSA_NO_NA <- read.csv("/Users/kindenhk/Documents/TFMSA/TFMSA_ANOVA_DATA_NO_NA.csv")
TFMSA_NO_NA$Treat <- as.factor(TFMSA_NO_NA$Treat)
TFMSA_NO_NA$Maturity <- as.factor(TFMSA_NO_NA$Maturity)
TFMSA_NO_NA$Pedigree <- as.factor(TFMSA_NO_NA$Pedigree)
TFMSA_NO_NA$Rep <- as.factor(TFMSA_NO_NA$Rep)
#Germination data
TFMSA_Germ <- read.csv("/Users/kindenhk/Documents/TFMSA/germination%_2023.csv")
TFMSA_Germ$Rep <- factor(TFMSA_Germ$Rep)
TFMSA_Germ$Treat <- factor(TFMSA_Germ$Treat)
TFMSA_Germ$Pedigree <- factor(TFMSA_Germ$Pedigree)
TFMSA_Germ$Germ <- as.numeric(TFMSA_Germ$Germ)
TFMSA_Germ <- na.omit(TFMSA_Germ)#Linear mixed-effects model
# Fit a Linear mixed-effects model
model_ML_DTF <- lmer(DTF ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
#DTF
summary(model_ML_DTF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DTF ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## REML criterion at convergence: 2166.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2802 -0.5810 -0.1228 0.6053 4.0373
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 21.086 4.592
## Rep (Intercept) 2.364 1.538
## Residual 20.605 4.539
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 58.4778 1.3686 22.0572 42.727 < 2e-16 ***
## Treat15 1.3000 0.6767 335.0000 1.921 0.055563 .
## Treat45 2.5889 0.6767 335.0000 3.826 0.000155 ***
## Treat75 4.1889 0.6767 335.0000 6.190 1.75e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.247
## Treat45 -0.247 0.500
## Treat75 -0.247 0.500 0.500plot(model_ML_DTF)
anova(model_ML_DTF)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treat 866.39 288.8 3 335 14.015 1.253e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VarCorr(model_ML_DTF)## Groups Name Std.Dev.
## Pedigree (Intercept) 4.5919
## Rep (Intercept) 1.5375
## Residual 4.5393# Estimated marginal means (treatment differences)
posthoc_DTF <- emmeans(model_ML_DTF, pairwise ~ Treat)
summary(posthoc_DTF)## $emmeans
## Treat emmean SE df lower.CL upper.CL
## 0 58.5 1.37 22.1 55.6 61.3
## 15 59.8 1.37 22.1 56.9 62.6
## 45 61.1 1.37 22.1 58.2 63.9
## 75 62.7 1.37 22.1 59.8 65.5
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Treat0 - Treat15 -1.30 0.677 335 -1.921 0.2211
## Treat0 - Treat45 -2.59 0.677 335 -3.826 0.0009
## Treat0 - Treat75 -4.19 0.677 335 -6.190 <.0001
## Treat15 - Treat45 -1.29 0.677 335 -1.905 0.2280
## Treat15 - Treat75 -2.89 0.677 335 -4.269 0.0001
## Treat45 - Treat75 -1.60 0.677 335 -2.364 0.0860
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates#PHT
model_ML_PHT <- lmer(PHT ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
summary(model_ML_PHT)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PHT ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## REML criterion at convergence: 3164
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6208 -0.4605 -0.0225 0.4160 7.4174
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 1697.67 41.203
## Rep (Intercept) 39.37 6.275
## Residual 313.95 17.719
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 148.100 10.280 20.227 14.407 4.23e-12 ***
## Treat15 -9.624 2.641 335.000 -3.644 0.000311 ***
## Treat45 -9.649 2.641 335.000 -3.653 0.000301 ***
## Treat75 -16.542 2.641 335.000 -6.263 1.16e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.128
## Treat45 -0.128 0.500
## Treat75 -0.128 0.500 0.500plot(model_ML_PHT)
anova(model_ML_PHT)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treat 12482 4160.6 3 335 13.253 3.392e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VarCorr(model_ML_PHT)## Groups Name Std.Dev.
## Pedigree (Intercept) 41.2028
## Rep (Intercept) 6.2747
## Residual 17.7185posthoc_PHT <- emmeans(model_ML_PHT, pairwise ~ Treat)
summary(posthoc_PHT)## $emmeans
## Treat emmean SE df lower.CL upper.CL
## 0 148 10.3 20.2 127 170
## 15 138 10.3 20.2 117 160
## 45 138 10.3 20.2 117 160
## 75 132 10.3 20.2 110 153
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Treat0 - Treat15 9.6244 2.64 335 3.644 0.0018
## Treat0 - Treat45 9.6489 2.64 335 3.653 0.0017
## Treat0 - Treat75 16.5422 2.64 335 6.263 <.0001
## Treat15 - Treat45 0.0244 2.64 335 0.009 1.0000
## Treat15 - Treat75 6.9178 2.64 335 2.619 0.0453
## Treat45 - Treat75 6.8933 2.64 335 2.610 0.0465
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates#Seed Set
# Rescale the data scoring from 0–10 to 1–9 to standardize and to avoid 0 data. 1 bieng low seed set and 9 being full seed set, anything lessthan 1 converted to 1 and all greater than 9 converted to 9
TFMSA_NO_NA$SeedSet_round<- round(TFMSA_NO_NA$SeedSet)
TFMSA_NO_NA$SeedSet_rescaled <- ifelse(TFMSA_NO_NA$SeedSet_round < 1, 1,
ifelse(TFMSA_NO_NA$SeedSet_round > 9, 9,
TFMSA_NO_NA$SeedSet_round))
# Check the new range
#range(TFMSA_NO_NA$SeedSet_rescaled)
#summary(TFMSA_NO_NA$SeedSet_rescaled)
model_ML_SS <- lmer(SeedSet_rescaled ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
summary(model_ML_SS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SeedSet_rescaled ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## REML criterion at convergence: 1352
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5786 -0.3210 -0.0390 0.2426 4.3542
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 0.29019 0.5387
## Rep (Intercept) 0.09879 0.3143
## Residual 2.30102 1.5169
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.5556 0.2479 18.7606 22.41 5.31e-15 ***
## Treat15 -4.0333 0.2261 335.0000 -17.84 < 2e-16 ***
## Treat45 -4.4889 0.2261 335.0000 -19.85 < 2e-16 ***
## Treat75 -4.5556 0.2261 335.0000 -20.15 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.456
## Treat45 -0.456 0.500
## Treat75 -0.456 0.500 0.500plot(model_ML_SS)
anova(model_ML_SS)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treat 1297.2 432.42 3 335 187.92 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VarCorr(model_ML_SS)## Groups Name Std.Dev.
## Pedigree (Intercept) 0.53869
## Rep (Intercept) 0.31431
## Residual 1.51691posthoc_SS <- emmeans(model_ML_SS, pairwise ~ Treat)
summary(posthoc_SS)## $emmeans
## Treat emmean SE df lower.CL upper.CL
## 0 5.56 0.248 18.8 5.036 6.07
## 15 1.52 0.248 18.8 1.003 2.04
## 45 1.07 0.248 18.8 0.547 1.59
## 75 1.00 0.248 18.8 0.481 1.52
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Treat0 - Treat15 4.0333 0.226 335 17.837 <.0001
## Treat0 - Treat45 4.4889 0.226 335 19.851 <.0001
## Treat0 - Treat75 4.5556 0.226 335 20.146 <.0001
## Treat15 - Treat45 0.4556 0.226 335 2.015 0.1846
## Treat15 - Treat75 0.5222 0.226 335 2.309 0.0979
## Treat45 - Treat75 0.0667 0.226 335 0.295 0.9911
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates#Germination
model_ML_Germ <- lmer(Germ ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_Germ)
summary(model_ML_Germ)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Germ ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_Germ
##
## REML criterion at convergence: 2559.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2181 -0.4156 0.0804 0.5911 2.5888
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 147.02 12.125
## Rep (Intercept) 18.49 4.301
## Residual 165.09 12.849
## Number of obs: 317, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 83.798 3.704 20.706 22.621 4.44e-16 ***
## Treat15 -6.366 1.981 291.885 -3.213 0.001459 **
## Treat45 -7.399 2.019 292.132 -3.665 0.000294 ***
## Treat75 -10.016 2.079 293.337 -4.817 2.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.253
## Treat45 -0.248 0.470
## Treat75 -0.241 0.464 0.458plot(model_ML_Germ)
anova(model_ML_Germ)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treat 4334.9 1445 3 292.25 8.7525 1.411e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VarCorr(model_ML_Germ)## Groups Name Std.Dev.
## Pedigree (Intercept) 12.1253
## Rep (Intercept) 4.3005
## Residual 12.8488posthoc_Germ <- emmeans(model_ML_Germ, pairwise ~ Treat)
summary(posthoc_Germ)## $emmeans
## Treat emmean SE df lower.CL upper.CL
## 0 83.8 3.70 21.1 76.1 91.5
## 15 77.4 3.73 21.7 69.7 85.2
## 45 76.4 3.75 22.2 68.6 84.2
## 75 73.8 3.79 22.9 65.9 81.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Treat0 - Treat15 6.37 1.98 292 3.213 0.0079
## Treat0 - Treat45 7.40 2.02 293 3.664 0.0017
## Treat0 - Treat75 10.02 2.08 294 4.815 <.0001
## Treat15 - Treat45 1.03 2.06 292 0.502 0.9586
## Treat15 - Treat75 3.65 2.10 293 1.735 0.3075
## Treat45 - Treat75 2.62 2.13 293 1.226 0.6111
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates#Diagnosis: The plots on linear mixed model suggest the data strongly violates key assumptions #(e.g., normality and homoscedasticity) so we want to be careful on p and used robust method that accounts outlayers and coorlation of variances
#robust linear mixed model -rlmer
#DTF
model_robust_DTF <- rlmer(DTF ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
summary(model_robust_DTF)## Robust linear mixed model fit by DAStau
## Formula: DTF ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8448 -0.5773 -0.0726 0.7213 5.1691
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 25.641 5.064
## Rep (Intercept) 2.035 1.426
## Residual 15.509 3.938
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 58.2702 1.4518 40.14
## Treat15 1.5264 0.6021 2.54
## Treat45 2.7532 0.6021 4.57
## Treat75 3.8785 0.6021 6.44
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.207
## Treat45 -0.207 0.500
## Treat75 -0.207 0.500 0.500
##
## Robustness weights for the residuals:
## 290 weights are ~= 1. The remaining 70 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.260 0.626 0.829 0.753 0.943 0.998
##
## Robustness weights for the random effects:
## 21 weights are ~= 1. The remaining 2 ones are
## 6 17
## 0.782 0.781
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 1 (Pedigree):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 2 (Rep):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)plot(model_robust_DTF)
posthoc_robust_DTF <- emmeans(model_robust_DTF, pairwise ~ Treat)
summary(posthoc_robust_DTF)## $emmeans
## Treat emmean SE df asymp.LCL asymp.UCL
## 0 58.3 1.45 Inf 55.4 61.1
## 15 59.8 1.45 Inf 57.0 62.6
## 45 61.0 1.45 Inf 58.2 63.9
## 75 62.1 1.45 Inf 59.3 65.0
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Treat0 - Treat15 -1.53 0.602 Inf -2.535 0.0547
## Treat0 - Treat45 -2.75 0.602 Inf -4.573 <.0001
## Treat0 - Treat75 -3.88 0.602 Inf -6.441 <.0001
## Treat15 - Treat45 -1.23 0.602 Inf -2.038 0.1742
## Treat15 - Treat75 -2.35 0.602 Inf -3.906 0.0005
## Treat45 - Treat75 -1.13 0.602 Inf -1.869 0.2415
##
## P value adjustment: tukey method for comparing a family of 4 estimates#PHT
model_robust_PHT <- rlmer(PHT ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
summary(model_robust_PHT)## Robust linear mixed model fit by DAStau
## Formula: PHT ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.7879 -0.5857 0.0055 0.5693 12.5177
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 1997.92 44.698
## Rep (Intercept) 19.88 4.458
## Residual 129.94 11.399
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 144.446 11.066 13.053
## Treat15 -7.171 1.743 -4.115
## Treat45 -9.037 1.743 -5.185
## Treat75 -14.458 1.743 -8.296
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.079
## Treat45 -0.079 0.500
## Treat75 -0.079 0.500 0.500
##
## Robustness weights for the residuals:
## 286 weights are ~= 1. The remaining 74 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.107 0.535 0.714 0.680 0.902 0.999
##
## Robustness weights for the random effects:
## 18 weights are ~= 1. The remaining 5 ones are
## 4 7 11 15 19
## 0.998 0.684 0.854 0.897 0.771
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 1 (Pedigree):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 2 (Rep):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)plot(model_robust_PHT)
posthoc_robust_PHT <- emmeans(model_robust_PHT, pairwise ~ Treat)
summary(posthoc_robust_PHT)## $emmeans
## Treat emmean SE df asymp.LCL asymp.UCL
## 0 144 11.1 Inf 123 166
## 15 137 11.1 Inf 116 159
## 45 135 11.1 Inf 114 157
## 75 130 11.1 Inf 108 152
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Treat0 - Treat15 7.17 1.74 Inf 4.115 0.0002
## Treat0 - Treat45 9.04 1.74 Inf 5.185 <.0001
## Treat0 - Treat75 14.46 1.74 Inf 8.296 <.0001
## Treat15 - Treat45 1.87 1.74 Inf 1.071 0.7075
## Treat15 - Treat75 7.29 1.74 Inf 4.181 0.0002
## Treat45 - Treat75 5.42 1.74 Inf 3.111 0.0101
##
## P value adjustment: tukey method for comparing a family of 4 estimates#Seed set
model_robust_rescaledSS <- rlmer(SeedSet_rescaled ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA)
summary(model_robust_rescaledSS)## Robust linear mixed model fit by DAStau
## Formula: SeedSet_rescaled ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_NO_NA
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3509.0 -0.3 0.0 0.2 5611.2
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 1.427e-07 0.0003778
## Rep (Intercept) 8.736e-08 0.0002956
## Residual 2.032e-06 0.0014255
## Number of obs: 360, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.0009896 0.0002247 26712
## Treat15 -5.0006671 0.0002179 -22945
## Treat45 -5.0009504 0.0002179 -22946
## Treat75 -5.0010333 0.0002179 -22946
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.485
## Treat45 -0.485 0.500
## Treat75 -0.485 0.500 0.500
##
## Robustness weights for the residuals:
## 3 observations c(126,129,209) are outliers with |weight| <= 0.00027 ( < 0.00028);
## 266 weights are ~= 1. The remaining 91 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000320 0.000639 0.000958 0.000979 0.000960 0.001920
##
## Robustness weights for the random effects:
## 19 weights are ~= 1. The remaining 4 ones are
## 7 11 18 19
## 0.373 0.719 0.997 0.987
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 1 (Pedigree):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 2 (Rep):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)plot(model_robust_rescaledSS)
posthoc_robust_rescaledSS <- emmeans(model_robust_rescaledSS, pairwise ~ Treat)
summary(posthoc_robust_rescaledSS)## $emmeans
## Treat emmean SE df asymp.LCL asymp.UCL
## 0 6.00099 0.00022465 Inf 6.00055 6.0014
## 15 1.00032 0.00022465 Inf 0.99988 1.0008
## 45 1.00004 0.00022465 Inf 0.99960 1.0005
## 75 0.99996 0.00022465 Inf 0.99952 1.0004
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Treat0 - Treat15 5.000667 0.00021795 Inf 22944.614 <.0001
## Treat0 - Treat45 5.000950 0.00021795 Inf 22945.914 <.0001
## Treat0 - Treat75 5.001033 0.00021795 Inf 22946.294 <.0001
## Treat15 - Treat45 0.000283 0.00021795 Inf 1.300 0.5631
## Treat15 - Treat75 0.000366 0.00021795 Inf 1.680 0.3343
## Treat45 - Treat75 0.000083 0.00021795 Inf 0.380 0.9813
##
## P value adjustment: tukey method for comparing a family of 4 estimates#Germination
model_robust_Germ <- rlmer(Germ ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_Germ)
summary(model_robust_Germ)## Robust linear mixed model fit by DAStau
## Formula: Germ ~ Treat + (1 | Pedigree) + (1 | Rep)
## Data: TFMSA_Germ
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4568 -0.5955 -0.0015 0.6319 3.0447
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pedigree (Intercept) 146.58 12.11
## Rep (Intercept) 14.52 3.81
## Residual 112.04 10.58
## Number of obs: 317, groups: Pedigree, 18; Rep, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 84.850 3.598 23.582
## Treat15 -5.338 1.674 -3.188
## Treat45 -6.970 1.706 -4.085
## Treat75 -8.593 1.758 -4.889
##
## Correlation of Fixed Effects:
## (Intr) Tret15 Tret45
## Treat15 -0.220
## Treat45 -0.216 0.470
## Treat75 -0.210 0.464 0.458
##
## Robustness weights for the residuals:
## 255 weights are ~= 1. The remaining 62 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.246 0.526 0.734 0.707 0.929 0.996
##
## Robustness weights for the random effects:
## 19 weights are ~= 1. The remaining 4 ones are
## 2 13 18 19
## 0.997 0.727 0.743 0.905
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 1 (Pedigree):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 2 (Rep):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)plot(model_robust_Germ)
posthoc_robust_Germ <- emmeans(model_robust_Germ, pairwise ~ Treat)
summary(posthoc_robust_Germ)## $emmeans
## Treat emmean SE df asymp.LCL asymp.UCL
## 0 84.8 3.60 Inf 77.8 91.9
## 15 79.5 3.62 Inf 72.4 86.6
## 45 77.9 3.63 Inf 70.8 85.0
## 75 76.3 3.66 Inf 69.1 83.4
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Treat0 - Treat15 5.34 1.67 Inf 3.188 0.0078
## Treat0 - Treat45 6.97 1.71 Inf 4.085 0.0003
## Treat0 - Treat75 8.59 1.76 Inf 4.889 <.0001
## Treat15 - Treat45 1.63 1.74 Inf 0.938 0.7845
## Treat15 - Treat75 3.26 1.78 Inf 1.831 0.2585
## Treat45 - Treat75 1.62 1.80 Inf 0.900 0.8047
##
## P value adjustment: tukey method for comparing a family of 4 estimates
# Create a list of models
models <- list(
SeedSet = rlmer(SeedSet_rescaled ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA),
PHT = rlmer(PHT ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA),
DTF = rlmer(DTF ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA),
Germ = rlmer(Germ ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_Germ)
)
# Initialize an empty data frame for storing results
results_table <- data.frame(
Trait = character(),
Effect = character(),
Estimate = numeric(),
Std.Error = numeric(),
t.value = numeric(),
p.value = character(),
Random_Effect_Variance = numeric(),
Residual_Variance = numeric(),
stringsAsFactors = FALSE
)
# Function to extract model results
extract_model_results <- function(model, trait) {
# Fixed effects
fixed <- summary(model)$coefficients
random <- summary(model)$varcor
# Extract random effect variance (ensure correct handling of structure)
random_var_Pedigree <- as.numeric(random$Pedigree[1]) # Random intercept variance
random_var_Rep <- as.numeric(random$Rep[1])
residual_var <- as.numeric(attr(random, "sc")^2) # Residual variance
# Extract fixed effects
fixed_results <- data.frame(
Trait = trait,
Effect = rownames(fixed),
Estimate = fixed[, "Estimate"],
Std.Error = fixed[, "Std. Error"],
t.value = fixed[, "t value"],
p.value = ifelse(abs(fixed[, "t value"]) > 2, "<0.05", ">=0.05"), # Approximate significance
Random_Effect_Variance_Ped = random_var_Pedigree,
Random_Effect_Variance_Rep = random_var_Rep,
Residual_Variance = residual_var
)
return(fixed_results)
}
# Loop through the models and extract results
for (trait in names(models)) {
model <- models[[trait]]
trait_results <- extract_model_results(model, trait)
results_table <- rbind(results_table, trait_results)
}
# Save the results as a CSV file
write.csv(results_table, "robust Mixed_Model_Results_2023.csv", row.names = FALSE)
# Print the results
results_table
models <- list(
SeedSet = lmer(SeedSet_rescaled ~ Treat + (1 | Pedigree)+ (1 | Rep), data = TFMSA_NO_NA),
PHT = lmer(PHT ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_NO_NA),
DTF = lmer(DTF ~ Treat + (1 | Pedigree)+ (1 | Rep), data = TFMSA_NO_NA),
Germ = rlmer(Germ ~ Treat + (1 | Pedigree) + (1 | Rep), data = TFMSA_Germ)
)
# Initialize an empty data frame for storing results
results_table <- data.frame(
Trait = character(),
Effect = character(),
Estimate = numeric(),
Std.Error = numeric(),
t.value = numeric(),
p.value = character(),
Random_Effect_Variance = numeric(),
Residual_Variance = numeric(),
stringsAsFactors = FALSE
)
# Function to extract model results
extract_model_results <- function(model, trait) {
# Fixed effects
fixed <- summary(model)$coefficients
random <- summary(model)$varcor
# Extract random effect variance (ensure correct handling of structure)
random_var_Pedigree <- as.numeric(random$Pedigree[1]) # Random intercept variance
random_var_Rep <- as.numeric(random$Rep[1])
residual_var <- as.numeric(attr(random, "sc")^2) # Residual variance
# Extract fixed effects
fixed_results <- data.frame(
Trait = trait,
Effect = rownames(fixed),
Estimate = fixed[, "Estimate"],
Std.Error = fixed[, "Std. Error"],
t.value = fixed[, "t value"],
p.value = ifelse(abs(fixed[, "t value"]) > 2, "<0.05", ">=0.05"), # Approximate significance
Random_Effect_Variance_Ped = random_var_Pedigree,
Random_Effect_Variance_Rep = random_var_Rep,
Residual_Variance = residual_var
)
return(fixed_results)
}
# Loop through the models and extract results
for (trait in names(models)) {
model <- models[[trait]]
trait_results <- extract_model_results(model, trait)
results_table <- rbind(results_table, trait_results)
}
# Save the results as a CSV file
write.csv(results_table, "linear_Mixed_Model_Results_2023.csv", row.names = FALSE)
# Print the results
results_table