代码之家  ›  专栏  ›  技术社区  ›  Jake

更新函数中的模型因子级别时工作不正常

  •  1
  • Jake  · 技术社区  · 6 年前

    我在R中运行了一些GLM模型,这些模型与我正在进行的喂养试验有关。我用两个预测因子回归我感兴趣的变量:一个因素有三个水平,一个连续变量。我想将每一级的拦截因素相互比较,以确定它们是否不同。为此,我编写了一个函数(在下面的可复制代码中称为interceptcarbire),它释放因子并更新模型,然后保存每个模型的结果。这是我对截获的所有成对比较的快速方法。

    问题是,当我运行函数时,它似乎没有正确地更新模型。 返回的列表中的每个项都是相同的,当它们应该改变时,每个项都有一个不同的因子级别,与其他级别正在比较的“(截距)”不同。我怀疑这与函数的环境有关,但我不确定。我还没有在StackOverflow或Google上找到类似的例子。

    下面是一个可重复的例子:

    food <- as.factor(rep(c("a", "b", "c"), each = 20))
    variable <- rbinom(60, 1, 0.7)
    movement <- rgamma(60, 10, 2)
    binomial.model <- glm(variable ~ food,
             family = "binomial")
    gamma.model <- glm(movement ~ food,
                   family = Gamma)
    interceptCompare <- function(model, factor) {
      results <- list() # empty list to store results
      for (i in unique(factor)) {
        factor <- relevel(factor, ref = i)
        model <- update(model)
        results[[i]] <- summary(model)$coefficients[1:3, ]
      }
      results <- lapply(results, function(x) round(x, 4))
      return(results)
    }
    
    interceptCompare(binomial.model, food)
    interceptCompare(gamma.model, food)
    
    2 回复  |  直到 6 年前
        1
  •  1
  •   Onyambu    6 年前

    您需要添加一行,以便更改数据,并在更新中使用它:

    interceptCompare <- function(model, factor) {
      results <- list() # empty list to store results
    
       s <- deparse(substitute(factor))#ADD THIS LINE
    
      for (i in unique(factor)) {
        factor <- relevel(factor, ref = i)
        model[["model"]][[s]] <- factor #CHANGE THE DATA IN THE MODEL
        model <- update(model,data=model[["model"]])# UPDATE THE MODEL
        results[[i]] <- summary(model)$coefficients[1:3, ]
      }
      results <- lapply(results, function(x) round(x, 4))
      return(results)
    }
    
    interceptCompare(binomial.model, food)
    $a
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)   1.3863     0.5590  2.4799   0.0131
    foodb        -0.7673     0.7296 -1.0516   0.2930
    foodc        -0.2877     0.7610 -0.3780   0.7054
    
    $b
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)   0.6190     0.4688  1.3205   0.1867
    fooda         0.7673     0.7296  1.0516   0.2930
    foodc         0.4796     0.6975  0.6876   0.4917
    
    $c
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)   1.0986     0.5164  2.1275   0.0334
    foodb        -0.4796     0.6975 -0.6876   0.4917
    fooda         0.2877     0.7610  0.3780   0.7054
    
    
    interceptCompare(gamma.model, food)
    $a
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)   0.2246     0.0156 14.3919   0.0000
    foodb        -0.0170     0.0213 -0.8022   0.4257
    foodc        -0.0057     0.0218 -0.2608   0.7952
    
    $b
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)   0.2076     0.0144 14.3919   0.0000
    fooda         0.0170     0.0213  0.8022   0.4257
    foodc         0.0114     0.0210  0.5421   0.5898
    
    $c
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)   0.2189     0.0152 14.3919   0.0000
    foodb        -0.0114     0.0210 -0.5421   0.5898
    fooda         0.0057     0.0218  0.2608   0.7952
    
        2
  •  0
  •   MrFlick    6 年前

    当你试图从你正在做的公式中替换符号时,你需要更加小心。你需要用R语言能理解的术语来表达它。你要以“食物”的名义传递,而不是像现在这样存储在“食物”向量中的值。这是一个更新,它似乎可以完成你想做的事情

    interceptCompare <- function(model, factor) {
      sym <- substitute(factor)
      results <- list() # empty list to store results
      for (i in unique(factor)) {
        change <- eval(bquote(~.-.(sym)+relevel(.(sym), ref=.(i))))
        new_model <- update(model, change)
        results[[i]] <- summary(new_model)$coefficients[1:3, ]
      }
      results <- lapply(results, function(x) round(x, 4))
      return(results)
    }
    

    在这里,我们用 substitute .然后我们使用 bquote() 建立一个新的公式,该公式将删除原始值,并使用特定的引用重新调整因子变量。然后我们把这个保存到一个新的对象中,这样我们就不会一直更新同一个模型。对于 binomial.model ,这个返回

    $`a`
                              Estimate Std. Error z value Pr(>|z|)
    (Intercept)                 0.8473     0.4879  1.7364   0.0825
    relevel(food, ref = "a")b   0.0000     0.6901  0.0000   1.0000
    relevel(food, ref = "a")c  -0.8473     0.6619 -1.2801   0.2005
    
    $b
                              Estimate Std. Error z value Pr(>|z|)
    (Intercept)                 0.8473     0.4879  1.7364   0.0825
    relevel(food, ref = "b")a   0.0000     0.6901  0.0000   1.0000
    relevel(food, ref = "b")c  -0.8473     0.6619 -1.2801   0.2005
    
    $c
                              Estimate Std. Error z value Pr(>|z|)
    (Intercept)                 0.0000     0.4472  0.0000   1.0000
    relevel(food, ref = "c")a   0.8473     0.6619  1.2801   0.2005
    relevel(food, ref = "c")b   0.8473     0.6619  1.2801   0.2005
    

    你可以看到它是如何改变 ref= 每次迭代时