Triad sou.

MacとWindowsを両方使うようになって作った環境構築メモ

PC R TeX

近年MacとWindowsを両方使うようになり、今後も同じ設定を使いそうなのでメモにしておきました。

R環境

R StudioはMac版とWindows版で使用感に差がないですね。

WebサイトとWebアプリ作成環境

未だにHTMLとJavaScriptを手書きしてWebサイトとWebアプリを作っており、この作業にはVisual Studio Codeを使っています。
Visual Studio Codeも使用感に差がないです（ログインしておくと設定を同期できるので便利）。
FTPクライアントはWinSCPが好きだったのですが、Macでは使えないのでFileZillaを使っています。
最近はWebサイトのサーバーにファイル転送をする程度なので、どっちでも問題ない感じでした。

TeX環境

Windowsはインストーラーから入れて、Macは brew install texlive から入れられたのでとても簡単でした。
文書編集にはTeXstudioを使っていたのですが、Visual Studio Codeに変えました。
一時期はOverleafをよく使っていたのですが、最近はファイル同期で問題も起きなくなったので、ブラウザよりもかなりレスポンスの早いローカル環境の方が快適になりました。
Visual Studio Code内で使うLaTeX Workshopの設定は、以下のサイトを参考にしました。
LaTeX環境構築の手順

テキストエディタ

EmEditorは巨大ファイルの読み込みが軽く、Windowsではこちらを使っています。
MacではCotEditorを使っています。

メーラー

Windows版とMac版があって、両方同じ感覚で使えました、ありがたいです。
メールのスレッド表示が読みやすいのと、アーカイブ機能が使えるのが良いです。
その昔Thunderbirdが徐々に重くなっていき使いにくくなりeM Clientに乗り換え、それからずっと愛用しています。
機能制限ありの無料版はメールアカウントが2つまで設定でき、カレンダーも使えます。
機能制限がない無期限アップグレード可能な永久ライセンスもあり、端末3つまでアクティベートできます（自分はこちらを使ってます）。

アーカイバー

圧縮ファイルの中身を確認してから解凍したいので、ファイルマネージャー付きのソフトウェアを使っています。
比較的似たようなことができます。

Word/Excel/PowerPoint

マクロとかを使わない限り、ほぼ気にならないレベルで互換性があると感じます。
論文を書く、表を作る、発表資料作成ぐらいの用途であれば、おおむね支障はないと思います。

画面共有

Mac→Macは画面共有アプリから vnc://IPアドレス で接続するのみで、快適に使えます。
Windows→Macは今のところrealVNCを使っています（ややラグいです）。

ファイルマネージャー

  • Windows 11でExplorerにタブが追加され、使いやすくなりました。
  • MacのFinderもCommand+Tでタブを追加でき、これも使いやすいです。

Mac用の設定

Windows 3.1からずっと使っていてWindowsに慣れきっているため、MacをWindowsに寄せる設定をいくつか行いました。

  • キー配置変更

ウィンドウ切り替え（次のウィンドウを操作対象にする）: 一画面派なので、キーボード設定からOption+Tabにする（Command+Tabの隣に配置）。
MacBook: Caps LockをCommandに変更
MacStudio: Caps LockをControlに、ControlをCommandに、WindowsキーをOptionに変更、半角/全角キーが機能するように変更（Macで半角/全角キーによる快適な入力切り替えを実現する方法 | しあ）。これはWindows用のキーボードを使っているため必要な設定で、普通にMac用のキーボードを買ったほうが絶対よい。

  • マウス・トラックパッド

LinearMouse でスクロール加速をなくしました（スクロールモードを行に変更すると加速度なしです）。
Scroll Reverser でマウスは逆に、トラックパッドはそのままにしました。

  • 日本語入力

日本語変換の挙動をWindowsに近づけるため、Google日本語入力 を入れました。
インストールするだけでは自動でGoogle日本語入力を使う設定にはならないので、システム設定のキーボード設定で入力ソースを追加する作業が必須でした。
Macでの日本語入力を快適にするGoogle日本語入力 | 社会保険労務士北村事務所

  • バックスラッシュと¥記号

Google日本語入力の方は、ツールバーアイコンから環境設定を開いて簡単に設定できました。
英字入力は、Macのデフォルト日本語入力の入力モードに英字を追加し、¥を入力するとバックスラッシュが入力できるように設定できます。
ややこしいですが、↑を設定後に、デフォルト英字入力（ABCで表示されているもの）を入力ソースから削除すると、Google日本語入力とバックスラッシュ設定済みの英字入力が使える様になりました。
デフォルトの"¥"キー入力を「\」に変更する方法 | NonNonLife

  • Fnキー

Windowsでは、Fn6〜Fn10にひらがな、半角英数等が割り当てられていてよく使うのですが、Macでは挙動が違うのでこれも変更しました。
MacとWindowsでは逆になっていて、キーボード設定から変えるだけでした。
MacでファンクションキーのF7・F8・F9・F10で変換したい | Office01

  • 仮想マシン

Parallels Desktop 26 の永久ライセンスを買ってみました。
Windows Proまでは不要だったので、やや割安なWindows Homeを買いました。
Windows Proが推奨されているそうですが、Windows Homeも使えるようでしたので、以下の手順でインストールしました。
Parallels Desktopをインストールすると最初に勝手にWindows Proがインストールされるので、完了するまで待ちます。
インストール完了後に、Parallels Desktopのコントロールセンターを開いて+マークを押すと、ISOイメージからOSをインストールするメニューが選択できます。
Windows 11 Arm64版のISOイメージ をMicrosoftのサイトからダウンロードし、これを指定するとWindows Home版のインストールを選択できます。
Windows Homeのインストールが完了したら、コントロールセンターでWindows Proを削除して完了でした。

A one sample test for the median survival time

統計

I couldn't find a concise summary of the test for the one sample test for median survival time (MST), so I'm writing this post.

Hypohteses

Let's consider a hypothesis testing problem of the MST, $M$, is equals to $M_0$:
\[
\left\{
\begin{align*}
H_0: M = M_0 \\
H_1: M \ne M_0 \\
\end{align*}
\right. \,.
\]

Test statistic

This hypothesis testing problem has been proven to be equivalent to a one-sample test based on the Kaplan–Meier estimator, $\hat{S}(t)$ (Brookmeyer and Crowley. Biometrics, 1982). The above hypothesis can be tested with the following test statistic:
\[
Z = \frac{g\{\hat{S}(M_0)\} - g\{S_0(M_0)\}}{\hat{\tau}(M_0)/\sqrt{n}} \to_L N(0,1) \; \mathrm{under} \; H_0,
\] where the function $g$ is a transformation for improving the normal approximation in small-sample settings, $\hat{\tau}^2(t)$ is an asymptotic variance estimator of the transformed Kaplan–Meier estimator at $t$, and $S_0$ is the survival function of the null hypothesis and $S_0(M_0) = 0.5$ holds. Thus, this statistic is equivalent to
\[
Z = \frac{g\{\hat{S}(M_0)\} - g(0.5)}{\hat{\tau}(M_0)/\sqrt{n}}.
\] For $g$, the arcsine square-root or log-minus-log transformations (i.e., $\arcsin\sqrt{S(t)}$ or $\log\{-\log S(t)\}$) are known to be appropriate (Borgan and Listøl. Scandinavian Journal of Statistics, 1990). The variance estimator for $g$ is given by the following formula:
\[
\hat{\tau}^2(t)=n g'\{\hat{S}(t)\}^2 \{\hat{S}(t)\}^2 \int_0^t \frac{1}{Y(s)\{Y(s)-\Delta N(s)\}} \mathrm{d}N(s),
\] where $N(t)=\sum_{i=1}^{n} N_i(t)$, $Y(t)=\sum_{i=1}^{n} Y_i(t)$, $\Delta N(t) = N(t) - N(t-)$, and $N_i(t)$ and $Y_i(t)$ is the counting process and at risk process for subject $i$ at $t$, respectively. For the arcsine square-root transformation,
\[
g'\{S(t)\}=\{4S(t)(1-S(t))\}^{-1/2}
\] and for the log-minus-log transformation,
\[
g'\{S(t)\}=\{S(t)\log S(t)\}^{-1}
\] holds.

R code

Here is an example of R code for the arcsin transformation ($M_0=8$).

library(survival)
library(ggfortify)
library(dplyr)

# dummy data
set.seed(135)
m1 <- 12
fol <- 12
n <- 90
t1 <- rexp(n, -log(0.5)/m1)
u1 <- runif(n, 0, fol)
f1 <- u1 + fol
delta1 <- ifelse(t1 < f1, 1, 0)
y1 <- ifelse(t1 < f1, t1, f1)

# KM estimator at M0
m0 <- 8
sfit1 <- survfit(Surv(t1, delta1) ~ 1, conf.type = "plain")
shat <- fortify(sfit1) %>% filter(time < m0) %>% slice(n())

# test statistic and p-value
gs0 <- asin(sqrt(0.5))
gs1 <- asin(sqrt(shat$surv))
se <- shat$std.err*sqrt(shat$surv/(1 - shat$surv)*0.25)
z <- (gs1 - gs0)/se
pval <- 2*(1 - pnorm(abs(z)))

median(sfit1)
shat$surv
z
pval

> median(sfit1)
      50 
12.34418 
> shat$surv
[1] 0.6666667
> z
[1] 3.223976
> pval
[1] 0.00126424

References

  • Brookmeyer R, Crowley J. A confidence interval for the median survival time. Biometrics 1982; 38(1): 9–41.
  • Borgan Ø, Liestøl K. A note on confidence intervals and bands for the survival function based on transformations. Scandinavian Journal of Statistics 1990; 17(1): 35–41.

I referred to this blog post for calculating variance:
dominicmagirr.github.io

Fine–Grayモデルのscaled Schoenfeld residualsをRで計算してみよう

R 統計

はじめに

RでFine–Grayモデルのscaled Schoenfeld residualsをRで計算したかったのですが、cmprskパッケージではSchoenfeld residualsしか計算することができませんでした。
今後もたまに使いそうなので、scaled Schoenfeld residualsを計算した結果をメモしておくことにしました。

scaled Schoenfeld residuals

Schoenfeld residualsはCox回帰モデルの比例ハザード性の評価によく用いられる残差（対数部分尤度の偏導関数に対する各被験者の寄与に基づいた残差）で、一般的な統計ソフトウェアで簡単に計算することができます。
Scaled Schoenfeld residualsはGrambsch & Therneau (1994)によって提案された残差で、Schoenfeld residualsの期待値が分散に依存しないように調整したもので、通常のSchoenfeld residualsよりも診断力が良いと言われています。
できれば、scaled Schoenfeld residualsの方を使いたいのですが、Rのcmprskパッケージのcrr関数では、これを計算できないようでした。

Cox回帰モデルのSchoenfeld residuals, $r_k(\boldsymbol{\beta})$, はイベント時間$t_k$（$k=1,2,\ldots,d$）に対して
\[
\mathbf{r}_k(\boldsymbol{\beta})=\mathbf{Z}_k(t_k) - \mathbf{M}(\boldsymbol{\beta}, t_k),
\]

\[
\mathbf{M}(\boldsymbol{\beta}, t_k)=\mathbf{S}^{(1)}(\boldsymbol{\beta}, t)/S^{(0)}(\boldsymbol{\beta}, t),
\]

\[
\mathbf{V}(\boldsymbol{\beta}, t_k)=
\frac{\mathbf{S}^{(2)}(\boldsymbol{\beta}, t)}{S^{(0)}(\boldsymbol{\beta}, t)}-
\left\{ \frac{\mathbf{S}^{(1)}(\boldsymbol{\beta}, t)}{S^{(0)}(\boldsymbol{\beta}, t)} \right\}^{\otimes 2},
\]

\[
\mathbf{S}^{(r)}(\boldsymbol{\beta}, t)=
\sum_{i=1}^n Y_i(t) \mathbf{Z}_i(t)^{\otimes r} \exp\{\boldsymbol{\beta}' \mathbf{Z}_i(t) \},
\]

で定義されています。
ただし、$\mathbf{Z}_i(t)$は個体$i$の時間$t$における共変量ベクトル、$\boldsymbol{\beta}$は$p$次元の回帰パラメータベクトル、$Y_i(t)$は個体$i$の時間$t$における非イベントを表わす$\{0,1\}$をとる確率過程とします。
残差の総和は$0$となり、非イベント例の残差は必ず$0$になるので、イベント時間のみを考えれば良いタイプの残差です。
Schoenfeld residualsの推定値は、回帰パラメータの推定値を代入することにより

\[
\hat{\mathbf{r}}_k(\hat{\boldsymbol{\beta}})=\mathbf{Z}_k(t_k) - \mathbf{M}(\hat{\boldsymbol{\beta}}, t_k),
\]

として得られます。

Scaled Schoenfeld residualsは、これを$\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}, t_k)$の逆行列でスケーリングした

\[
\hat{\mathbf{r}}_k^*(\hat{\boldsymbol{\beta}})=\{\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}, t_k) \}^{-1} \hat{\mathbf{r}}_k(\hat{\boldsymbol{\beta}}),
\]

で定義されています。
$\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}, t_k)$は時間の関数でまあまあ計算が大変ですが、時間変化に対する変動があまり大きくないことが知られており、Grambsch & Therneau (1994) で近似式

\[
\{\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}, t_k) \}^{-1} \approx d \widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}),
\]

が同時に提案されてます。
大体のソフトウェアはこれを使って近似したscaled Schoenfeld residuals

\[
\tilde{\mathbf{r}}_k^*(\hat{\boldsymbol{\beta}}) = d \widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}) \hat{\mathbf{r}}_k(\hat{\boldsymbol{\beta}}),
\]

を出しているようです。
スケールを適宜変換したイベント時間$t_k$（変換なし、対数、Kaplan–Meier推定値、ランク変換などがよく使われる）と$\tilde{\mathbf{r}}_k^*(\hat{\boldsymbol{\beta}}) + \hat{\boldsymbol{\beta}}$の散布図を書いてスムージングしたものを使用して、比例ハザード性の評価を行うことができます（スコア検定も一緒に提案されています）。
スムージングしたプロットが直線的であれば概ね比例ハザード性が成立していて、トレンドがある場合に比例ハザード性からの乖離が疑われます。
「生存時間解析入門 原書第2版」の6章にも、もう少し詳しい説明が記載してあります。

Fine–Grayモデルのscaled Schoenfeld residuals

Fine–Grayモデルも基本的な性質はCox回帰モデルと同じであり、同様の計算でscaled Schoenfeld residualsが計算可能と考えられます。
Rのcmprskパッケージのcrr関数ではSchoenfeld residualsまでは計算してくれるので、これを少し変換すれば

\[
\tilde{\mathbf{r}}_k^*(\hat{\boldsymbol{\beta}}) = d \widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}) \hat{\mathbf{r}}_k(\hat{\boldsymbol{\beta}}),
\]

が容易に計算できると思われました。

Rのコード

cmprskパッケージのダミーデータを使い、変数x1のscaled Schoenfeld residualsのプロットを書いたコードです（コードがミスっていたので修正しました）。

library(cmprsk)
library(splines)
set.seed(10)
ftime <- rexp(200)
fstatus <- sample(0:2, 200, replace = TRUE)
cov <- matrix(runif(600), nrow = 200)
dimnames(cov)[[2]] <- c('x1','x2','x3')

p <- dim(cov)[2]
nevent <- sum(fstatus == 1)
tk <- ftime[fstatus == 1]
z <- crr(ftime, fstatus, cov)
# z$resがSchoenfeld residuals
sz <- nevent*t(z$var %*% t(z$res)) + matrix(rep(z$coef, nevent), ncol = p, byrow = TRUE)
sz <- data.frame(tk = tk[order(tk)], sz)

# splineをしつつ、散布図を描く
sfit <- lm(X1 ~ ns(tk, df = 3), data = sz)
x <- seq(0, 5, length.out = 200)
plot(sz$tk, sz$X1, cex = 1, pch = 16, xlab = "Time", ylab = "beta(t) for x1")
lines(x, predict(sfit, data.frame(tk = x)), lwd = 2)
abline(a = 0, b = 0, lty = 2)

タイがあるときはz$resの値が少し怪しく（平均化される？）自力で計算が必要になりそうです。

参考文献

GitHub Pagesとpkgdownの使い方を勉強したまとめ

R

パッケージ紹介ページを作成したかった

最近Rパッケージを作り、せっかくだからpkgdownのページを作ってみたいなと思ってちょっと試した結果をメモしておきます。

パッケージを作るときの（追加）作業

RStudioのプロジェクト機能を使用して、GitHubと連携してパッケージを作成済みであることを想定します。
happygitwithr.com


Reference manualのファイル（Rdファイル）はroxygen2の機能を使用して作成しておきます（Rのソースファイルにコメントとしてreference manualファイルの生成に必要な内容を所定のタグを使って定義する）。
roxygen2.r-lib.org


pkgdownにはreference manualファイルの内容を見栄えの良いページに変換してくれるので、とてもありがたいです。
pkgdownでは 色々できる ようです（他の機能も試して、感想などをまとめたいところです）。
pkgdown.r-lib.org

とりあえずreference manualファイルとindex.mdまたはREADME.mdがあればpkgdownでページを作成できるようだったので、README.mdも書いておきます。
rmarkdown.rstudio.com

pkgdownの設定と実行

以下のページに作業手順が書いてあります。
pkgdown.r-lib.org

RStudioでパッケージのプロジェクト内で、初回のみ以下を実行するようです。

# Run this once to publish your site regularly
usethis::use_pkgdown_github_pages()
# Run once
# Remove docs/ from gitignore to ensure it is checked into git.
usethis::use_pkgdown()

以下と全く同じエラーが出たので、git remote add origin [URL] をしておきました。
forum.posit.co

Reference manualファイルやコンテンツをアップデートしたときは、以下で再生成できます。

# Run everytime you want to update your site
pkgdown::build_site()

生成したHTMLファイル等はパッケージディレクトリの /docs に格納されます。
デフォルト設定では .gitignore に /docs が追記され、gh-pages というbranchが作られます。
GitHub Pagesでは gh-pages という名前のbranchをドキュメント公開元にする仕様だったそうですが、今はオプションから別のbranchを使う設定ができるようになっています。
切り替えがやや面倒くさいので、自分はmaster branch（以下ではmainという名前のbranch）から公開する設定を試しました。
.gitignore から /docs を消し、生成された /docs のファイルは適当なタイミングでcommitしてpushしておきます。

GitHub Pagesの設定

パッケージを載せているrepositoryのPagesの設定を行います。
Settingsから

Pagesを開き

SourceとBranchの所を適切に設定します。

gh-pages のbranchを使う場合は、↑の画像とは異なる設定にする必要があります（これは main という名前のbranchから公開する設定です）。

ページ公開ができていればチェックマークがつくようです。


今回はGitHub Pagesの機能を使ってトップページも一緒に作りました。
nshi-stat.github.io

数式を使用したい場合

公式ページには_pkgdown.ymlのtemplateの所にmath-rendering: katexを追記して、pkgdown::build_site()を再実行するとKaTeXを使用して数式をレンダリングしてくれると書いてありました。

template:
  bootstrap: 5
  math-rendering: katex

しかし、私の環境ではこれは動作せず、以下のようにCDNから読み込むように直接書く必要がありました（将来的には修正されるかもしれません）。

template:
  bootstrap: 5
  includes:
    in_header: |
      <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.16.21/dist/katex.min.css" integrity="sha384-zh0CIslj+VczCZtlzBcjt5ppRcsAmDnRem7ESsYwWwg3m/OaJ2l4x7YBZl9Kxxib" crossorigin="anonymous">
      <script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.21/dist/katex.min.js" integrity="sha384-Rma6DA2IPUwhNxmrB/7S3Tno0YY7sFu9WSYMCuulLhIqYSGZ2gKCJWIqhBWqMQfh" crossorigin="anonymous"></script>
      <script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.21/dist/contrib/auto-render.min.js" integrity="sha384-hCXGrW6PitJEwbkoStFjeJxv+fSOOQKOPbJxSfM6G5sWZjAyWhXiTIIAmQqnlLlh" crossorigin="anonymous" onload="renderMathInElement(document.body);"></script>

KaTeXのCDNを読み込むコードはバージョンが上がると変わるので、ちょっと面倒です。
最新版 を確認して使うと良いのかなと思います。

GitHub Pagesのサイト公開の設定

GitHubに username.github.io という名前のrepositoryを作って↑の公開設定を行うと、https://username.github.io/ のURLから見れるページが作成できました。
Repository username.github.io のサブディレクトリ AAA にはファイルを配置しても意味はなく、https://username.github.io/AAA/ からはアクセスできませんでした。
Repository username.github.io のサブディレクトリ docs/AAA にファイルを配置すると、https://username.github.io/docs/AAA/ からアクセスできるようになります。
他の名前の BBB repository（パッケージのrepositoryなど）で↑の公開設定を行うと、https://username.github.io/BBB/ のURLから見れるページが作成できました。

おわりに

RStudioとGitHubさえ使えれば、Rパッケージの開発とパッケージWebサイトの公開まで超簡単にできる時代になり、隔世の感があります。
個人のHPの方はHTMLを手書きしておりますので、これからも古式ゆかしい文化を守っていく所存でございます。
ご降誕祭、おめでとうございます。

ガンマ分布のハザード関数が収束する話

統計

ガンマ分布のハザード関数

先日調べ物をしていたときに、ガンマ分布のハザード関数が$t \to \infty$で定数$\lambda$ (rate parameter) に収束するという話を知りました。


一応簡単に導出をしているものがないか探してみたのですが、見つからなかったので確認してみることにしました。

計算してみよう

ガンマ分布のハザード関数は、
\[
h(t)=
\frac{f(t)}{1-F(t)}=
\frac{\frac{1}{\Gamma (k)} \lambda^{k} t^{{k}-1} \exp(-\lambda t) }{1 - \frac{ \gamma \left(k, \lambda t \right)}{\Gamma (k)}}=
\frac{\lambda^{k} t^{{k}-1} \exp(-\lambda t)}{\Gamma \left(k, \lambda t \right)},
\] です。
ただし、$f(t)$はガンマ分布の確率密度関数、$F(t)$はガンマ分布の累積分布関数、$\Gamma (k)=\int_0^{\infty} x^{{k}-1}\exp(-x) \mathrm{d}x$、$\Gamma (k, z)=\int_z^{\infty} x^{{k}-1}\exp(-x) \mathrm{d}x$、$\gamma (k,z)=\int_0^{z} x^{{k}-1}\exp(-x) \mathrm{d}x$ です。
あとは、これについて$t \to \infty$を考えれば良いだけです。
ここで、
\[
\lim_{t \to \infty} h(t) =
\lim_{t \to \infty} \frac{\lambda^{k} t^{{k}-1} \exp(-\lambda t)}{ \int_{\lambda t}^{\infty} x^{{k}-1}\exp(-x) \mathrm{d}x } =
\frac{0}{0},
\] となるので、ロピタルの定理を用いると（一応、$y=1/t$に置換して、前提条件が満たされることは確認しました）、
\[
\begin{aligned}
\lim_{t \to \infty} \frac{\lambda^{k} t^{{k}-1} \exp(-\lambda t)}{ \int_{\lambda t}^{\infty} x^{{k}-1}\exp(-x) \mathrm{d}x }
&=
\lim_{t \to \infty} \frac{\frac{\mathrm{d}}{\mathrm{d}t}\left\{ \lambda^{k} t^{{k}-1} \exp(-\lambda t) \right\}}{\frac{\mathrm{d}}{\mathrm{d}t} \left\{\int_{\lambda t}^{\infty} x^{{k}-1}\exp(-x) \mathrm{d}x \right\}}
\\ &=
\lim_{t \to \infty} \frac{\lambda^k ({k}-1)t^{{k}-2} \exp(-\lambda t) - \lambda^{{k}+1} t^{{k}-1} \exp(-\lambda t) }{0 - \lambda \{ (\lambda t)^{{k}-1} \exp(-\lambda t) \}}
\\ &=
\lim_{t \to \infty} \lambda + \frac{1-k}{t}
\\ &=
\lambda,
\end{aligned}
\] が得られます（定積分の微分は普段使わないので良い復習になりました）。

プロットを確認

Rでは、expintパッケージのgammainc関数を用いれば、不完全ガンマ関数の計算ができるようです。
$k$が大きいと収束が遅くなるので、小さい値を設定してグラフを書いてみます。

k <- 1.05
lambda <- 2
t <- seq(0.01, 20, by = 0.01)
y <- lambda^k*t^(k - 1)*exp(-lambda*t)/expint::gammainc(k, lambda*t)
plot(t, y)
abline(h = lambda)

こんな感じになりました。

補足

江村先生から同様の証明が書いてある論文があると教えていただきました（江村先生、ありがとうございました！）。

  • 武冨奈菜美, 山本和嬉. 生存時間解析・信頼性解析のための統計モデル. 日本統計学会誌 2023; 52(2): 69–112. DOI: 10.11329/jjssj.52.69.

survivalパッケージのtt関数

R 統計

はじめに

ちょっと色々とやっていて R の survival パッケージの tt 関数の挙動が気になるという話になりました。
survival パッケージでは時間変動共変量の解析が実施可能であり、時間変動共変量と時間変動効果（回帰係数）についての vignette には詳細な解説が書かれています。
vignette 中のコードには色々な事例が書かれているのですが、tt 関数を使っているものがあります。
tt 関数は、時間変換関数 を組み込むための機能です。
なんでこんなものが必要になるかというと、モデル中に共変量と時間の交互作用項（例えば、$\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ のようなもの）が含まれる場合には、部分尤度関数を計算する際にリスク集合毎に分母も含めてイチから全計算をする必要があり、単純な Cox 回帰モデルの場合と計算方法が変わるためです（おまけでも触れていますが、時間交互作用項がある場合でもデータをイベント時間毎に層を持つ形に再構成して、層化統合解析をするのと同じ計算になります）。
モデル中に共変量と時間の交互作用項がないときは、分母の計算にはリスク集合に含まれる部分までの和 $\sum_{j \in R(t_i)} \exp(\mathbf{X}_j^T \boldsymbol{\beta})$ が必要なのでソートして累積和を1セット計算するだけでよいのですが、交互作用項がある場合は $\sum_{j \in R(t_i)} \exp\{\mathbf{X}_j(t_i)^T \boldsymbol{\beta}\}$ になって毎回全計算が必要です（これを自動化するために tt 関数が用意されているということかと思います）。
この違いは、自力で Cox 回帰の部分尤度関数の計算を実装してみるとよく分かるので、やってみると良いでしょう（簡単です）。

この tt 関数は、時間関数形を

coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung,
     tt = function(x, t, ...) pspline(x + t/365.25))

こんな感じに定義して使うのですが、関数の定義を書かずに

coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung)

のようにしても結果が出てきます。
このように定義を書かない場合には一体何が出てきているのか、という話題になり、これを確認した結果をまとめてみようかなと思いました。

tt 関数のデフォルトの挙動

これについてはヘルプファイルには書いていないです（確認した限り、そのような記載は見つけられませんでした）。
正解は vignette をよく読むと書いてありました。

There are other interesting uses for the time-transform capability. One example is O'Brien's logit-rank test procedure [6]. He proposed replacing the covariate at each event time with a logit transform of its ranks. This removes the influence of any outliers in the predictor x. For this case we ignore the event time argument and concentrate on the groupings to create the following tt function.

> function(x, t, riskset, weights){
obrien <- function(x) {
r <- rank(x)
(r-.5)/(.5+length(r)-r)
}
unlist(tapply(x, riskset, obrien))
}

This relies on the fact that the input arguments to tt() are ordered by the event number or risk set. This function is used as a default if no tt argument is present in the coxph call, but there are tt terms in the model formula. (Doing so allowed me to depreciate the survobrien function).

ざっくりいうと「tt 関数は時間関数形の指定とは別の面白い使い方ができる、例えば、O'Brien's logit-rank testの実装にも使えて、これは coxph の tt 引数に何も指定しない場合の tt 関数のデフォルトになっている（これにより、私は survobrien 関数を非推奨にできた）」という感じです。
どうやら、tt 引数の指定をしないで実行すると、リスク集合毎に共変量の順位を変換した変数を生成するため、通常よくやられているような時間関数形を指定した解析とは目的が違う結果が出てくるということが分かります（Therneau 先生...、これはちょっとやめてほしいかも...）。
これで計算される回帰係数は共変量の順位に基づいた結果になっていて時間は無視されるので、時間関数形については評価できていません！

tt (time-transform) 関数のデフォルトが、全然 time-transform と関係がない仕様というのは若干思うところがあります。

通常の tt 関数の使い方

これも念のため紹介しておきます。
一番よくある使い方は、比例ハザード性の検定だったり、比較的簡単な構造を持つ時間変動効果のモデリングではないかと思います。
例えば、対数時間のスケールにおける時間変動効果 $\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ を持つモデルは tt 関数を使って当てはめることが可能です。
また、これと同じモデルにおいて、時間変動する部分の回帰係数のスコア検定を実施できる cox.zph という関数もあります。

fit1 <- coxph(Surv(time, status) ~ karno + tt(karno), data = veteran, tt = function(x, t, ...) x*log(t))
fit2 <- coxph(Surv(time, status) ~ karno, data = veteran)
summary(fit1)
cox.zph(fit2, transform = "log")

両者の実行結果は、

> summary(fit1)
Call:
coxph(formula = Surv(time, status) ~ karno + tt(karno), data = veteran, 
    tt = function(x, t, ...) x * log(t))

  n= 137, number of events= 128 

               coef exp(coef)  se(coef)      z Pr(>|z|)    
karno     -0.083723  0.919686  0.016783 -4.989 6.08e-07 ***
tt(karno)  0.013408  1.013498  0.004196  3.195   0.0014 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
karno        0.9197     1.0873    0.8899    0.9504
tt(karno)    1.0135     0.9867    1.0052    1.0219

Concordance= 0.709  (se = 0.024 )
Likelihood ratio test= 52.99  on 2 df,   p=3e-12
Wald test            = 49.81  on 2 df,   p=2e-11
Score (logrank) test = 56.97  on 2 df,   p=4e-13

> cox.zph(fit2, transform = "log")
       chisq df      p
karno   10.5  1 0.0012
GLOBAL  10.5  1 0.0012

となりました。
tt 関数を使った検定は Wald 検定（$P = 0.0014$）で、zph 検定はスコア検定（$P = 0.0012$）のため若干値は異なります。

参考文献

  • Therneau T, Crowson C, Atkinson E. Using time dependent covariates and time dependent coefficients in the Cox model. [Link]
  • O'Brien P. A nonparametric test for association with censored data. Biometrics 1978;34:243–250. DOI: 10.2307/2530014.

おまけ

一応 tt 関数の中身を見てみて、何を計算しているのかを確認しました。

library(survival)

temp <- subset(pbc, id <= 312, select = c(id:sex, stage))
pbc2 <- tmerge(temp, temp, id = id, death = event(time, status))
pbc2 <- tmerge(pbc2, pbcseq, id = id, lbili = tdc(day, log(bili)))
cox_model0 <- coxph(Surv(tstart, tstop, death == 2) ~ lbili, pbc2)
cox_model1 <- coxph(Surv(tstart, tstop, death == 2) ~ tt(lbili), pbc2)

tt <- function(x, riskset) {
  obrien <- function(x) {
    r <- rank(x)
    (r - 0.5)/(0.5 + length(r) - r)
  }
  unlist(tapply(x, riskset, obrien))
}
Y <- cbind(pbc2$tstart, pbc2$tstop, as.numeric(pbc2$death == 2))
sort.end <- order(-Y[, 2], Y[, 3])
sort.start <- order(-Y[, 1])
newstrat <- c(1L, rep(0, nrow(Y) - 1))
counts <- .Call(survival:::Ccoxcount2, Y, as.integer(sort.start - 1L), as.integer(sort.end - 1L), as.integer(newstrat))
tindex <- counts$index
istrat <- rep(1:length(counts$nrisk), counts$nrisk)
Y <- Surv(rep(counts$time, counts$nrisk), counts$status)
tv <- pbc2$lbili[tindex]
ttv <- tt(tv, istrat)
pbc3 <- data.frame(time = Y[, 1], status = counts$status, ttv, istrat)
cox_modelttv <- coxph(Surv(time, status) ~ ttv + strata(istrat), data = pbc3)

cox_model1
cox_modelttv

これが tt 関数をそのまま実行した結果で、

> cox_model1
Call:
coxph(formula = Surv(tstart, tstop, death == 2) ~ tt(lbili), 
    data = pbc2)

               coef exp(coef)  se(coef)    z      p
tt(lbili) 0.0073115 1.0073383 0.0005712 12.8 <2e-16

Likelihood ratio test=74.29  on 1 df, p=< 2.2e-16
n= 1807, number of events= 125

こっちは tt 関数が作っているデータを確認するために検算した方です。

> cox_modelttv
Call:
coxph(formula = Surv(time, status) ~ ttv + strata(istrat), data = pbc3)

         coef exp(coef)  se(coef)    z      p
ttv 0.0073115 1.0073383 0.0005712 12.8 <2e-16

Likelihood ratio test=74.29  on 1 df, p=< 2.2e-16
n= 24422, number of events= 125

データ数以外は一致していることが分かります。
検算用のデータ（pbc3）を見ると、リスク集合毎に層が作られていることが確認できます。

> pbc3
    time status         ttv istrat
1   4191      0  1.36363636      1
2   4191      0  0.36842105      1
3   4191      0  0.23809524      1
4   4191      0  2.71428571      1
5   4191      0  0.08333333      1
6   4191      0  0.62500000      1
7   4191      0  0.62500000      1
8   4191      0  1.88888889      1
9   4191      0  4.20000000      1
10  4191      0  1.00000000      1
11  4191      0  0.08333333      1
12  4191      0 25.00000000      1
13  4191      1  7.66666667      1
14  4079      0  0.78947368      2
15  4079      0  0.13333333      2
16  4079      0  0.54545455      2
17  4079      0  3.85714286      2
18  4079      0  0.30769231      2
19  4079      0  1.00000000      2
20  4079      0  1.61538462      2
21  4079      0  2.09090909      2
22  4079      0  5.80000000      2
23  4079      0  1.26666667      2
24  4079      0  0.30769231      2
25  4079      0 10.33333333      2
26  4079      0  0.54545455      2
27  4079      0  0.03030303      2
28  4079      0  0.13333333      2
29  4079      0  2.77777778      2
30  4079      1 33.00000000      2
31  3853      0  0.61290323      3

『宇宙怪人しまりす統計よりも重要なことを学ぶ』を読みました

書籍 統計

楽しく読ませていただきました。

いくつか印象に残った点についての感想です。

第1話（としまりすコラム1）、報告時に探索的研究と検証的研究を区別しないといけない、という話題が出てきていました。やっぱり、これは大事だなと思いました。完全に切り分けるのも難しいのかもしれませんが、論文に読者の判断の助けになるような記号をつけられるような仕組みがほしいと思いました（P値を使うことを禁止しないんだったらこれぐらいはやってもいいのでは）。事前に主要評価項目が定義されている、サンプルサイズ設計がちゃんとされている、解析計画が事前に公表されている等。某系列のようにsupplementaryにチェックリストへの報告義務を課すところもありますが、これはどのぐらい読まれているんだろうといつも思います。

第2話、裁判記録のお話は聞いたことがなかったので個人的に衝撃を受けました、みんな読んだ方が良いなと思いました。誠実であり続けることを誓います。

第3話、合流の話も出てきて非常に参考になりました。

第4話、層別解析とサブグループ解析の話はあるあるですね。自分も授業では触れるようにしています。確かに層化統合解析とか別の呼び方にしたほうが間違われにくそうです。ターゲット集団の話が出てきて、授業で話すんだったらどこまで話すのがいいんだろう、とふと思いました。

第5話、審査会でしまりすくんと先生はさらっと正しいコメントしていますね（カッコいい）。現状では同様の指摘に該当する研究がどのぐらいあるんだろうかと思います。このようなコメントが来ないことが常識になるような世界を整備できるようにしないといけないな、と思いました。

自分が関わる医学研究の範囲に限っての感想ですが、因果推論を導入するメリットは因果関係を示すことができるからではなく、研究計画の段階で効果の定義をめちゃくちゃ意識できるところかな、とより強く感じました。

追記

COI開示です。