Passing & Bablok 法による回帰直線のパラメータ     Last modified: Sep 03, 2009

目的

Passing & Bablok 法による回帰直線のパラメータを求める。

使用法

PassingBablok(x, y, n=1, PB=TRUE)

# print メソッド
print.passing.bablok(obj, digits=5, sig=0.95)

# plot メソッド
plot.passing.bablok(obj, posx="topleft", posy=NULL, xlab=obj$names.xy[1], ylab=obj$names.xy[2],	hist=FALSE, ...)

引数

x            独立変数ベクトル
y            従属変数ベクトル
n            ブートストラップ法で信頼区間を求めるときの回数
             デフォルトは 1。この場合には信頼区間を求めない
PB            Passing & Bablck 法に厳格に従わないときに FALSE と指定する(使用例参照)
obj           "passing.bablok" オブジェクト
digits        表示桁数
sig           ブートストラップ法による信頼区間を求めるときの信頼度(デフォルトで 95% 信頼限界)
posx, posy    legend 関数のための位置引数("topleft", "bottomright" なども可)
xlab, ylab    軸の名前
hist          ブートストラップ法の結果のヒストグラム表示をするかどうか
...           plot 関数に渡されるその他の引数

ソース

インストールは,以下の 1 行をコピーし,R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/PassingBablok.R", encoding="euc-jp")

# Passing & Bablok 法による回帰直線のパラメータ推定
PassingBablok <- function(   x,                                      # 独立変数ベクトル
                                y,                                      # 従属変数ベクトル
                                n=1,                                    # ブートストラップ法で信頼区間を求めるときの回数
                                PB=TRUE)                                # Passing & Bablock 法に厳格に従わないときに FALSE と指定する
{
        PassingBablok0 <- function(x, y)                             # 1 組のデータについて,切片と傾きの推定値を計算する
        {
                suffix <- combn(length(x), 2)
                slope <- diff(matrix(y[suffix], 2))/diff(matrix(x[suffix], 2)) # すべての二点の組合せの傾きを計算
                slope <- slope[!is.nan(slope)]                               # 傾きが垂直(計算上は NaN になる)のものを除く
                if (PB) {                                               # Passing & Bablok 法では傾きが -1 になるものを除き,
                        slope <- slope[slope != -1]
                        k <- sum(slope < -1)                              # 傾きが -1 未満のものの個数が必要
                }
                else {
                        k <- 0                                               # 負の相関の様な場合にも適切?な推定値を得るためには PB=FALSE を指定する
                }
                slope <- sort(slope)                                 # メディアンを推定値にするが,ちょっと特殊なことをする
                n <- length(slope)
                if (n %% 2 == 0) {                                      # 普通のメディアンの計算より,k 番目に大きいものを使う
                        Slope <- (slope[(n+1)/2+k]+slope[(n+1)/2+k+1])/2
                }
                else {
                        Slope <- slope[(n+1)/2+k]
                }
                return(c(Intercept=median(y-Slope*x), Slope=Slope))     # 結果を返す(2 要素を持つベクトル)
        }
        Driver <- function(x, y)                                     # ブートストラップ法のためのドライバー
        {
                n <- length(x)
                suffix <- sample(n, n, replace=TRUE)                 # リサンプリング
                return(PassingBablok0(x[suffix], y[suffix]))            # リサンプリングしたデータについてパラメータを推定
        }
        names.xy <- c(deparse(substitute(x)), deparse(substitute(y)))        # 変数名を控えておく
        OK <- complete.cases(x, y)                                   # 欠損値を持つケースを除く
        x <- x[OK]
        y <- y[OK]
        ans <- list(coefficients=PassingBablok0(x, y),                       # 引数に対してパラメータを推定する
                    names.xy=names.xy, x=x, y=y)
        if (n > 1) {
                ans2 <- replicate(n, Driver(x, y))                   # ブートストラップを n 回実行
                ans <- append(ans, list(intercepts=ans2[1,],  slopes=ans2[2,]))
        }
        class(ans) <- "passing.bablok"                                       # print, plot メソッドのためにクラス名をつけておく
        return(ans)
}
# print メソッド
print.passing.bablok <- function(    obj,                            # "passing.bablok" オブジェクト
                                        digits=5,                       # 表示桁数
                                        sig=0.95)                       # 信頼度
{
        if (length(obj) == 4) {
                cat("Intercept:", round(obj$coefficients[1], digits),
                    "    Slope:", round(obj$coefficients[2], digits), "\n")
        }
        else {
                alpha <- (1-sig)/2
                LCL <- c(quantile(obj$intercepts, alpha), quantile(obj$slopes, alpha))
                UCL <- c(quantile(obj$intercepts, 1-alpha), quantile(obj$slopes, 1-alpha))
                ans <- data.frame(obj$coefficients, LCL, UCL)
                dimnames(ans) <- list(c("Intercept:", "Slope:"),
                                      c("Estimate", paste(c(alpha, 1-alpha), "%", sep="")))
                print(ans, digits=digits)               
        }
}
# plot メソッド
plot.passing.bablok <- function(     obj,                            # "passing.bablok" オブジェクト
                                        posx="topleft", posy=NULL,      # legend 関数のための位置引数
                                        xlab=obj$names.xy[1],           # x 軸の名前
                                        ylab=obj$names.xy[2],           # y 軸の名前
                                        hist=FALSE,                     # ヒストグラムを描くとき TRUE にする
                                        ...)                            # その他の任意の plot 関数の引数
{
        if (hist && length(obj) == 6) {                                   # ブートストラップの結果を,hist=TRUE のときに,ヒストグラムで表示する
                layout(matrix(1:2, 2))
                hist(obj$intercepts, xlab="Intercept", main="", right=FALSE)
                hist(obj$slopes, xlab="Slope", main="", right=FALSE)
                layout(1)
        }
        else {                                                          # 散布図と Passing & Bablok 法の回帰直線と直線回帰式を表示する
                plot(obj$x, obj$y, xlab=xlab, ylab=ylab, ...)
                abline(obj$coefficients)
                abline(lm(obj$y~obj$x), lty=2, col=2)
                legend(posx, posy, legend=c("Passing-Bablok", "linear regression"), lty=1:2, col=1:2)
        }
}


使用例

> set.seed(123)

> d <- gendat2(100, 0.7)

> (a <- PassingBablok(d[,1], d[,2])) # 点推定値のみを求める

Intercept: 0.02861     Slope: 1.05069 

> print(a, digits=7)

Intercept: 0.0286131     Slope: 1.050686 

> plot(a, xlab="x", ylab="y", pch=19) # 散布図と回帰直線

fig

> (b <- PassingBablok(d[,1], d[,2], n=100)) # 点推定値と信頼区間を求める

           Estimate   0.025%  0.975%
Intercept: 0.028613 -0.18854 0.20789
Slope:     1.050686  0.87304 1.27921

> plot(b, hist=TRUE) # ブートストラップパラメータの分布図

fig

本来 Passing & Bablok 法は,同じ対象についての 2 方法による測定値の相関を見るものであり,
下図のような負の相関を持つようなデータには対応しない。

fig

PB 引数を FALSE とすれば,負の相関を持つデータの回帰もできる。
(原著者はこの点については言及がない。正しいものかどうかは保証しない)

fig


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh