【问题标题】:How to get a data.frame with cases from a contingency table in r?如何从r中的列联表中获取带有案例的data.frame?
【发布时间】:2016-07-29 16:49:56
【问题描述】:

我想从书中重现一些计算(logit 回归)。这本书给出了一个列联表和结果。

这是表格:

    .


    example <- matrix(c(21,22,6,51), nrow = 2, byrow = TRUE)
    #Labels:
    rownames(example) <- c("Present","Absent")
    colnames(example) <- c(">= 55", "<55")

它给了我这个:

            >= 55 <55
    Present    21  22
    Absent      6  51

但要使用 glm() 函数,数据必须采用以下方式:

(两列,一列“Age”,一列“Present”,用0/1填充)

    age <- c(rep(c(0),27), rep(c(1),73))
    present <- c(rep(c(0),21), rep(c(1),6), rep(c(0),22), rep(c(1),51))

    data <- data.frame(present, age)

    > data
        present age
    1         0   0
    2         0   0
    3         0   0
    .         .   .
    .         .   .
    .         .   .
    100       1   1

有没有一种简单的方法可以从表/矩阵中获取这种结构?

【问题讨论】:

  • 手工编码似乎是一次性的最简单方法。如果您认为您会多次执行此操作,您可以构建一个函数,将矩阵/表的元素作为上述rep() 函数的参数。如果矩阵的尺寸保持为 2X2,这将特别容易。如果它们发生变化,那么它将涉及更多的编码。
  • age &lt;- c(rep(c(0),27), rep(c(1),73)) 中定义的年龄变量与表定义不匹配。 55岁以上的27人用1编码,55岁以下的73人用0编码。
  • 谢谢!我没注意到。

标签: r glm contingency


【解决方案1】:
reshape2::melt(example)

这会给你,

     Var1  Var2 value
1 Present >= 55    21
2  Absent >= 55     6
3 Present   <55    22
4  Absent   <55    51

您可以轻松地将其用于glm

【讨论】:

  • 您也可以使用as.data.frame(as.table(example)) 获取此信息(但value 列是Freq
  • 另外,我可能错了,但不确定它是否适合在 glm 中轻松使用。我写了一个答案,显示了stackoverflow.com/a/36508547/4598520
【解决方案2】:

您也许可以将countsToCases 函数用作defined here

countsToCases(as.data.frame(as.table(example))) 
#        Var1  Var2
#1    Present >= 55
#1.1  Present >= 55
#1.2  Present >= 55
#1.3  Present >= 55
#1.4  Present >= 55
#1.5  Present >= 55
# ...

如果您愿意,您可以随时将变量重新编码为数字。

【讨论】:

    【解决方案3】:

    我会去:

    library(data.table)
    tab <- data.table(AGED = c(1, 1, 0, 0),
                      CHD = c(1, 0, 1, 0),
                      Count = c(21, 6, 22, 51))
    
    tabExp <- tab[rep(1:.N, Count), .(AGED, CHD)]
    

    编辑:快速解释,因为我花了一些时间才弄明白:

    data.table 对象中.N 存储组的行数(如果与by 分组)或只是整个data.table 的行数,所以在这个例子中

    tab[rep(1:.N, Count)]
    

    tab[rep(1:4, Count)]
    

    最后

    tab[rep(1:4, c(21, 6, 22, 51)]
    

    是等价的。

    与基础 R 相同:

    tab2 <- data.frame(AGED = c(1, 1, 0, 0),
                       CHD = c(1, 0, 1, 0),
                       Count = c(21, 6, 22, 51))
    
    tabExp2 <- tab2[rep(1:nrow(tab2), tab2$Count), c("AGED", "CHD")]
    

    【讨论】:

      【解决方案4】:

      下面的代码可能看起来很长,但只有 group_by()do() 指令处理扩展数据。剩下的就是以长格式更改数据并将字符变量编码为 0 和 1。我尝试从您在问题中给出的确切矩阵开始。

      加载数据操作包

      library(tidyr)
      library(dplyr)
      

      创建数据框

      按照您的示例创建一个矩阵,但避免在列名中使用“>”符号

      example <- matrix(c(21,22,6,51), nrow = 2, byrow = TRUE)
      rownames(example) <- c("Present","Absent")
      colnames(example) <- c("above55", "below55")
      

      将矩阵转换为数据框

      example <- data.frame(example) %>%
          add_rownames("chd")    
      

      或者干脆直接创建一个数据框

      data.frame(chd = c("Present", "Absent"),
                 above55 = c(21,6),
                 below55 = c(22,51))
      

      重塑数据

      data2 <- example %>% 
          gather(age, nrow, -chd) %>%
          # Encode chd and age as 0 or 1
          mutate(chd = ifelse(chd=="Present",1,0),
                 age = ifelse(age=="above55",1,0)) %>%
          group_by(chd, age) %>%
          # Expand each variable by nrow
          do(data.frame(chd = rep(.$chd,.$nrow),
                        age = rep(.$age,.$nrow)))
      
      head(data2)
      # Source: local data frame [6 x 2]
      # Groups: chd, age [1]
      # 
      #     chd   age
      #   (dbl) (dbl)
      # 1     0     0
      # 2     0     0
      # 3     0     0
      # 4     0     0
      # 5     0     0
      # 6     0     0
      tail(data2)
      # Source: local data frame [6 x 2]
      # Groups: chd, age [1]
      # 
      #     chd   age
      #   (dbl) (dbl)
      # 1     1     1
      # 2     1     1
      # 3     1     1
      # 4     1     1
      # 5     1     1
      # 6     1     1
      
      table(data2)
      #        age
      # chd  0  1
      #   0 51  6
      #   1 22 21
      

      与您的示例相同,但年龄编码除外 我在上面的评论中提到的问题。

      【讨论】:

        【解决方案5】:

        所以,glm 并不是那么死板。部分?glm读取

         For ‘binomial’ and ‘quasibinomial’ families the response can also
         be specified as a ‘factor’ (when the first level denotes failure
         and all others success) or as a two-column matrix with the columns
         giving the numbers of successes and failures.
        

        我假设您想测试年龄对Present/Absent 的影响。 关键是要指定响应,如(在伪代码中)c(success, failure)

        所以你需要像data.frame(Age= ..., Present = ..., Absent) 这样的数据。从example 执行此操作的最简单方法是转置,然后强制转换为data.frame,并添加一列:

        example_t <- as.data.frame(t(example))
        example_df <- data.frame(example_t, Age=factor(row.names(example_t)))
        

        给你

              Present Absent   Age
        >= 55      21      6 >= 55
        <55        22     51   <55
        

        然后,您可以运行 glm:

        glm(cbind(Present, Absent) ~ Age, example_df, family = 'binomial')
        

        得到

        Call:  glm(formula = cbind(Present, Absent) ~ Age, family = "binomial",
            data = example_for_glm)
        
        Coefficients:
        (Intercept)       Age<55
              1.253       -2.094
        
        Degrees of Freedom: 1 Total (i.e. Null);  0 Residual
        Null Deviance:      18.7
        Residual Deviance: -1.332e-15   AIC: 11.99
        

        附录

        您也可以通过@therimalaya 的回答到达这里。但这只是第一步

        as.data.frame(as.table(example))
        

        (只会让你走到那里)

            Var1  Var2 Freq
        1 Present >= 55   21
        2  Absent >= 55    6
        3 Present   <55   22
        4  Absent   <55   51
        

        但要真正拥有成功和失败的专栏,您需要做更多的事情。您可以使用tidyr 到达那里

        as.data.frame(as.table(example)) %>% tidyr::spread(Var1, Freq)
        

        类似于我上面的example_df

          Var2 Present Absent
        1 >= 55      21      6
        2   <55      22     51
        

        【讨论】:

          猜你喜欢
          • 2016-05-15
          • 1970-01-01
          • 2010-11-19
          • 2021-08-03
          • 1970-01-01
          • 1970-01-01
          • 2021-07-16
          • 1970-01-01
          相关资源
          最近更新 更多