如何从nlme包的系统发育控制PGLS ANOVA获取各变量效应量？

nlme包的anova.gls()确实不会直接输出平方和或效应量指标，以下是几种可行的解决方法：

1. 手动计算边际R² #

PGLS的R²需要基于加权平方和计算（因为广义最小二乘引入了系统发育相关权重），步骤如下：

# 提取模型响应变量与权重矩阵
y <- pgls.model$model[[all.vars(pgls.model$terms)[1]]]
V <- getVarCov(pgls.model)

# 计算加权总平方和
y_mean <- mean(y)
WSS_total <- t(y - y_mean) %*% solve(V) %*% (y - y_mean)

# 计算加权残差平方和
resid <- residuals(pgls.model, type = "response")
WSS_resid <- t(resid) %*% solve(V) %*% resid

# 计算边际R²（仅固定效应解释的变异比例）
pgls_r2 <- 1 - (WSS_resid / WSS_total)
cat("PGLS 边际R²:", round(as.numeric(pgls_r2), 3), "\n")

2. 使用MuMIn包快速计算R² #

MuMIn包的r.squaredGLMM()函数支持gls模型，可直接输出边际R²：

library(MuMIn)

# 计算R²
r2_results <- r.squaredGLMM(pgls.model)
print(r2_results)

输出中的R2m为边际R²（仅固定效应贡献），由于gls模型无随机效应，R2c会与R2m一致。

3. 计算Eta²与Cohen's f #

这两个指标需要通过嵌套模型对比，计算单个变量的加权平方和贡献：

# 复用方法1中的总加权平方和WSS_total

# 拟合不含Lnalt的嵌套模型
pgls_no_lnalt <- gls(Lnote ~ wing, data = Barbet, correlation = corPagel(value = 1, form = ~ species, phy = Tree))
# Lnalt的加权平方和
WSS_lnalt <- t(residuals(pgls_no_lnalt)) %*% solve(getVarCov(pgls_no_lnalt)) %*% residuals(pgls_no_lnalt) - WSS_resid
# Lnalt的Eta²与Cohen's f
eta2_lnalt <- WSS_lnalt / WSS_total
f_lnalt <- sqrt(eta2_lnalt / (1 - eta2_lnalt))

# 同理计算wing的效应量
pgls_no_wing <- gls(Lnote ~ Lnalt, data = Barbet, correlation = corPagel(value = 1, form = ~ species, phy = Tree))
WSS_wing <- t(residuals(pgls_no_wing)) %*% solve(getVarCov(pgls_no_wing)) %*% residuals(pgls_no_wing) - WSS_resid
eta2_wing <- WSS_wing / WSS_total
f_wing <- sqrt(eta2_wing / (1 - eta2_wing))

# 输出结果
cat("Lnalt的Eta²:", round(as.numeric(eta2_lnalt), 3), "\n")
cat("Lnalt的Cohen's f:", round(as.numeric(f_lnalt), 3), "\n")
cat("wing的Eta²:", round(as.numeric(eta2_wing), 3), "\n")
cat("wing的Cohen's f:", round(as.numeric(f_wing), 3), "\n")

内容的提问来源于stack exchange，提问作者User777