我收到了一位审阅者的评论,他希望在人口统计学特征表(表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)
我需要具有Bonferroni校正的变量status
,edema
和stage
的每个级别的p值。我没有成功通过documentation。另外,使用卡方比较行之间的样本大小是否正确?
UPDATE:
我不确定我的方法是否正确,但是我想与您分享。我为变量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
中提取pValue,有一种调用对象属性的方法(首先说明),还有一种快速而肮脏的方法(在底部)。
要提取它们,首先我复制您的代码以创建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()
查看“ tableOne”对象的内容>
attributes(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变量的pValue。如之前设置的那样,您可以同时提取
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
话虽如此,并且提取了您的pValues,但我认为有一种更便捷,更快捷的方法来实现相同目的:如果在tableOne vignette中查看“导出”部分,则可以使用
print(tableOne, quote = TRUE)
,然后只需复制并粘贴到电子表格(如LibreOffice,Excel ...)。然后,我将选择带有pValue的列,对其进行转置,然后将其返回给R,以使用p.adjust()
计算调整后的p值并将其复制回电子表格中以进行日记提交