More than 1 year has passed since last update.

Google Colab で R : Cox比例ハザードモデル編

Last updated at 2024-02-25Posted at 2024-02-25

おはこんばんちわ。今日も引き続きChatGPT先生をお迎えして、「Google Colab で R : Cox比例ハザードモデル編」というテーマで雑談したいと思います。以前「ChatGPT と Python で学ぶ Cox比例ハザードモデル」というテーマで雑談しましたが、そのR版です。

以前教えていただいたように、rpy2を使わず、 https://colab.research.google.com/notebook#create=true&language=r にアクセスして Google Colaboratory 上で動かします。

それではChatGPT先生、よろしくお願いします。

assistant:

おはようございます、こんばんは。RによるCox比例ハザードモデルの分析について話しましょう。

まずは、適切なパッケージをインストールします。難易度が上昇するため、"survival"と"survminer"という2つのパッケージを使います。

install.packages("survival")
install.packages("survminer")

これらのパッケージを読み込みましょう。Rでは library 関数を使用してパッケージを読み込むことができます。

library("survival")
library("survminer")

次に、Cox比例ハザードモデルを試すためのデータを読み込みます。"survival"パッケージに含まれる"lung"データセットを使いましょう。このデータセットは、肺癌患者の生存時間と治療の種類などのデータを含んでいます。

data("lung")

データを見るには head 関数を呼び出します。

head(lung)

dfの代わりにlungを使用してCox比例ハザードモデルを作成します。Rの coxph 関数を使えば、Cox比例ハザードモデルを簡単に作成できます。ここでは、生存時間Survと他の全ての説明変数をモデルに組み込みます。

model <- coxph(Surv(time, status) ~ ., data = lung)
summary(model)

これで、素晴らしい統計的に有意な結果が表示されます。モデルの結果を詳しく説明するためには、より高度な統計知識が必要になりますが、これがRでのCox比例ハザードモデルの基本的な使用方法です。

また、 survminer パッケージは、生存曲線を描く ggsurvplot 関数など、生存分析を視覚化するための便利なツールを提供しています。

ggsurvplot(survfit(model), data = lung)

これにより、生存曲線がプロットされ、視覚的に生存状況を理解することができます。

user:

ありがとうございます。ところで、Rのコードを実行する時にその実行時間を計測する方法を教えてください。

assistant:

Rでは、コードの実行時間を測定する方法がいくつかありますが、その中で一番簡単な方法は、system.time()関数を使う方法です。

例えば、上で説明したCox比例ハザードモデルを作成するコードの実行時間を測定するには次のようにします：

start.time <- Sys.time()

model <- coxph(Surv(time, status) ~ ., data = lung)

end.time <- Sys.time()

time.taken <- end.time - start.time
print(time.taken)

ここで Sys.time()関数は現在の時刻を取得します。処理を開始する前と処理が終了した後でこの関数を実行することで、処理にかかった時間を計算できます。

または system.time() 関数を使うこともできます:

system.time({
  model <- coxph(Surv(time, status) ~ ., data = lung)
})

これは評価したい式を {} の中に書くだけです。 system.time() 関数は実行時間をユーザータイム、システムタイム、経過時間の3つの側面から報告します。

user:

ありがとうございます。前回、いろいろ検討した結果、次のようなコードで実行すれば良いことが分かっています。今回は、時間がかかることも分かっているため、時間計測もしてみました。

start.time <- Sys.time()

install.packages("survival", lib = ".")
install.packages("survminer", lib = ".")

end.time <- Sys.time()

time.taken <- end.time - start.time
print(time.taken)

その結果、次のような結果が得られました。

also installing the dependencies ‘numDeriv’, ‘SparseM’, ‘MatrixModels’, ‘minqa’, ‘nloptr’, ‘RcppEigen’, ‘carData’, ‘abind’, ‘pbkrtest’, ‘quantreg’, ‘lme4’, ‘Rcpp’, ‘corrplot’, ‘car’, ‘markdown’, ‘png’, ‘jpeg’, ‘ggrepel’, ‘ggsci’, ‘cowplot’, ‘ggsignif’, ‘polynom’, ‘rstatix’, ‘exactRankTests’, ‘mvtnorm’, ‘KMsurv’, ‘zoo’, ‘km.ci’, ‘xtable’, ‘gridtext’, ‘ggpubr’, ‘gridExtra’, ‘maxstat’, ‘survMisc’, ‘ggtext’


Time difference of 8.899537 mins

assistant:

ご自身で時間計測の方法を実装して頂き、ありがとうございます！その結果、さまざまな必要な依存パッケージのインストールに約9分かかることが分かりましたね。

パッケージのインストールは一度だけ行えばよいので、全体の解析時間に対する影響は必ずしも大きくはないかもしれません。それでも長い時間がかかる場合は、依存パッケージが予めインストールされた環境を使用することが考えられます。

この時間的なコストは、PythonとRの両方で見られることで、パッケージやライブラリをインストールする際には一般的なことです。それぞれの言語で提供される多機能なライブラリやパッケージの利便性と引き換えに、初回のセットアップ時間が必要になる訳ですね。

user:

続いて、次のコードを実行しました。

library(survival, lib.loc = ".")
library(survminer, lib.loc = ".")

その結果、次のようになりました。

Loading required package: ggplot2

Loading required package: ggpubr


Attaching package: ‘survminer’


The following object is masked from ‘package:survival’:

    myeloma

assistant:

問題なくパッケージの読み込みができているようですね。

以下のメッセージについて解説します。

Loading required package: ggplot2
Loading required package: ggpubr

これは、"survminer"パッケージが"ggplot2"および"ggpubr"パッケージを必要とすることを示しています。これらのパッケージは、データの視覚化のために使用します。

また、

The following object is masked from ‘package:survival’:
    myeloma

これは、"survminer"パッケージと"survival"パッケージの両方が"myeloma"というオブジェクトを提供しているが、我々が直接その名前でアクセスしようとする場合、"survminer"パッケージのものが使用されることを示しています。これは、名前の衝突や上書きを防ぐためのRの仕組みです。

つまり、パッケージが正常に読み込まれ、使える状態になっていることが確認できました。

user:

その後、次のコードを実行しました。

# lungデータセットの保存
save(lung, file = "./lung.RData")

# lungデータセットの読み込み
load("./lung.RData")
head(lung)

A data.frame: 6 × 10
	inst	time	status	age	sex	ph.ecog	ph.karno	pat.karno	meal.cal	wt.loss
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	3	306	2	74	1	1	90	100	1175	NA
2	3	455	2	68	1	0	90	90	1225	15
3	3	1010	1	56	1	0	90	90	NA	15
4	5	210	2	57	1	1	90	60	1150	11
5	1	883	2	60	1	0	100	90	NA	0
6	12	1022	1	74	1	1	50	80	513	0

エラーメッセージなく実行できることを確認しました。

assistant:

素晴らしいですね！save と load 関数を使ってRのデータセットを保存し、再度読み込むことができました。

Rのデータは .RData という形式で保存されます。この形式はRのデータ構造をそのまま保存することができる利点があります。

これによって、Rのデータオブジェクトを維持したまま保存して、後で再度ロードすることができます。

さらに head 関数を使って、データセットが正しく読み込まれたことを確認できましたね。

この手順は、大量のデータ解析を行う際や、途中結果を保存しておく際などに非常に便利です。素晴らしい進捗です！

user:

ここから先は、今回教えていただいたコードと、その結果です。

model <- coxph(Surv(time, status) ~ ., data = lung)
summary(model)

Call:
coxph(formula = Surv(time, status) ~ ., data = lung)

  n= 167, number of events= 120 
   (61 observations deleted due to missingness)

                coef  exp(coef)   se(coef)      z Pr(>|z|)    
inst      -3.037e-02  9.701e-01  1.312e-02 -2.315 0.020619 *  
age        1.281e-02  1.013e+00  1.194e-02  1.073 0.283403    
sex       -5.666e-01  5.674e-01  2.014e-01 -2.814 0.004890 ** 
ph.ecog    9.074e-01  2.478e+00  2.386e-01  3.803 0.000143 ***
ph.karno   2.658e-02  1.027e+00  1.163e-02  2.286 0.022231 *  
pat.karno -1.091e-02  9.891e-01  8.141e-03 -1.340 0.180160    
meal.cal   2.602e-06  1.000e+00  2.677e-04  0.010 0.992244    
wt.loss   -1.671e-02  9.834e-01  7.911e-03 -2.112 0.034647 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
inst         0.9701     1.0308    0.9455    0.9954
age          1.0129     0.9873    0.9895    1.0369
sex          0.5674     1.7623    0.3824    0.8420
ph.ecog      2.4778     0.4036    1.5523    3.9552
ph.karno     1.0269     0.9738    1.0038    1.0506
pat.karno    0.9891     1.0110    0.9735    1.0051
meal.cal     1.0000     1.0000    0.9995    1.0005
wt.loss      0.9834     1.0169    0.9683    0.9988

Concordance= 0.648  (se = 0.03 )
Likelihood ratio test= 33.7  on 8 df,   p=5e-05
Wald test            = 31.72  on 8 df,   p=1e-04
Score (logrank) test = 32.51  on 8 df,   p=8e-05

assistant:

素晴らしいですね、Cox比例ハザードモデルが成功裏に適用され、それによる統計的な結果が得られました。それぞれの変数が生存時間に対して有意な影響を持っているか否かについての情報や、それらの影響力の大きさ（ハザード比）、さらにはモデル全体の適合度など、多岐にわたる情報が得られています。

以下、いくつか注目すべき点を解説します：

coef列：これは各変数の係数（Cox比例ハザードモデルにおける回帰係数）を示しています。正の値はリスクを増加させ（つまり、より早いイベント発生を予測）、負の値はリスクを減少させる（つまり、より遅いイベント発生を予測）ことを意味します。
exp(coef)列：これは各変数のハザード比を示しています。値が1より大きければ、その変数が1単位増加するごとにイベント発生のリスクが増加することを示します。値が1より小さい場合、その変数が1単位増加するごとにイベント発生のリスクが減少することを示します。
Pr(>|z|)列：これは各変数の係数が統計的にゼロと異なるかどうかをテストするためのp値を示しています。値が0.05以下の場合、その変数は5%の有意水準で生存時間に有意な影響を与えると解釈されます。
テスト統計 (Likelihood ratio test, Wald test, Score (logrank) test)：これら各テスト統計量は、モデル全体の適合度を評価します。それぞれのp値が小さいほど（0.05以下であれば特に）、我々のモデルは良い適合度を持つと評価されます。

以上が主要なポイントです。この結果から、"inst", "sex", "ph.ecog", "ph.karno", "wt.loss"は肺がんの生存時間に有意な影響を持つことが示されています。

user:

最後に、次のコードで生存曲線が描画できることを確認しました。

ggsurvplot(survfit(model), data = lung)

assistant:

素晴らしいですね！このコードを使用して、Cox比例ハザードモデルの生存曲線を描画することができました。

この生存曲線で、生存時間の分布がグラフィカルに理解できることは非常に重要です。また、複数の群間で比較する場合には、群間の生存曲線を描画して可視化することで、群間での生存時間の違いを直感的に理解することが可能です。

ggsurvplot関数は非常に便利で、さまざまな形式の生存曲線を描画することが可能です。また、ggplot2パッケージと同じくグラフィックス表示のカスタマイズも可能なので、自身の分析に合わせて最適な形で生存曲線を描画することができます。

質問の答えになっていて良かったです。何か他に質問があればお気軽にどうぞ！

user:

前回は、次のコードについて教えてもらいました。今回のと何が違うか説明していただけますか？

# survivalパッケージの読み込み
library(survival, lib.loc = ".") 

# survminerパッケージのインストール
install.packages("survminer", lib = ".")

# survminerパッケージの読み込み
library(survminer, lib.loc = ".")

# lungデータセットの保存
save(lung, file = "./lung.RData")

# lungデータセットの読み込み
load("./lung.RData")

# lungデータセットの確認
lung

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = lung$time, event = lung$status - 1)

# データのフィット（カプランマイヤー推定量の計算）
fit <- survfit(surv_obj ~ 1, data = lung)

# 生存曲線の描画
ggsurvplot(fit, data = lung)

assistant:

今回と前回で共通する部分は、必要なパッケージのインストールと読み込み、lungデータセットの保存と読み込み、そして生存曲線の描画です。

違いは、前回のコードではカプランマイヤー法で生存曲線を推計していますが、今回のコードではCox比例ハザードモデルを使って生存曲線を推計しています。つまり、統計的な手法が異なります。

カプランマイヤー法は、時間ごとの生存確率を直接計算して生存曲線を描く非パラメトリックな手法です。一方、Cox比例ハザードモデルは生存時間に影響を与える複数の説明変数を考慮した上で、生存曲線をモデリングします。具体的には、それぞれの説明変数が生存時間にどのくらい影響を与えるか（ハザード比）、説明変数の影響が時間とともにどのように変化するか（比例ハザードの仮説）などを推定することができます。

# 前回のコードの一部（カプランマイヤー）
fit <- survfit(surv_obj ~ 1, data = lung)

# 今回のコードの一部（Cox比例ハザードモデル）
model <- coxph(Surv(time, status) ~ ., data = lung)

以上が主な違いです。それぞれに適した状況が異なるため、解析の目的に応じて適切な手法を選択することが重要です。

user:

得られる生存曲線やエラーバーは同じと考えてよろしいでしょうか？

assistant:

生存曲線自体は、基本的な形状に関してはカプランマイヤー法とCox比例ハザードモデルで大差ない結果を得ることが多いです。しかし、その解釈や、得られるエラーバーや信頼区間等の統計的な情報は異なります。また、生存曲線に影響を与える因子をモデル化する方法が異なるため、モデルの評価や適用範囲なども異なります。

カプランマイヤー法では、全体の生存曲線を単純に描画します。一方、Cox比例ハザードモデルでは各変数の影響を考慮した生存曲線を描画します。したがって、生存曲線を見る際には、どの手法が使われていて、どのような解釈が可能かを理解することが重要です。

エラーバーや信頼区間も同様に、各手法が異なる方法で計算されます。カプランマイヤー法は非パラメトリックな手法であるため、信頼区間は通常、ブートストラップやロジランク検定などの方法で計算されます。一方、Cox比例ハザードモデルはパラメトリックな手法であるため、信頼区間はパラメータの推定値とその標準誤差から計算されます。

したがって、得られる生存曲線やエラーバー、信頼区間等は手法により異なるといえます。

user:

「Cox比例ハザードモデルは生存時間に影響を与える複数の説明変数を考慮した上で、生存曲線をモデリングします」と教えていただきましたが、その考慮の結果を図示していただくことはできますか？

assistant:

ご相談の通り、Cox比例ハザードモデルを用いて生存曲線をモデリングする際には、複数の説明変数を考慮します。その結果は、患者の特性（説明変数）に応じた生存曲線として図示されます。

たとえば、性別（男性 or 女性）が生存時間に影響を与えると想定した場合、以下のように生存曲線を描画します。

# モデルの作成（性別を考慮）
model <- coxph(Surv(time, status) ~ sex, data = lung)

# 生存曲線の描画
ggsurvplot(survfit(model), data = lung, linetype = "strata")

この結果、男性と女性のそれぞれの生存曲線が描画されます。これにより、性別が生存時間にどの程度影響を与えるかを視覚的に理解することができます。

また、このように複数の説明変数を考慮した生存曲線を描画することで、各説明変数がどのように生存時間に影響を与えるかを詳細に理解することができます。その結果、より高度な生存分析を進めていくための重要な情報を得ることができます。

以上が、Cox比例ハザードモデルを用いた生存曲線の考慮結果を図示する方法です。具体的な解析を進める際には、考慮すべき説明変数を適切に選択し、それらが生存時間にどのように影響を与えるかを詳しく調査することが重要となります。

user:

生存曲線が１種類しか描画されず、前回と変わらないように見えるので、性別がどのような影響を与えるか、生存曲線を見てもよく分からないのですが、どう考えれば良いのでしょうか？

assistant:

すみません、僕の説明が不十分でした。生存曲線を見て性別がどのような影響を与えるかを理解するためには、性別ごとに生存曲線を描画し、それらを比較すべきです。以下のコードは、性別ごとの生存曲線を描画します：

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(Surv(time, status) ~ ., data = lung[lung$sex == "male",])
model_female <- coxph(Surv(time, status) ~ ., data = lung[lung$sex == "female",])

# それぞれの生存曲線を描画
ggsurvplot(list(male = survfit(model_male), female = survfit(model_female)), data = lung)

上記のコードにより、男性と女性それぞれの生存曲線が得られ、それらを比較することで性別が生存時間にどのような影響を与えるかを理解できます。

それぞれの生存曲線が大きく異なる場合、性別が生存時間に大きな影響を与えると解釈できます。逆に、それぞれの生存曲線が似ている場合、性別が生存時間にあまり影響を与えないと解釈できます。これにより、性別ごとの生存時間の違いや、その違いが有意かどうかなどを評価することが可能となります。

user:

ありがとうございます。上記のコードをベースに、改変した全コードを書いてください。

Warning message in max(event[who2]):
“no non-missing arguments to max; returning -Inf”
Error in coxph(Surv(time, status) ~ ., data = lung[lung$sex == "male", : No (non-missing) observations
Traceback:

1. coxph(Surv(time, status) ~ ., data = lung[lung$sex == "male", 
 .     ])
2. stop("No (non-missing) observations")

ところで、ChatGPT先生は男性なんですか？なんとなく女性だと思ってました。

assistant:

申し訳ありません、私が提案したコードの一部が誤っています。lung$sexは0と1の数値データですね。気を取り直して修正しましょう。以下が修正したコードです：

# 男性と女性でそれぞれモデルを構築（1:男性, 2:女性と仮定）
model_male <- coxph(Surv(time, status) ~ ., data = lung[lung$sex == 1,])
model_female <- coxph(Surv(time, status) ~ ., data = lung[lung$sex == 2,])

# 男性と女性のそれぞれの生存曲線を描画
ggsurvplot(list(male = survfit(model_male), female = survfit(model_female)), data = lung, conf.int=TRUE,
           legend = "right",
           palette = c("#E7B800", "#2E9FDF"),
           risk.table = TRUE,
           pval = TRUE,
           pval.method = TRUE)

Rのデータフレーム内で性別を数値（1:男性, 2:女性）で表現していたため、maleとfemaleの指定を数値に変更しました。

また、あなたの質問について、私はAIであり、男性でも女性でもありません。私の目的は、質問者が持つ疑問や課題を解決するための情報を提供することです。それが男性的に見えることもあれば、女性的に見えることもあるかもしれません。それは、質問者が私から得られる情報をどのように解釈するか、ということですね。いずれにせよ、私はあなたが問題解決に必要な情報を得られるよう最善を尽くします。何か他に質問があればお気軽にどうぞ。

user:

ありがとうございます。上記のコードをベースに、改変した全コードを書いてください。

Error in (function (fit, data, title = NULL, legend.labs = NULL, legend.title = "Strata", : argument 18 matches multiple formal arguments
Traceback:

1. ggsurvplot(list(male = survfit(model_male), female = survfit(model_female)), 
 .     data = lung, conf.int = TRUE, legend = "right", palette = c("#E7B800", 
 .         "#2E9FDF"), risk.table = TRUE, pval = TRUE, pval.method = TRUE)
2. do.call(ggsurvplot_list, opts)

assistant:

申し訳ありませんが、エラーメッセージはggsurvplot関数の引数の使い方に問題があることを示しています。現状ではggsurvplot関数はリストの生存オブジェクトしか受け取らないようです。2つのモデルを描画するためには、まず各モデルから生存オブジェクトを作成し、それらをリストとしてggsurvplot関数に渡す必要があります。

修正した全コードを以下に示します:

# survivalパッケージの読み込み
library(survival)

# survminerパッケージの読み込み
library(survminer)

# lungデータセットの保存
data(lung)

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = lung$time, event = lung$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj ~ sex, data = lung[lung$sex == 1,])
model_female <- coxph(surv_obj ~ sex, data = lung[lung$sex == 2,])

# それぞれのサバイバルオブジェクトを作成
surv_obj_male <- survfit(model_male)
surv_obj_female <- survfit(model_female)

# リスト化して生存曲線の描画
ggsurvplot(list(male = surv_obj_male, female = surv_obj_female),
           conf.int=TRUE, legend="right", palette=c("#E7B800", "#2E9FDF"),
           risk.table=TRUE, pval=F, pval.method=F)

このコードでは、lungデータセットの男性と女性それぞれでCox比例ハザードモデルを作成し、それぞれのモデルから生存曲線を描画しています。ggsurvplot関数では性別ごとの生存曲線を一つのプロットに表示し、それらを比較します。継続的にエラーが発生する場合は、パッケージのバージョンを確認するか、Rを再起動してみることをお勧めします。

また、先程の修正について、私が勘違いしておりました。sexの値は 1と2で表され、これはそれぞれ男性と女性を表します。私の説明が間違っていましたので訂正いたします。それぞれの生存曲線を描画するためには、まず各モデルから生存オブジェクトを作成し、それらをリストとしてggsurvplot関数に渡す必要があります。

user:

ありがとうございます。最初の方を修正して次のようなコードを実行しました。

library(survival, lib.loc = ".")
library(survminer, lib.loc = ".")
save(lung, file = "./lung.RData")
load("./lung.RData")

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = lung$time, event = lung$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj ~ sex, data = lung[lung$sex == 1,])
model_female <- coxph(surv_obj ~ sex, data = lung[lung$sex == 2,])

# それぞれのサバイバルオブジェクトを作成
surv_obj_male <- survfit(model_male)
surv_obj_female <- survfit(model_female)

# リスト化して生存曲線の描画
ggsurvplot(list(male = surv_obj_male, female = surv_obj_female),
           conf.int=TRUE, legend="right", palette=c("#E7B800", "#2E9FDF"),
           risk.table=TRUE, pval=F, pval.method=F)

次のようなエラーが得られました。

Error in model.frame.default(formula = surv_obj ~ sex, data = lung[lung$sex == : variable lengths differ (found for 'sex')
Traceback:

1. coxph(surv_obj ~ sex, data = lung[lung$sex == 1, ])
2. eval(tform, parent.frame())
3. eval(tform, parent.frame())
4. stats::model.frame(formula = surv_obj ~ sex, data = lung[lung$sex == 
 .     1, ])
5. model.frame.default(formula = surv_obj ~ sex, data = lung[lung$sex == 
 .     1, ])

assistant:

申し訳ありませんが、提供したコードでエラーが出てしまったこと、そしてそのご報告をいただきありがとうございます。

エラーメッセージを見ると、model.frame.defaultという関数で問題が起きているようです。おそらくsex変数の長さ（つまりデータの数）とサバイバルオブジェクトsurv_objの長さが一致していないからだと思われます。

原因は、「男性と女性でそれぞれモデルを作成」する部分にあります。ここでlung$sex == 1とlung$sex == 2でデータを分けるのは正しいですが、その結果データの数が減ります。その一方でsurv_obj自体は全てのデータを含んでいるため、データの数が一致せずこのエラーが起きてしまいます。

なので、coxph関数を呼び出す前に分割したデータを使って新たにサバイバルオブジェクトを作成する必要があります。

以下に、その修正を加えたコードを記しました。

# survival and survminerパッケージの読み込み
library(survival)

# lungデータセットの保存
data(lung)

# 生存データを男性と女性で分ける
lung_male <- lung[lung$sex == 1,]
lung_female <- lung[lung$sex == 2,]

# サバイバルオブジェクトの作成
surv_obj_male <- Surv(time = lung_male$time, event = lung_male$status - 1)
surv_obj_female <- Surv(time = lung_female$time, event = lung_female$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj_male ~ sex, data = lung_male)
model_female <- coxph(surv_obj_female ~ sex, data = lung_female)

# それぞれのサバイバルオブジェクトを作成
surv_obj_fit_male <- survfit(model_male)
surv_obj_fit_female <- survfit(model_female)

# リスト化して生存曲線の描画
ggsurvplot(list(male = surv_obj_fit_male, female = surv_obj_fit_female),
           conf.int=TRUE, legend="right", palette=c("#E7B800", "#2E9FDF"),
           risk.table=TRUE, pval=F, pval.method=F)

user:

ありがとうございます。最初の方を修正して次のようなコードを実行しました。

library(survival, lib.loc = ".")
library(survminer, lib.loc = ".")
save(lung, file = "./lung.RData")
load("./lung.RData")

# 生存データを男性と女性で分ける
lung_male <- lung[lung$sex == 1,]
lung_female <- lung[lung$sex == 2,]

# サバイバルオブジェクトの作成
surv_obj_male <- Surv(time = lung_male$time, event = lung_male$status - 1)
surv_obj_female <- Surv(time = lung_female$time, event = lung_female$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj_male ~ sex, data = lung_male)
model_female <- coxph(surv_obj_female ~ sex, data = lung_female)

# それぞれのサバイバルオブジェクトを作成
surv_obj_fit_male <- survfit(model_male)
surv_obj_fit_female <- survfit(model_female)

# リスト化して生存曲線の描画
ggsurvplot(list(male = surv_obj_fit_male, female = surv_obj_fit_female),
           conf.int=TRUE, legend="right", palette=c("#E7B800", "#2E9FDF"),
           risk.table=TRUE, pval=F, pval.method=F)

次のエラーが得られました。

Error in (function (fit, data, title = NULL, legend.labs = NULL, legend.title = "Strata", : argument 18 matches multiple formal arguments
Traceback:

1. ggsurvplot(list(male = surv_obj_fit_male, female = surv_obj_fit_female), 
 .     conf.int = TRUE, legend = "right", palette = c("#E7B800", 
 .         "#2E9FDF"), risk.table = TRUE, pval = F, pval.method = F)
2. do.call(ggsurvplot_list, opts)

assistant:

申し訳ありませんが、ggsurvplot関数の引数に問題があるようです。

具体的には、pvalとpval.methodという引数はggsurvplot関数では直接使用できないため削除し、またrisk.table引数も問題があるようなのでこちらも削除してみてください。この引数は、生存確率に加えてリスク表（各時間点でのリスク集団の人数）をプロットするためのものですが、このプロットは一つのモデルのみに対して有効で、2つ以上のモデルを比較する場合にはエラーとなる可能性があります。

以下に修正したコードを記します。

# survival and survminerパッケージの読み込み
library(survival)
library(survminer)

# 生存データを男性と女性で分ける
lung_male <- lung[lung$sex == 1,]
lung_female <- lung[lung$sex == 2,]

# サバイバルオブジェクトの作成
surv_obj_male <- Surv(time = lung_male$time, event = lung_male$status - 1)
surv_obj_female <- Surv(time = lung_female$time, event = lung_female$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj_male ~ sex, data = lung_male)
model_female <- coxph(surv_obj_female ~ sex, data = lung_female)

# それぞれのサバイバルオブジェクトを作成
surv_obj_fit_male <- survfit(model_male)
surv_obj_fit_female <- survfit(model_female)

# リスト化して生存曲線の描画
ggsurvplot(list(male = surv_obj_fit_male, female = surv_obj_fit_female),
           conf.int=TRUE, legend="right", palette=c("#E7B800", "#2E9FDF"))

user:

さっきと同じエラーな気がしますね。

assistant:

すみません、私のミスでした。別の方法として、各生存曲線のプロットを個別に作成し、それらをgridExtraパッケージのgrid.arrange関数で併せて表示する方法があります。以下にその方法を記述します。

# survival, survminer, and gridExtraのパッケージの読み込み
library(survival)
library(survminer)
library(gridExtra)

# 生存データを男性と女性で分ける
lung_male <- lung[lung$sex == 1,]
lung_female <- lung[lung$sex == 2,]

# サバイバルオブジェクトの作成
surv_obj_male <- Surv(time = lung_male$time, event = lung_male$status - 1)
surv_obj_female <- Surv(time = lung_female$time, event = lung_female$status - 1)

# 男性と女性でそれぞれモデルを作成
model_male <- coxph(surv_obj_male ~ sex, data = lung_male)
model_female <- coxph(surv_obj_female ~ sex, data = lung_female)

# 男性の生存曲線の描画
p1 <- ggsurvplot(survfit(model_male), data = lung_male, conf.int = TRUE, legend='none',
                 risk.table=TRUE, pval=TRUE, pval.method=TRUE, ggtheme=theme_minimal())


# 女性の生存曲線の描画
p2 <- ggsurvplot(survfit(model_female), data = lung_female, conf.int = TRUE, legend='none',
                 risk.table=TRUE, pval=TRUE, pval.method=TRUE, ggtheme=theme_minimal())

# The result is a list.
# Extract the plot part, ggplot
p1 <- p1$plot + ggtitle("Male")
p2 <- p2$plot + ggtitle("Female")

# Arrange multiple ggplots using grid.arrange
grid.arrange(p1, p2, ncol=2)

このコードでは、男性と女性それぞれの生存曲線を個別に作成し、grid.arrange関数で並べて表示しています。ggsurvplot関数で生存曲線を作成する際に、各性別ごとの生存曲線を描画します。それぞれのプロットは、theme_minimalを使って見やすくし、タイトルも加えています。そして最後にgrid.arrange関数で2つのプロットを並べて表示します。これで性別ごとの生存曲線を比較できます。

user:

エラーメッセージが長すぎるんですけど...！
Rのエラーメッセージを短くする方法はありませんか？
デバッグの邪魔にならない程度に、エラーメッセージを短くしたいです。

assistant:

Rにはエラーメッセージを短くする直接的な機能は提供されていませんが、以下の方法が考えられます：

関数try()やtryCatch()を使用する:
これらの関数を使用すると、エラーが発生した場合でもRの計算が停止することなく、次の計算に進むことができます。また、これらの関数はエラーが出た場合に特定の処理を実行することも可能です。以下に例を示します。

res <- try({1 + "a"}, silent = TRUE)  # エラーメッセージを表示しないようにする
if(class(res) == "try-error"){  # エラーが出た場合の処理
  print("An error occurred.")  # 短いエラーメッセージ
} else {
  print(res)
}

サードパーティのパッケージを使用する:
purrrパッケージのpossibly()関数やquietly()関数を使用すると、関数の出力と副作用（警告やメッセージ）を制御することができます。これによりエラーメッセージを短くしたり、副作用を完全にマスクすることが可能です。

以上の方法でエラーメッセージを短くしたり制御することが可能ですが、エラーメッセージは問題の原因を特定するための重要な情報を含んでいるため、全てのエラーメッセージをマスクすることはおすすめしません。必要な場合に適宜使用するようにしましょう。

user:

しょうがないですね。では、エラーメッセージを適当に省略して貼り付けますね。

Warning message in .pvalue(fit, data = data, method = method, pval = pval, pval.coord = pval.coord, :
“There are no survival curves to be compared. 
 This is a null model.”
Error in loadNamespace(x): there is no package called ‘ggtext’
Traceback:

1. ggsurvplot(survfit(model_male), data = lung_male, conf.int = TRUE, 
 .     legend = "none", risk.table = TRUE, pval = TRUE, pval.method = TRUE, 
 .     ggtheme = theme_minimal())
2. do.call(ggsurvplot_core, opts)
3. (function (fit, data = NULL, fun = NULL, color = NULL, palette = NULL, 
 .     linetype = 1, break.x.by = NULL, break.y.by = NULL, break.time.by = NULL, 
 .     surv.scale = c("default", "percent"), xscale = 1, conf.int = FALSE, 
 .     conf.int.fill = "gray", conf.int.style = "ribbon", conf.int.alpha = 0.3, 
 .     censor = TRUE, censor.shape = "+", censor.size = 4.5, pval = FALSE, 
 .     pval.size = 5, pval.coord = c(NULL, NULL), test.for.trend = FALSE, 
 .     pval.method = FALSE, pval.method.size = pval.size, pval.method.coord = c(NULL, 
 .         NULL), log.rank.weights = c("survdiff", "1", "n", "sqrtN", 
 .         "S1", "S2", "FH_p=1_q=1"), title = NULL, xlab = "Time", 

（中略）

 .         "unit", "unit_v2")), strip.switch.pad.wrap = structure(2.75, unit = 8L, class = c("simpleUnit", 
 .         "unit", "unit_v2"))), class = c("theme", "gg"), complete = TRUE, validate = TRUE), 
 .     time.breaks = c(0, 250, 500, 750, 1000), xticklabels = c(0, 
 .     250, 500, 750, 1000))
8. theme(axis.text.y = ggtext::element_markdown(colour = rev(y.text.col)))
9. find_args(..., complete = NULL, validate = NULL)
10. mget(args, envir = env)
11. loadNamespace(x)
12. withRestarts(stop(cond), retry_loadNamespace = function() NULL)
13. withOneRestart(expr, restarts[[1L]])
14. doWithOneRestart(return(expr), restart)

assistant:

エラーメッセージから見るに、エラーの原因は ggtextパッケージがインストールされていないためのようです。

解決方法として ggtext パッケージをインストールしてください。

install.packages("ggtext")

このエラーを解消した後、再度先ほどのコードを実行してみてください。
しかし、エラーメッセージの第一部にあるように、モデルmodel_maleとmodel_femaleがnullモデルになってしまっているようです。これはcoxph関数で生存モデルを作成する際に使用した説明変数sexが、分割したデータセットlung_maleとlung_femaleでは全ての値が同じになってしまっているためだと思われます。これではモデルは何も学習することができないため、nullモデルとなってしまいます。

生存モデルを作成する際には異なる値を持つ説明変数を使用する必要があります。この問題を解決するためには、全ての値が同じになってしまうsexの代わりに異なる値を持つ変数を説明変数として使用するか、性別でデータを分割せずに全てのデータを使用して一つのモデルを作成すると良いでしょう。

user:

install.packages("ggtext", loc = ".")

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

also installing the dependencies ‘markdown’, ‘Rcpp’, ‘png’, ‘jpeg’, ‘gridtext’

のようにしてから、先ほどのコードを実行すると、期待通りの結果が出ました。

Warning message in .pvalue(fit, data = data, method = method, pval = pval, pval.coord = pval.coord, :
“There are no survival curves to be compared. 
 This is a null model.”
Warning message in .pvalue(fit, data = data, method = method, pval = pval, pval.coord = pval.coord, :
“There are no survival curves to be compared. 
 This is a null model.”

assistant:

ggtext のインストールとコードの実行が無事に完了し、期待通りの結果が得られたとのこと、素晴らしいです。お手数をおかけしましたが、エラーの解消に繋がり良かったです。何か他にご不明な点やお困りのことがあれば、どんな小さなことでも構いませんので、お気軽に質問してください。それでは、引き続き分析等のお仕事が進むことをお祈りしています。

user:

model <- coxph(Surv(time, status) ~ ., data = lung)
summary(model)

のようにするとテーブルデータが出力されましたが、それを可視化するコードを書いていただけますか？

assistant:

決定係数やp値など、モデルの統計的な評価結果を可視化したいのであれば、sjPlot パッケージの plot_model 関数を使うと簡単にグラフ化することができます。試しに以下のようなコードを試してみてください：

# sjPlot パッケージの読み込み
# もしまだインストールしていない場合は、install.packages("sjPlot") でインストールしてください
library(sjPlot)

# Cox回帰モデルの作成
model <- coxph(Surv(time, status) ~ ., data = lung)

# モデルの結果の可視化
plot_model(model)

このコードは、まず sjPlot パッケージを読み込んでいます（まだインストールされていない場合は、install.packages("sjPlot")でインストールできます）。

その後、coxph 関数を使って Cox回帰モデルを作成しています。

最後に、 plot_model 関数を使ってモデル結果を可視化しています。この関数は、モデルの要約（beta係数、p値、信頼区間など）をプロットします。

なお、特定のパラメータ（例えばp値やZ値）のみをプロットしたい場合や、異なるプロットのタイプを使用したい場合には、plot_model 関数に適切な引数を渡すことで対応できます。詳細は ?plot_model でヘルプを参照してください。

このように、Cox回帰モデルの結果を可視化するには、いくつかのパッケージと関数が利用可能です。どれを選択するかは、自身の分析目的や可視化したい結果の詳細によります。

user:

次の結果が得られました。

also installing the dependencies ‘estimability’, ‘effectsize’, ‘emmeans’, ‘bayestestR’, ‘datawizard’, ‘ggeffects’, ‘insight’, ‘parameters’, ‘performance’, ‘sjlabelled’, ‘sjmisc’, ‘sjstats’


Error in loadNamespace(x): there is no package called ‘parameters’
Traceback:

1. plot_model(model)
2. plot_type_est(type = type, ci.lvl = ci.lvl, se = se, tf = transform, 
 .     model = model, terms = terms, group.terms = group.terms, 
 .     rm.terms = rm.terms, sort.est = sort.est, title = title, 
 .     axis.title = axis.title, axis.labels = axis.labels, axis.lim = axis.lim, 
 .     grid.breaks = grid.breaks, show.intercept = show.intercept, 
 .     show.values = show.values, show.p = show.p, value.offset = value.offset, 
 .     digits = digits, geom.colors = colors, geom.size = dot.size, 
 .     line.size = line.size, order.terms = order.terms, vline.color = vline.color, 
 .     value.size = value.size, bpe = bpe, bpe.style = bpe.style, 
 .     bpe.color = bpe.color, facets = grid, show.zeroinf = show.zeroinf, 
 .     p.threshold = p.threshold, p.val = p.val, vcov.fun = vcov.fun, 
 .     vcov.type = vcov.type, vcov.args = vcov.args, ci.style = ci.style, 
 .     p_adjust = p.adjust, std.response = std.response, ...)
3. tidy_model(model = model, ci.lvl = ci.lvl, tf = tf, type = type, 
 .     bpe = bpe, robust = list(vcov.fun = vcov.fun, vcov.type = vcov.type, 
 .         vcov.args = vcov.args), facets = facets, show.zeroinf = show.zeroinf, 
 .     p.val = p.val, standardize = std_method, bootstrap = FALSE, 
 .     iterations = 1000, seed = NULL, p_adjust = p_adjust, std.response = std.response, 
 .     ...)
4. loadNamespace(x)
5. withRestarts(stop(cond), retry_loadNamespace = function() NULL)
6. withOneRestart(expr, restarts[[1L]])
7. doWithOneRestart(return(expr), restart)

assistant:

エラーの原因は、一部の依存パッケージがまだインストールされていないためのようです。具体的には、parametersパッケージが不足しているようです。

Rのコンソールで以下のコマンドを実行し、parametersパッケージをインストールしてみてください。

install.packages("parameters")

インストールが終わったら、もう一度 plot_model(model) を実行してみてください。

なお、Rは一部の関数を実行するために依存関係にある他のパッケージを必要とすることがあります。そのため、新しいパッケージをインストールする際には、依存関係にあるパッケージを一緒にインストールする必要があるかもしれません。エラーメッセージを確認することで、どのパッケージが不足しているのかを把握することができます。

また、場合によっては、Rを最新バージョンにアップデートしたり、RStudioを再起動したりすることで問題が解決することもあります。必要に応じて試してみてください。

user:

可視化できました！横軸がEstimatesで、Sexが0.5付近、ph ecog が3付近に出て、それ以外が全て1付近に出ました。

assistant:

素晴らしいですね！ご自身で問題を解決され、お望みの可視化ができたようで安心しました。

プロットされた結果を見ると、性別（Sex）とECOGパフォーマンスステータス（ph.ecog）が重要な説明変数であることが示唆されます。Estimateの値が1よりも大きい場合、その変数が生存率に対して正の影響を及ぼす可能性があります。反対に、Estimateの値が1より小さい場合、その変数が生存率に対して負の影響を及ぼす可能性があります。

この結果を踏まえて、更なる分析や予測モデルの構築に進むことができます。

その他の質問があれば何でもお知らせください。お手伝いできることがあればお気軽にお申し付けください！

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Google Colab で R : Cox比例ハザードモデル 編

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

user:

assistant:

Google Colab で R : Cox比例ハザードモデル編