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.frame和model.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




