在R的ProfileLikelihood包中为GLM模型添加偏移量的正确方法

背景 #

我在做Agresti《分类数据分析》的习题4.10，题目给出的表格行对应年龄组、列对应心脏瓣膜类型，单元格内是死亡数（上方）和总时长（下方）。要求使用的模型为：
$$
\log(\mu_{ij}) = \alpha + 1 \times \log(t_{ij}) + \beta_1 a_i + \beta_2 \nu_j
$$
其中$t_{ij}$是单元格总时长，$a_i$为“55+”指示变量，$\nu_j$为“Mitral”指示变量，目标是求解$\beta_1$的95%轮廓似然置信区间。

遇到的问题 #

使用ProfileLikelihood包的profilelike.glm函数时，多次尝试添加偏移量log(tot_time)均报错：

• 传入offset.glm = log(tot_time)，报错“object 'tot_time' not found”
• 传入offset.glm = "log(tot_time)"，报错“invalid type (list) for variable 'offset(glm.off)'”
• 在公式中使用offset.glm(log(tot_time))，报错“could not find function "offset.glm"”
• 在公式中使用offset(log(tot_time))可以运行，但不符合文档“应使用offset.glm而非offset”的要求

正确解决方法 #

根据ProfileLikelihood包的文档，profilelike.glm的offset.glm参数需要传入预计算完成的数值向量，而非表达式或字符串。具体步骤如下：

1. 提前在数据框中计算好偏移量列：

# 假设你的数据框名为df，总时长列名为tot_time
df$offset_col <- log(df$tot_time)

2. 调用profilelike.glm时，通过offset.glm参数传入这个预计算的向量，同时公式中不再使用offset()函数：

# 假设死亡数列是death，age55_plus是a_i对应的变量，mitral是ν_j对应的变量
library(ProfileLikelihood)

# 拟合模型并计算β1的轮廓似然置信区间
pl_output <- profilelike.glm(
  formula = death ~ age55_plus + mitral,
  family = poisson(),
  data = df,
  offset.glm = df$offset_col,
  which.par = "age55_plus"  # 指定要计算置信区间的参数对应变量
)

# 查看95%轮廓似然置信区间
confint(pl_output)

关键说明 #

文档中“应使用offset.glm而非offset”的意思是：不要在公式里用offset()函数声明偏移量，而是通过profilelike.glm的专属参数offset.glm传入已经计算好的数值向量。之前的错误都是因为参数类型不匹配——要么直接传未计算的表达式，要么传字符串，都不符合参数要求。

内容的提问来源于stack exchange，提问作者William M.