4

我有一些数据包含大约 2000 个人的组变量 (0/1) 和个人分数。数据集如下所示：

ID group score  
A1 1 3.5  
A2 1 3.2  
A3 0 2.8  
A4 0 2.5

我想测试是否可以通过分数预测组变量，并在 R 中使用了以下内容： 

glm(group~score,family=binomial)

现在我想通过改组组变量来测试我的 p 值，然后再次执行 glm。我想这样做至少 10,000 次，甚至可能更多，每次打印文件中分数的 p 值，以便每个排列有一行。我已经查看了 sample()，但我很难将它与 glm() 以及如何仅输出 p 值结合起来。在脚本/公式中，我想轻松更改排列的数量，如果我选择添加协变量，还可以更改 glm 公式。

感谢您的任何帮助！

4

1 回答 1

3

你在正确的轨道上。

示例（我添加了一个值来抑制有关“数字拟合概率为 0 或 1”的警告）

ex <- read.table(textConnection(
"ID group score  
A1 1 3.5  
A2 1 3.2  
A3 0 2.8  
A4 0 2.5
A5 1 2.4"),header=TRUE)

g0 <- glm(group~score,data=ex,family=binomial)

现在您需要一个函数来计算汇总 p 值（您可以在 中即时执行此操作replicate，但这种方式更简洁）。

pvalfun <- function() {
   g <- update(g0,data=transform(ex,group=sample(group)))
   coef(summary(g))["score","Pr(>|z|)"]
}
res <- replicate(1000,pvalfun())

或者

library(plyr)
res <- raply(1000,pvalfun(),.progress="text")

或者

library(glmperm)
ptest2 <- prr.test(group~score,"score",data=ex,family=binomial)
summary(ptest2)
于 2013-10-07T14:33:07.907 回答