傾向性の検定

• Rの傾向性の検定関数prop.trend.test()を見てみよう
  • 以下に示すように、xというベクトルとnというベクトルとscoreというカテゴリの「値ベクトル」が入力
  • いわゆる分割表的に言えば、x,n-xがケースとコントロールのベクトルで、nはケース・コントロールの和のベクトル
  • scoreは、カテゴリに与えた値であり、「まっすぐなトレンド」のときは、"1,2,3,4,..."と言う値付けをする
    • SNP 2x3 で優性モデルにするときは、"1,2,2"(とか、もっとわかりやすく"0,1,1"とか)にする
  • 内部ではカテゴリごとのケースの全体に対する比率freqとカテゴリの値scoreとの間でlm()線形回帰している
  • 線形回帰をするときにweights = w = n/(p(1-p))を与えている。ここでpはケース数がケース＋コントロールの総数に占める比率である
  • lm()では、sum(w*e^2)が最小になるように最小二乗法を用いている
  • このlm()の結果をanova()してカイ二乗分布に従う統計量chisqを出している
  • トレンド検定では、「カテゴリの値、傾け方」はベクトルで指定(score)していて、その「傾け方」の形をそのままに急傾斜にしたり緩傾斜にしたりしているので、自由度は１

> prop.trend.test
function (x, n, score = seq_along(x)) 
{
    method <- "Chi-squared Test for Trend in Proportions"
    dname <- paste(deparse(substitute(x)), "out of", deparse(substitute(n)))
    dname <- paste(dname, ",\n using scores:", paste(score, collapse = " "))
    x <- as.vector(x)
    n <- as.vector(n)
    score <- as.vector(score)
    freq <- x/n
    p <- sum(x)/sum(n)
    w <- n/p/(1 - p)
    a <- anova(lm(freq ~ score, weights = w))
    chisq <- a["score", "Sum Sq"]
    names(chisq) <- "X-squared"
    df <- c(df = 1)
    pval <- pchisq(chisq, 1, lower.tail = FALSE)
    rval <- list(statistic = chisq, parameter = df, p.value = as.numeric(pval), 
        method = method, data.name = dname)
    class(rval) <- "htest"
    return(rval)
}
<bytecode: 0000000009DF3480>
<environment: namespace:stats>