R的Cochran-Armitage趋势检验

Question

我有这样的数据：(encoding_0、encoding_1和encoding_2列是对每个snp_id:值的is_severe列进行编码0,1,2的次数的计数

X            snp_id         is_severe encoding_1 encoding_2 encoding_0
1     0  GL000191.1-37698         0          0          1          7
2     1  GL000191.1-37698         1          0          2         11
3     3 GL000192.1-100085         0          5          3          0
4     4 GL000192.1-100085         1          3         10          0

当二项式值为is_severe列，有序数据(已在频率表中(为以下列时，我希望执行Cochran-Armitage趋势测试：encoding_0、encoding_1、encoding _2。我发现了这个执行测试的功能：

catt_2 <-
function(y, x, score = c(0, 1, 2)) {
miss <- unique(c(which(is.na(y)), which(is.na(x))))
n.miss <- length(miss)
if(n.miss > 0) {
y <- y[-miss]
x <- x[-miss]
}
if(!all((y == 0) | (y == 1))) 
stop("y should be only 0 or 1.")
if(!all((x == 0) | (x == 1) |(x == 2))) 
stop("x should be only 0, 1 or 2.")
ca <- x [y == 1]
co <- x [y == 0]
htca <- table(ca)
htco <- table(co)
A <- matrix(0, 2, 3)
colnames(A) <- c(0, 1, 2)
rownames(A) <- c(0, 1)
A[1, names(htca)] <- htca
A[2, names(htco)] <- htco
ptt <- prop.trend.test(A[1, ], colSums(A), score = score)
res <- list(
chisq = as.numeric(ptt$statistic), 

p.value = as.numeric(ptt$p.value)
)
return(res)
}

但是自变量X是一个向量，它不包含频率，只包含编码类型(0,1,2(，我已经有了频率表。正如我所理解的频率计算和测试性能(使用prop.trend.test-https://search.r-project.org/R/refmans/stats/html/prop.trend.test.html(发生在这里：

htca <- table(ca)
htco <- table(co)
A <- matrix( 2, 2)
colnames(A) <- c(1,2)
rownames(A) <- c(0, 1)
A[1, names(htca)] <- htca
A[2, names(htco)] <- htco
ptt <- prop.trend.test(A[1, ], colSums(A), score = score)

我想知道是否有可能改变这个功能，并使其适应我已经有了频率的情况。(我必须为每个snp_id执行该函数，为此我将使用by((函数(例如，我尝试将一些数据转换为矩阵：

ï..is_severe encod_0 encod_1 encod_2
[1,]            0       1       2       5
[2,]            1       3       2       8

然后做了：

prop.trend.test(my_mat,colSums(my_mat))

它返回了这个：

prop.trend.test(my_mat,colSums(my_mat))
Error in model.frame.default(formula = freq ~ score, data = list(freq = x/n,  : 
variable lengths differ (found for '(weights)')

谢谢：(

Answer 1

catt <-
function(y, x0,x1,x2, score = c(0, 1, 2)) {
h=data.frame(is_severe=c(0,1),
encod_0=c(0,0),
encod_1=c(0,0),
encod_2=c(0,0)
)
h$is_severe=y
h$encod_0=x0
h$encod_1=x1
h$encod_2=x2
my_mat <- as.matrix(h)
A <- matrix(0, 2, 3)
colnames(A) <- c(0, 1, 2)
rownames(A) <- c(0, 1)
A[1,] =my_mat[1,2:4]
A[2, ]=my_mat[2,2:4]
ptt <- prop.trend.test(A[1, ],colSums(A), score = score) 
p.value = as.numeric(ptt$p.value) 
res <- list(
chisq = as.numeric(ptt$statistic), 

p.value = as.numeric(ptt$p.value)
)
return(res)

}
catt_2 <-
function(y, score = c(0, 1, 2)) {

A <-y
colnames(A) <- c(0, 1, 2)
rownames(A) <- c(0, 1)
ptt <- prop.trend.test(A[1, ], colSums(A), score = score)
res <- list("2x3-table"=A,
chisq = as.numeric(ptt$statistic), 

p.value = as.numeric(ptt$p.value)
)
return(res)
}