使用 glm() 代替简单的卡方检验

机器算法验证 r 假设检验 广义线性模型 卡方检验 抵消
2022-01-27 01:22:21

我有兴趣在glm()R中更改零假设。

例如:

x = rbinom(100, 1, .7)  
summary(glm(x ~ 1, family = "binomial"))

的假设如果我想将 null 更改为 = 某个任意值怎么办? p=0.5pglm()

我知道这也可以用prop.test()and来完成chisq.test(),但我想探索使用glm()它来测试与分类数据相关的所有假设的想法。

3个回答

您可以使用偏移量:在对数赔率或 logit 标度上glm使用family="binomial"估计参数,因此对应于 0 的对数赔率或 0.5 的概率。如果要与的概率进行比较,则希望基线值为现在的统计模型β0=0pq=logit(p)=log(p/(1p))

YBinom(μ)μ=1/(1+exp(η))η=β0+q

只有最后一行与标准设置有所不同。在 R 代码中:

  • offset(q)在公式中使用
  • logit/log-odds 函数是qlogis(p)
  • 有点烦人的是,您必须为响应变量中的每个元素提供一个偏移值 - R 不会自动为您复制一个常量值。这是通过设置数据框在下面完成的,但您可以只使用rep(q,100).
x = rbinom(100, 1, .7)
dd <- data.frame(x, q = qlogis(0.7)) 
summary(glm(x ~ 1 + offset(q), data=dd, family = "binomial"))

查看 GLM 参数的置信区间:

> set.seed(1)
> x = rbinom(100, 1, .7)
> model<-glm(x ~ 1, family = "binomial")
> confint(model)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.3426412 1.1862042 

这是对数赔率的置信区间。

对于,我们有的假设相当于检查置信区间是否包含 0。这个不包含,所以假设被拒绝。p=0.5log(odds)=logp1p=log1=0p=0.5

现在,对于任意,您可以计算对数赔率并检查它是否在置信区间内。p

使用基于 glm.summary 函数中 z/t 值的 p 值作为假设检验是不(完全)正确/准确的。

  1. 这是令人困惑的语言。报告的值称为 z 值。但在这种情况下,他们使用估计的标准误差代替真实偏差。因此,实际上它们更接近 t 值比较以下三个输出:
    1)summary.glm
    2)t-test
    3)z-test

    > set.seed(1)
    > x = rbinom(100, 1, .7)
    
    > coef1 <- summary(glm(x ~ 1, offset=rep(qlogis(0.7),length(x)), family = "binomial"))$coefficients
    > coef2 <- summary(glm(x ~ 1, family = "binomial"))$coefficients
    
    > coef1[4]  # output from summary.glm
    [1] 0.6626359
    > 2*pt(-abs((qlogis(0.7)-coef2[1])/coef2[2]),99,ncp=0) # manual t-test
    [1] 0.6635858
    > 2*pnorm(-abs((qlogis(0.7)-coef2[1])/coef2[2]),0,1) # manual z-test
    [1] 0.6626359
    
  2. 它们不是精确的 p 值。使用二项分布精确计算 p 值会更好(以现在的计算能力,这不是问题)。假设误差为高斯分布,t 分布并不精确(它高估了 p,超过 alpha 水平在“现实”中发生的频率较低)。请参阅以下比较:

    # trying all 100 possible outcomes if the true value is p=0.7
    px <- dbinom(0:100,100,0.7)
    p_model = rep(0,101)
    for (i in 0:100) {
      xi = c(rep(1,i),rep(0,100-i))
      model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
      p_model[i+1] = 1-summary(model)$coefficients[4]
    }
    
    
    # plotting cumulative distribution of outcomes
    outcomes <- p_model[order(p_model)]
    cdf <- cumsum(px[order(p_model)])
    plot(1-outcomes,1-cdf, 
         ylab="cumulative probability", 
         xlab= "calculated glm p-value",
         xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
    lines(c(0.00001,1),c(0.00001,1))
    for (i in 1:100) {
      lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
    #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
    }
    
    title("probability for rejection as function of set alpha level")
    

    alpha 拒绝的 CDF

    黑色曲线代表平等。红色曲线在其下方。这意味着对于通过 glm 汇总函数计算的给定 p 值,我们在现实中发现这种情况(或更大的差异)的频率低于 p 值所指示的频率。