目的

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

使用法

PassingBablok(x, y, n=1)

# 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。この場合には信頼区間を求めない
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) # ブートストラップ法で信頼区間を求めるときの回数
{
  PassingBablok0 <- function(x, y) # 1 組のデータについて，切片と傾きの推定値を計算する
  {
    cor.sign2 = sign(cor(x, y, method="kendall"))
    suffix <- combn(length(x), 2)
    slope <- diff(matrix(y[suffix], 2))/diff(matrix(x[suffix], 2)) # すべての二点の組合せの傾きを計算
    slope <- slope[!is.infinite(slope)] # 傾きが垂直（計算上は Inf になる）のものを除く
    if (cor.sign2 == cor.sign) { # Passing & Bablok 法では傾きが -1 になるものを除き，
      slope <- slope[slope != -1]
      k <- sum(slope < -1) # 傾きが -1 未満のものの個数が必要
      slope <- sort(slope) # メディアンを推定値にするが，ちょっと特殊なことをする
    }
    else {
      slope <- slope[slope != 1]
      k <- sum(slope > 1)
      slope <- sort(slope, decreasing=TRUE) # メディアンを推定値にするが，ちょっと特殊なことをする
    }
    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]
  cor.sign = sign(cor(x, y, method="kendall"))
  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, na.rm=TRUE), quantile(obj$slopes, alpha, na.rm=TRUE))
    UCL <- c(quantile(obj$intercepts, 1-alpha, na.rm=TRUE), quantile(obj$slopes, 1-alpha, na.rm=TRUE))
    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) # 散布図と回帰直線



> (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) # ブートストラップパラメータの分布図