目的

ED50  や LD50 を求める

使用法

ed50(x, n, r, transform=c("none", "ln", "log"))

引数

x	各群の用量（ベクトル）
n	各群の例数（ベクトル）
r	各群の反応数（ベクトル）
transform
	用量をそのまま使うとき（デフォールト）は "none"
	自然対数変換するときは "ln"
	常用対数変換するときは "log"

ソース

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

# ED50 や LD50 を求める
ed50 <- function(    x,                                      # 各群の用量ベクトル
                        n,                                      # 各群の例数ベクトル
                        r,                                      # 各群の反応数ベクトル
                        transform=c("none", "ln", "log"))       # 対数変換をするかどうか
{
        # 第一近似
        estimation1 <- function(x, n, r)
        {
                select <- r > 0 & r < n                                # 条件に合うものだけを取り出す
                x <- x[select]
                n <- n[select]
                r <- r[select]
                p9 <- r/n                                    # 標本比率
                y <- qnorm(p9)+5                             # プロビット変換
                w <- 1
                wx2 <- sum(w*x^2)
                wy <- sum(w*y)
                wx <- sum(w*x)
                wyx <- sum(w*y*x)
                ww <- length(x)
                temp <- ww*wx2-wx^2
                alpha <- (wx2*wy-wx*wyx)/temp                        # 切片
                beta <- (ww*wyx-wx*wy)/temp                  # 傾き
                ollf <- sum(r*log(p9)+(n-r)*log(1-p9))               # 観察値の対数尤度
                list(ollf=ollf, alpha=alpha, beta=beta)
        }

        # 第二近似（第一近似から出発して収束計算）
        estimation2 <- function(x, n, r, alpha, beta, epsilon=1e-5)
        {
                sqrt.pi2 <- 1/sqrt(2*pi)
                for (ii in 1:100) {
                        y <- alpha+beta*x
                        x9 <- y-5
                        p9 <- pnorm(x9)
                        z <- exp(-x9^2/2)*sqrt.pi2
                        y <- y+ (r/n-p9)/z
                        w <- n*z^2/((1-p9)*p9)
                        wx2 <- sum(w*x^2)
                        wy <- sum(w*y)
                        wx <- sum(w*x)
                        wyx <- sum(w*y*x)
                        ww <- sum(w)
                        ssnwx <- wx2-wx^2/ww
                        beta1 <- (wyx-wx*wy/ww)/ssnwx
                        xbar <- wx/ww
                        alpha1 <- (wy/ww)-beta1*xbar
                        if (abs(alpha-alpha1)/alpha < epsilon && abs(beta-beta1)/beta < epsilon) {
                                ed50 <- (5-alpha1)/beta1     # ED50
                                se <- sqrt(1/ww+(ed50-xbar)^2/ssnwx)/abs(beta1)
                                g <- (qnorm(0.975)/beta1)^2/ssnwx
                                cl <- ed50+g*(ed50-xbar)/(1-g)
                                pm = qnorm(0.975)/beta1/(1-g)*sqrt((1-g)/ww+(ed50-xbar)^2/ssnwx)
                                return(list(alpha=alpha1, beta=beta1, ED50=ed50, SE=se, cll =cl-pm, clu=cl+pm))
                        }
                        alpha <- alpha1
                        beta <- beta1
                }
                stop("収束しませんでした")
        }

        # 適合度の検定
        GoodnessOfFitness <- function(x, n, r, alpha, beta)
        {
                p9 <- pnorm(alpha+beta*x-5)
                p0 <- r/n
                cllf <- sum(r*log(p9)+(n-r)*log(1-p9))               # あてはめ後の対数尤度
                q <- 1-p9
                chi <- sum(n*(p0-p9)^2/(q-q^2))                      # カイ二乗値
                df <- length(x)-2                            # 自由度
                p <- pchisq(chi, df, lower.tail=FALSE)               # P 値
                list(tbl=cbind(x, n, r, n*p9, p0, p9), cllf=cllf, chi=chi, df=df, p=p)
        }

        # 書式付きプリント関数
        printf <- function(fmt, ...) cat(sprintf(fmt, ...))

# ed50 関数本体開始

        transform <- match.arg(transform)                    # 引数の補完
        dose <- "x"                                          # 用量のラベル
        inv <- function(x) x                                 # 無変換のときの逆関数（無駄だけど）
        if (transform != "none") {
                stopifnot(x > 0)
                if (transform == "ln") {
                        x <- log(x)                          # 自然対数変換
                        dose <- "ln(x)"                              # 用量のラベル
                        inv <- function(x) exp(x)            # 自然対数の逆関数
                }
                else if (transform == "log") {
                        x <- log10(x)                                # 常用対数変換
                        dose <- "log10(x)"                   # 用量のラベル
                        inv <- function(x) 10^x                      # 常用対数の逆関数
                }
        }

        result <- estimation1(x, n, r)                               # 第一近似解
        ollf <- result$ollf                                  # 観察値の対数尤度
        alpha <- result$alpha                                        # 切片
        beta <- result$beta                                  # 傾き

        result <- estimation2(x, n, r, alpha, beta,          # 第二近似解
                                epsilon=1e-13)
        alpha <- result$alpha                                        # 切片
        beta <- result$beta                                  # 傾き
        printf("P_hat = %g + %g * %s\n", alpha, beta, dose)
        printf("ED50 = %g\n", inv(result$ED50))
        printf("95%% 信頼区間 = [ %g, %g ]\n\n", inv(result$cll), inv(result$clu))

        result <- GoodnessOfFitness(x, n, r, alpha, beta)
        printf("%7s %7s %7s %10s %10s %10s\n", "用量", "n", "r", "e", "r/n", "e/n")
        for (i in 1:nrow(result$tbl)) {
                temp <- result$tbl[i,]
                printf("%7.5g %7.0f %7.0f %10.5g %10.5f %10.5f\n", temp[1], temp[2], temp[3], temp[4], temp[5], temp[6])
        }
        printf("\n")
        printf("カイ二乗値 = %g, 自由度 = %d, P 値 = %.3f\n\n", result$chi, as.integer(result$df), result$p)
        printf("　　観察値の対数尤度 = %g\n", ollf)
        printf("あてはめ後の対数尤度 = %g\n", result$cllf)
}


使用例

x <- c(10, 20, 30)
n <- rep(10, 3)
r <- c(1, 5, 8)
ed50(x, n, r)
ed50(x, n, r, tr="ln")

出力結果例

> ed50(x, n, r)	# 用量をそのまま使う場合
P_hat = 2.78185 + 0.104719 * x	# α，βの推定値
ED50 = 21.1818	# ED50 の推定値
95% 信頼区間 = [ 15.2076, 28.0387 ]	# ED50 の 95% 信頼区間

   用量       n       r          e        r/n        e/n	# サマリー
     10      10       1     1.2081    0.10000    0.12081
     20      10       5     4.5075    0.50000    0.45075
     30      10       8     8.2211    0.80000    0.82211

カイ二乗値 = 0.172154, 自由度 = 1, P 値 = 0.678	# 適合度の検定の結果

　　観察値の対数尤度 = -15.1863
あてはめ後の対数尤度 = -15.2728

> ed50(x, n, r, tr="ln")	# 用量を自然対数変換する場合
P_hat = -0.745754 + 1.92892 * ln(x)
ED50 = 19.663
95% 信頼区間 = [ 13.8884, 27.9769 ]

   用量       n       r          e        r/n        e/n
 2.3026      10       1    0.96075    0.10000    0.09608
 2.9957      10       5     5.1308    0.50000    0.51308
 3.4012      10       8     7.9243    0.80000    0.79243

カイ二乗値 = 0.0120987, 自由度 = 1, P 値 = 0.912

　　観察値の対数尤度 = -15.1863
あてはめ後の対数尤度 = -15.1924