基于Hat Matrix手动计算Beta Regression遇阻，请求技术协助

我尝试参考逻辑回归帽子矩阵的思路手动计算贝塔回归的帽子矩阵，但代码无法运行，也不确定方法是否正确，恳请协助。

原代码：

require(betareg)
df<-data("ReadingSkills")
y<-ReadingSkills$accuracy
n<-length(y)
x1<-rnorm(n,0,1)
x2<-rnorm(n,0,1)

X<-cbind(1,x1, x2)


bfit <- solve(t(X) %*% X + u * diag(3)) %*% t(X) %*% y
bfit1<-betareg(accuracy ~ x1+x2, data = ReadingSkills, link="logit")
bfit
bfit1

v <- 1/(1+exp(-X%*%bfit))
VX <- X*v 
H <- VX%*%solve(crossprod(VX,VX),t(VX))

问题分析 #

1. 未定义变量u：代码中直接使用了u但未赋值，这是导致运行报错的直接原因。
2. 模型逻辑混淆：贝塔回归和逻辑回归的参数估计、帽子矩阵计算逻辑完全不同，不能直接照搬逻辑回归的方式。贝塔回归针对(0,1)区间的响应变量，采用极大似然迭代估计参数，包含均值和精度两个部分，不是简单的线性回归形式。

修正方案与代码 #

贝塔回归的帽子矩阵需要结合模型的均值链接函数、精度参数来计算，步骤如下：

# 加载包与数据
library(betareg)
data("ReadingSkills")
y <- ReadingSkills$accuracy
n <- length(y)
# 生成自变量（固定随机种子保证结果可复现）
set.seed(123)
x1 <- rnorm(n, 0, 1)
x2 <- rnorm(n, 0, 1)
X <- cbind(1, x1, x2)

# 1. 先拟合正确的贝塔回归模型
bfit1 <- betareg(accuracy ~ x1 + x2, data = ReadingSkills, link = "logit")
# 提取均值部分的参数估计
beta_hat <- coef(bfit1)[1:3]
# 提取精度参数（phi）
phi_hat <- coef(bfit1)[4]

# 2. 计算线性预测值与均值（logit链接的逆变换）
eta <- X %*% beta_hat
mu <- plogis(eta)

# 3. 计算权重矩阵：结合链接函数导数与贝塔分布方差
W <- diag(n)
for (i in 1:n) {
  # 贝塔分布的方差
  var_mu <- mu[i] * (1 - mu[i]) / (1 + phi_hat)
  # logit链接函数的导数平方
  link_deriv_sq <- (mu[i] * (1 - mu[i]))^2
  # 权重计算
  W[i,i] <- phi_hat / var_mu * link_deriv_sq
}

# 4. 构造加权设计矩阵并计算帽子矩阵
X_weighted <- sqrt(W) %*% X
H <- X_weighted %*% solve(crossprod(X_weighted)) %*% t(X_weighted)

# 查看帽子矩阵维度（应为n×n）
dim(H)

说明 #

• 贝塔回归的帽子矩阵核心是加权最小二乘形式的设计矩阵，权重同时考虑了链接函数的导数特性和响应变量的方差结构。
• 由于贝塔回归参数没有解析解，必须先通过betareg拟合得到参数估计，才能手动推导帽子矩阵。

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