我收到一位审稿人的评论,他希望获得人口统计特征表(表 1)中每一行特定变量水平的所有 p 值。尽管这个要求对我来说显得很奇怪(而且不准确),但我愿意同意他的建议。
library(tableone)
## Load data
library(survival); data(pbc)
# drop ID from variable list
vars <- names(pbc)[-1]
## Create Table 1 stratified by trt (can add more stratifying variables)
tableOne <- CreateTableOne(vars = vars, strata = c("trt"), data = pbc, factorVars = c("status","edema","stage"))
print(tableOne, nonnormal = c("bili","chol","copper","alk.phos","trig"), exact = c("status","stage"), smd = TRUE)
我需要为变量 status
、edema
和 stage
的每个级别提供 p 值,并进行 Bonferroni 校正。我浏览了 documentation没有成功。
此外,使用卡方来比较跨行的样本大小是否正确?
更新:
我不确定我的方法是否正确,但我想与您分享。我为变量 status
生成了每个层的虚拟变量,然后计算了 chisq
。
library(tableone)
## Load data
library(survival); data(pbc)
d <- pbc[,c("status", "trt")]
# Convert dummy variables
d$status.0 <- ifelse(d$status==0, 1,0)
d$status.1 <- ifelse(d$status==1, 1,0)
d$status.2 <- ifelse(d$status==2, 1,0)
t <- rbind(
chisq.test(d$status.0, d$trt),
# p-value = 0.7202
chisq.test(d$status.1, d$trt),
# p-value = 1
chisq.test(d$status.2, d$trt)
#p-value = 0.7818
)
t
用于多重比较的 BONFERRONI ADJ:
p <- t[,"p.value"]
p.adjust(p, method = "bonferroni")
最佳答案
这个问题是前一段时间发布的,所以我想你已经回答了审稿人。
我真的不明白为什么只为三个变量计算调整后的 p 值。事实上,调整 p 值取决于进行比较的次数。如果您将 p.adjust()
与包含 3 个 p 值的向量一起使用,则结果实际上不会根据所进行的比较量进行“调整”(您确实进行了十几个半!)
我展示了如何提取所有 p 值,以便您计算调整后的值。
要从包 tableOne
中提取 pValues 有一种调用对象属性的方法(首先解释),和两种快速但肮脏的方法(在底部).
要提取它们,首先我复制您的代码以创建您的 tableOne:
library(tableone)
## Load data
library(survival); data(pbc)
# drop ID from variable list
vars <- names(pbc)[-1]
## Create Table 1 stratified by trt (can add more stratifying variables)
tableOne <- CreateTableOne(vars = vars, strata = c("trt"), data = pbc, factorVars = c("status","edema","stage"))
您可以通过 attributes()
attributes(tableOne)
你可以看到一个tableOne通常有一个表格用于连续变量和分类变量。您也可以在其中使用 attributes()
attributes(tableOne$CatTable)
# you can notice $pValues
现在您知道 pValue 的“位置”,您可以使用 attr()
attr(tableOne$CatTable, "pValues")
与数字变量类似的东西:
attributes(tableOne$ContTable)
# $pValues are there
attr(tableOne$ContTable, "pValues")
您有 Normal 和 NonNormal 变量的 pValues。 正如您之前设置的那样,您可以提取两者
mypCont <- attr(tableOne$ContTable, "pValues") # put them in an object
nonnormal = c("bili","chol","copper","alk.phos","trig") # copied from your code
mypCont[rownames(mypCont) %in% c(nonnormal), "pNonNormal"] # extract NonNormal
"%!in%" <- Negate("%in%")
mypCont[rownames(mypCont) %!in% c(nonnormal), "pNonNormal"] # extract Normal
综上所述,并且提取了您的 pValue,我认为有两种更方便快捷但又肮脏的方法可以完成同样的任务:
快速而肮脏的方法 A:使用 dput()
打印 tableOne。然后在控制台中搜索 pValues 所在的位置并将它们复制粘贴到脚本中,以将它们存储在一个对象中
快速而肮脏的方式 B:如果您查看 tableOne vignette有一个“导出”部分,您可以使用 print(tableOne, quote = TRUE)
然后只需复制并粘贴到电子表格(如 LibreOffice、Excel...)。
然后我会选择带有 pValue 的列,转置它,然后将其返回到 R,以使用 p.adjust()
计算调整后的 p 值,并将它们复制回电子表格以提交期刊
关于r - 库 "TableOne"多重比较。逐行计算 p 值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/52024849/