首页 > 解决方案 > for 循环中的多个 GLM

问题描述

我正在寻找遵循 NB 分布的一大组 (274) 相关响应变量的参数估计。目标是使用具有两级分类预测器的每个变量的广义线性模型的参数。我熟悉 R 中的 MASS 包,但似乎没有办法进行多变量分析。有没有一种方法可以在 R 中实现(例如 for 循环)。

提前感谢您的任何帮助/建议

标签: rglmmultivariate-testing

解决方案


你首先将你的模型构造成一个字符串,

options = c('Sex', 'Sex/(Age + Eth*Lrn)')
evaluation_str = c()
for(i in c(1:2)){
        evaluation_str <- c(evaluation_str, paste0('glm.nb(Days ~ ', options[i], ', data = quine)'))
}

> evaluation_str
[1] "glm.nb(Days ~ Sex, data = quine)"                 "glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)"

然后使用eval()和 parse() 函数运行一个循环,我更喜欢lapply(),所以我不需要先启动一个向量作为结果容器。

results <- lapply(c(1:2), function(x){
        eval(parse(text=evaluation_str[i]))
})

> results
[[1]]

Call:  glm.nb(formula = Days ~ Sex/(Age + Eth * Lrn), data = quine, 
    init.theta = 1.597990733, link = log)

Coefficients:
    (Intercept)             SexM       SexF:AgeF1       SexM:AgeF1       SexF:AgeF2       SexM:AgeF2       SexF:AgeF3  
        3.01919         -0.47541         -0.70887         -0.72373         -0.61486          0.62842         -0.34235  
     SexM:AgeF3        SexF:EthN        SexM:EthN       SexF:LrnSL       SexM:LrnSL  SexF:EthN:LrnSL  SexM:EthN:LrnSL  
        1.15084         -0.07312         -0.67899          0.94358          0.23891         -1.35849          0.76142  

Degrees of Freedom: 145 Total (i.e. Null);  132 Residual
Null Deviance:      234.6 
Residual Deviance: 167.6    AIC: 1093

[[2]]

Call:  glm.nb(formula = Days ~ Sex/(Age + Eth * Lrn), data = quine, 
    init.theta = 1.597990733, link = log)

Coefficients:
    (Intercept)             SexM       SexF:AgeF1       SexM:AgeF1       SexF:AgeF2       SexM:AgeF2       SexF:AgeF3  
        3.01919         -0.47541         -0.70887         -0.72373         -0.61486          0.62842         -0.34235  
     SexM:AgeF3        SexF:EthN        SexM:EthN       SexF:LrnSL       SexM:LrnSL  SexF:EthN:LrnSL  SexM:EthN:LrnSL  
        1.15084         -0.07312         -0.67899          0.94358          0.23891         -1.35849          0.76142  

Degrees of Freedom: 145 Total (i.e. Null);  132 Residual
Null Deviance:      234.6 
Residual Deviance: 167.6    AIC: 1093

如果要检查每个单独的模型。你可以这样做:

> coef <- results[[1]]$coefficients
> coef
    (Intercept)            SexM      SexF:AgeF1      SexM:AgeF1      SexF:AgeF2      SexM:AgeF2      SexF:AgeF3 
      3.0191882      -0.4754051      -0.7088737      -0.7237349      -0.6148641       0.6284201      -0.3423469 
     SexM:AgeF3       SexF:EthN       SexM:EthN      SexF:LrnSL      SexM:LrnSL SexF:EthN:LrnSL SexM:EthN:LrnSL 
      1.1508429      -0.0731245      -0.6789869       0.9435760       0.2389097      -1.3584906       0.7614156 

推荐阅读