Rで解析：信頼性・生存時間分析のためのL-moment統計「lmomco」パッケージ

保健医療や生物統計学分野では、生存時間データや信頼性分析に取り組む機会が多くあります。このような領域では、分布の特性やパラメータを推定することが重要な課題です。例えば、製品の故障率や患者の生存率を分析する際には、正確な分布モデリングとパラメータ推定が必要となります。

このような時におすすめなのが「lmomco」パッケージです。このパッケージでは、L-Momentsや分布のパラメータ推定、LM比率図作成、多変量Lコモメント解析などが可能です。特に、右側・左側打ち切り（censoring）に対応したL-momentの計算や、13種類の分位関数演算子（quantile function operators）により、信頼性・生存時間分析に対応しています。

他のL-moment関連パッケージ（lmomやLmomentsなど）と比較した強みは、対応する分布の種類の豊富さと、打ち切りやトリムL-momentといった拡張的な解析手法までを一貫したインターフェースでカバーする総合性です。目的に応じた分布フィッティングや信頼性分析を効率よく実行できるため、研究や業務の効率化に貢献します。

パッケージバージョンは2.5.7。Windows 11 x64 (build 26200)のR version 4.6.0で確認しています。

＜おすすめのRに関する書籍です＞

Rではじめるデータサイエンス 第2版 | Hadley Wickham, Mine Çetinkaya-Rundel, Garrett Grolemund, 大橋 真也

Rによる機械学習[第3版] (Programmer's SELECTION)

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("lmomco")
# パッケージの読み込み
library("lmomco")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

信頼性係数（L-alpha）の計算：headrick.sheng.lalphaコマンド

Headrick and Sheng (2013) の論文中のTable 1データを用いて、観測値どうしの相関構造からL-alpha（信頼性係数）を算出する例します。引数の指定方法（Lコモメント行列を直接与える、共分散パスウェイを使う、非超過確率を利用しないなど）によって、得られるalpha値がどう変わるかを確認できます。

### データ作成#####
# Headrick and Sheng (2013) の Table 1
# 1〜3列目が観測値、4〜6列目が推定非超過確率
TV1 <-
  c(2, 4, 3, 0.15, 0.45, 0.15,       5, 7, 7, 0.75, 0.95, 1.00,
    3, 5, 5, 0.35, 0.65, 0.40,       6, 6, 6, 0.90, 0.80, 0.75,
    7, 7, 6, 1.00, 0.95, 0.75,       5, 2, 6, 0.75, 0.10, 0.75,
    2, 3, 3, 0.15, 0.25, 0.15,       4, 3, 6, 0.55, 0.25, 0.75,
    3, 5, 5, 0.35, 0.65, 0.40,       4, 4, 5, 0.55, 0.45, 0.40)
T1 <- matrix(ncol=6, nrow=10)
for(r in seq(1,length(TV1), by=6)) T1[(r/6)+1, ] <- TV1[r:(r+5)]
colnames(T1) <- c("X1", "X2", "X3", "FX1", "FX2", "FX3"); T1 <- as.data.frame(T1)

lco2 <- matrix(nrow=3, ncol=3)
lco2[1,1] <- lmoms(T1$X1)$lambdas[2]
lco2[2,2] <- lmoms(T1$X2)$lambdas[2]
lco2[3,3] <- lmoms(T1$X3)$lambdas[2]
lco2[1,2] <- 2*cov(T1$X1, T1$FX2); lco2[1,3] <- 2*cov(T1$X1, T1$FX3)
lco2[2,1] <- 2*cov(T1$X2, T1$FX1); lco2[2,3] <- 2*cov(T1$X2, T1$FX3)
lco2[3,1] <- 2*cov(T1$X3, T1$FX1); lco2[3,2] <- 2*cov(T1$X3, T1$FX2)
########

# Headrick and Sheng (2013) の結果：alpha = 0.807
headrick.sheng.lalpha(lco2)$alpha
L-alpha 
0.8074766 

headrick.sheng.lalpha(Lcomoment.matrix(T1[,1:3], k=2)$matrix)$alpha
L-alpha 
0.7805825 

# FXsを使用しない場合：alpha = 0.781
headrick.sheng.lalpha(T1[,1:3])$alpha
L-alpha 
0.7805825

# a=0.5（デフォルトはHazen法）
headrick.sheng.lalpha(T1[,1:3], bycovFF=TRUE)$alpha
L-alpha 
0.7805825

headrick.sheng.lalpha(T1[,1:3], bycovFF=TRUE, a=0.5)$alpha
L-alpha 
0.7805825

正規分布に基づく正常性検定：tau34sq.normtestコマンド

# 正規分布の母平均と分散の検定を行う
HarriCoble <- tau34sq.normtest(rnorm(20), verbose=TRUE)
#        *** Harri-Coble Tau34-squared Test for Normality ***
#  --A normality test using sample L-skew (T3) and L-kurtosis (T4)--
#  Z-score(T3) = -0.4828 and Z-score(T4) = -0.114
#    (Mapping of T3 and T4 to standard normal variates)
#  Tau34-squared test statistic, T34sq = Z(T3)^2 + Z(T4)^2
#    (A squared Euclidean distance that is Chi-Squared(df=2) distributed.)
#  T34sq = 0.2461 |--> ChiSquared(T34sq, df=2) = pchisq(T34sq, df=2)
#  ChiSquared(0.2461, 2 degrees freedom) = 11.58 percentile
#    p-value = 0.8842 [complement of the pchisq(T34sq, df=2, lower.tail=TRUE)]
#    The p-value is > alpha: accept Ho, conclude 'normal data'.
#  Reference:
#   Harri, A. and Coble, K.H., 2011, Normality testing---
#     Two new tests using L-moments: J. Appl. Stat. 38(7), 1369-1379.

＜おすすめのRに関する書籍です＞

医学研究のためのRによる統計解析入門 | 加葉田 大志朗, 新谷 歩

Amazonで加葉田 大志朗, 新谷 歩の医学研究のためのRによる統計解析入門。アマゾンならポイント還元本が多数。

www.amazon.co.jp

累積分布関数を計算し、プロット：z.par2cdfコマンド

# 乱数を固定
set.seed(21)
the.gpa   <- vec2par(c(100,1000,0.1),type='gpa')

# データをシミュレーション
fake.data <- rlmomco(30,the.gpa)

# ゼロの観測値を追加
fake.data <- sort(c(fake.data,rep(0,10)))

# 目盛りを内側に向け、右軸にもタイトルを付けるためプロットパラメータを変更
par(mgp=c(3,0.5,0), mar=c(5,4,4,3))

# 正の値のデータについてパラメータを計算
gpa.all <- pargpa(lmoms(fake.data))
gpa.nzo <- pargpa(lmoms(fake.data[fake.data > 0]))

# サンプルサイズ
n <- length(fake.data)

# ゼロ値の推定確率
p <- length(fake.data[fake.data == 0])/n
F <- nonexceeds(sig6=TRUE); F <- sort(c(F,p)); qF <- qnorm(F)

# 以下のxベクトルには元々ゼロが含まれるため、挿入は不要
# x<0では不合理だが、実装のテストのため使用
x <- seq(-100, max(fake.data))

# シミュレーションデータのプロット位置を計算
PP <- pp(fake.data)

# サンプルをプロット
plot(fake.data, qnorm(PP), xlim=c(0,4000), yaxt="n", ylab="")
add.lmomco.axis(las=2, tcl=0.5, side=2, twoside=FALSE,
                                        side.type="NPP", otherside.type="SNV")

# 母集団分布（ゼロを除く）
lines(quagpa(F,gpa.all), qF)
cdf <- qnorm(z.par2cdf(x,p,gpa.nzo))

# 上記のヘルプを参照
cdf[! is.finite(cdf)] <- min(fake.data,qnorm(PP))

# ゼロを条件としたフィッティングモデル
lines(x, cdf,lwd=3)

# 上記のコードを繰り返し実行し、結果の変化を確認
# デフォルト設定に戻す
par(mgp=c(3,1,0), mar=c(5,4,4,2)+0.1)

四分位函数を計算しプロット：z.par2quaコマンド

# 母集団パラメータ（真の値に近いもの）を定義
the.gpa   <- vec2par(c(100,1000,0.1),type='gpa')

# データをシミュレーション
fake.data <- rlmomco(30,the.gpa)

# ゼロの観測値を追加
fake.data <- sort(c(fake.data, rep(0,10)))

# 目盛りを内側に向けるためパラメータを変更
par(mgp=c(3,0.5,0))

# 正の値のデータについてパラメータを計算
gpa.all <- pargpa(lmoms(fake.data))
gpa.nzo <- pargpa(lmoms(fake.data[fake.data > 0]))

# サンプルサイズ
n <- length(fake.data)

# ゼロ値の推定確率
p <- length(fake.data[fake.data == 0])/n
F <- nonexceeds(sig6=TRUE); F <- sort(c(F,p)); qF <- qnorm(F)

# シミュレーションデータのプロット位置を計算
PP <- pp(fake.data)

# サンプルをプロット
plot(qnorm(PP), fake.data, ylim=c(0,4000), xaxt="n", xlab="")
add.lmomco.axis(las=2, tcl=0.5, twoside=TRUE, side.type="SNV", otherside.type="NA")

# 母集団分布（ゼロを除く）
lines(qF,quagpa(F,gpa.all))

# ゼロを条件としたフィッティングモデル
lines(qF,z.par2qua(F,p,gpa.nzo),lwd=3)

# デフォルト設定に戻す
par(mgp=c(3,1,0))

lmomco軸を追加してプロット：add.lmomco.axisコマンド

# 目盛りを内側に向けるためパラメータを変更
par(mgp=c(3,0.5,0))

# 合成データを生成
X <- sort(rnorm(65)); pp <- pp(X)

# lmomco軸を追加してプロット
plot(qnorm(pp), X, xaxt="n", xlab="", ylab="QUANTILE", xlim=c(-2,3))
add.lmomco.axis(las=2, tcl=0.5, side.type="RI", otherside.type="NPP")

# デフォルト設定に戻す
par(mgp=c(3,1,0))

＜おすすめのRに関する書籍です＞

初心者でもすぐにできるフリー統計ソフトEZR(Easy R)で誰でも簡単統計解析(改訂第2版) | 神田善伸

Amazonで神田善伸の初心者でもすぐにできるフリー統計ソフトEZR(Easy R)で誰でも簡単統計解析(改訂第2版)。アマゾンならポイント還元本が多数。

www.amazon.co.jp

plotlmrdiaコマンドの使用例

# lmrdia()のプロットを行う
plotlmrdia(lmrdia())

plotlmrdia46コマンドの使用例

# プロット
plotlmrdia46(lmrdia46(), 
              nogld_byt5opt=FALSE, 
              nosymgdd=FALSE, 
              autolegend=TRUE, 
              xleg="topleft")

実行例

データの生成とL-momentの計算

# 乱数シード設定
set.seed(123)

# 正規分布からデータを生成
data <- rnorm(100, mean = 50, sd = 10)

# L-momentを計算
lmoments_data <- lmoms(data)

# L-momentの結果表示
lmoments_data

実行結果：

$lambdas
[1] 50.9040591  5.1781669  0.1381172  0.6560920 -0.0870208

$ratios
[1]          NA  0.10172405  0.02667299  0.12670353 -0.01680533

$trim
[1] 0

$leftrim
NULL

$rightrim
NULL

$source
[1] "lmoms"

＜おすすめのRに関する書籍です＞

RStudioではじめるRプログラミング入門 | ギャレット・グロ-ルマンド, 大橋真也, 長尾高弘

Amazonでギャレット・グロ-ルマンド, 大橋真也, 長尾高弘のRStudioではじめるRプログラミング入門。アマゾンならポイント還元本が多数。

www.amazon.co.jp

この記事が誰かの役に立ちますように。