【问题标题】:Extract estimates of GAM提取 GAM 的估计值
【发布时间】:2012-10-17 22:17:31
【问题描述】:

我对 R 相当陌生,目前正在阅读一本书“Generalized Additive Models”,Wood 的 R 简介(2006 年)并进行了一些练习,特别是关于空气污染和死亡的部分,这是我的领域兴趣。我使用 mgcv 包运行以下模型。

library(gamair) 
library(mgcv) 
data(chicago) 

ap1<-gam(death ~ pm10median + so2median + o3median +s(time,bs="cr",k=200)+ s(tmpd,bs="cr"), data=chicago,family=poisson)

如何提取 pm10median 和 x 的 95% CI 的效果估计值并将输出导出为 CSV 或任何其他选项?

【问题讨论】:

    标签: r gam mgcv


    【解决方案1】:

    保存模型的摘要

    summary_model <- summary(ap1)
    

    您想要的部分(对于线性项)在 p.table 元素中

    summary_model$p.table
                    Estimate   Std. Error      z value    Pr(>|z|)
    (Intercept) 4.7457425965 1.480523e-03 3205.4510971 0.000000000
    pm10median  0.0002551498 9.384003e-05    2.7189871 0.006548217
    so2median   0.0008898646 5.543272e-04    1.6053056 0.108426561
    o3median    0.0002212612 2.248015e-04    0.9842516 0.324991826
    
    
    write.csv(summary_model$p.table, file = 'p_table.csv')
    

    如果你想要样条项,那么这是

    summary_model$s.table
                   edf     Ref.df    Chi.sq       p-value
    s(time) 167.327973 187.143378 1788.8201 4.948832e-259
    s(tmpd)   8.337121   8.875807  110.5231  1.412415e-19
    

    您可以手动计算 95% CI 并添加这些如果您愿意。 (由于高 DF 将使用 Z 分数)

    p_table <- data.frame(summary_model$p.table)
    p_table <- within(p_table, {lci <- Estimate - qnorm(0.975) * Std..Error
                                uci <- Estimate + qnorm(0.975) * Std..Error})
    p_table
                   Estimate   Std..Error      z.value    Pr...z..          uci           lci
    (Intercept) 4.7457425965 1.480523e-03 3205.4510971 0.000000000 4.7486443674  4.742841e+00
    pm10median  0.0002551498 9.384003e-05    2.7189871 0.006548217 0.0004390729  7.122675e-05
    so2median   0.0008898646 5.543272e-04    1.6053056 0.108426561 0.0019763260 -1.965968e-04
    o3median    0.0002212612 2.248015e-04    0.9842516 0.324991826 0.0006618641 -2.193416e-04\
    

    根据 cmets 编辑

    如果您有多个 gam 模型,例如 ap1ap2ap3,并且您希望系统地处理它们,Rish 的方法是将它们放在一个列表中并使用 @987654330 @

    # create list
    model_list <- list(ap1, ap2, ap3)
    # give the elements useful names
    names(model_list) <- c('ap1','ap2','ap3')
    
    # get the summaries using `lapply
    
    summary_list <- lapply(model_list, summary)
    
    # extract the coefficients from these summaries
    
     p.table_list <- lapply(summary_list, `[[`, 'p.table')
    
     s.table_list <- lapply(summary_list, `[[`, 's.table')
    

    您现在创建的列表是相关组件。

    【讨论】:

    • 非常感谢您的有用回答,这正是我一直在寻找的。​​span>
    • 我想重新提出这个问题,希望我不要违反论坛规则。如何系统地存储几个 gam 对象(mgcv 包)的系数?目前我正在单独存储每个估计值,无法弄清楚如何为多个对象执行此操作。
    • 尊敬的 mnel,谢谢,非常感谢您的帮助。系数可以在汇总列表中找到,我稍作修改如下:summary_list [[, 'p.table')
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2018-02-26
    • 1970-01-01
    • 1970-01-01
    • 2021-02-26
    相关资源
    最近更新 更多