You need to enable JavaScript to run this app.
最新活动
大模型
产品
解决方案
定价
生态与合作
支持与服务
开发者
了解我们

R语言中含Formula参数的Fisher算法GLM函数调用报错排查

解决GLM Fisher算法函数中的invalid type (closure) for variable 't'错误

嗨,我来帮你搞定这个报错!你遇到的invalid type (closure) for variable 't'错误,核心原因是数据里的自变量叫t,而R自带的转置函数也叫t()——当model.frame处理公式时,它把t当成了内置函数(类型是closure),而不是你数据里的变量,这就导致了报错。下面分步骤给你解决:

一、先搞定变量名冲突问题

有两种简单的方式解决这个冲突:

方式1:重命名数据里的t

直接把数据里的t改成一个不会和内置函数重名的名字,比如time

# 加载数据后修改列名
load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata"))
# base R写法,不需要额外包
data$time <- data$t
data$t <- NULL
# 更新公式
f <- formula(y ~ time)
# 调用函数
myglm(f, data)

方式2:在公式里用反引号包裹t

如果不想改数据,就在公式里用反引号强制R把t识别成变量:

load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata"))
# 用反引号把t包起来
f <- formula(y ~ `t`)
myglm(f, data)

二、优化你的Fisher算法函数(修复其他潜在问题)

除了变量名冲突,你的函数还有几个可能导致错误或结果不准的地方,我帮你做了优化:

library(stats) # 确保LogLik函数可用

myglm <- function(formula, data, start = NULL, epsilon = 0.01) {
  # 用model.frame正确提取数据,避免依赖列顺序
  mf <- model.frame(formula, data = data, drop.unused.levels = TRUE)
  X <- model.matrix(formula, mf)
  Y <- model.response(mf)
  
  n <- nrow(X)
  p <- ncol(X)
  
  # 优化beta初始化逻辑:优先用传入的start,否则用最小二乘估计
  if (is.null(start)) {
    beta <- solve(crossprod(X), crossprod(X, Y)) # 最小二乘初始化更合理
  } else {
    if (length(start) != p) stop("start的长度必须和参数个数一致")
    beta <- start
  }
  
  beta_0 <- beta + 10*epsilon # 确保第一次循环能正常执行
  # Fisher迭代过程
  while (norm(beta - beta_0, type = "2") / norm(beta_0, type = "2") > epsilon) {
    beta_0 <- beta
    eta <- X %*% beta
    lambda <- exp(eta)
    # 用crossprod替代t(X)%*%X,计算更高效
    F_mat <- crossprod(X, X * lambda)
    s <- crossprod(X, Y - lambda) # 得分函数
    beta <- beta + solve(F_mat) %*% s
  }
  
  vcov <- solve(F_mat)
  # 动态生成系数表行名,适配不同参数数量
  coef_table <- cbind(beta, sqrt(diag(vcov)))
  colnames(coef_table) <- c("Coefficients", "Standard error")
  rownames(coef_table) <- colnames(X)
  
  # 计算偏差:修正对数似然计算方式
  mod_sat <- glm(formula, family = poisson(link = "log"), data = mf)
  log_likelihood <- sum(Y * eta - lambda)
  deviance <- 2*(as.numeric(LogLik(mod_sat)) - log_likelihood)
  
  return(list(coef = coef_table, deviance = deviance, vcov = vcov))
}

优化点说明:

  • model.framemodel.response提取Y,不再依赖“数据第一列是响应变量”的假设,更鲁棒;
  • 优化beta初始化逻辑,支持传入自定义初始值,否则用最小二乘估计(比全1向量更合理);
  • crossprod替代t(X)%*%X,计算更高效且代码更简洁;
  • 动态生成系数表的行名,避免当参数个数不是2时出现行名不匹配的错误;
  • 修正了对数似然的计算方式(用sum而非矩阵乘法,避免维度问题);
  • epsilon增加了默认参数,方便后续调整迭代精度。

三、测试修正后的代码

现在用上面的任意一种变量名处理方式,加上优化后的函数,就能正常运行了:

# 示例:用反引号的方式
load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata"))
f <- formula(y ~ `t`)
result <- myglm(f, data)
# 查看系数表
print(result$coef)
# 查看偏差
print(result$deviance)

内容的提问来源于stack exchange,提问作者user9231414

火山引擎 最新活动