Triad sou.

WSLをインストールしてRStudio Serverを動かしてみた

ということで、インストールしてみた。


WSLのインストールとUbuntuの入手。
管理者権限で起動したPowerShell上でWSLをインストール。

Enable-WindowsOptionalFeature -Online -FeatureName Microsoft-Windows-Subsystem-Linux

再起動してMicrosoft StoreからUbuntu 18.04 LTSをインストール。


Ubuntuのターミナルを起動&初期設定し、Rの最新版をCRANからインストールするように設定(2020/04/25現在ではまだR-3.6.3だった、2020/04/29にR-4.0.0のバイナリが配置されたので更新)。

echo -e "\n## For R package"  | sudo tee -a /etc/apt/sources.list
echo "deb https://cloud.r-project.org/bin/linux/ubuntu $(lsb_release -cs)-cran40/" | sudo tee -a /etc/apt/sources.list
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9

3行目が失敗する場合は以下もご参照のこと。


make, g++, r-baseをインストール。

sudo apt update
sudo apt install make g++ r-base


Rstudio Serverをインストール(適宜最新のファイルに置き換えると良いかも; 2020/5/22に1.3.959に更新されていた)。

sudo apt-get install gdebi-core
wget https://download2.rstudio.org/server/bionic/amd64/rstudio-server-1.3.959-amd64.deb
sudo gdebi rstudio-server-1.3.959-amd64.deb


RStudio Severを起動。

sudo rstudio-server start

ブラウザで開く。

http://localhost:8787/


WindowsのCドライブは/mnt/c/にマウントされている。

パッケージ作成環境構築

パッケージ作成環境構築のため、devtoolsなどを入れておきたかったが、WSL上ではxml2がうまくインストールできない事があるようなので、対処法をメモしておく。
ターミナル上で、libxml2-dev, libssl-dev, libcurl4-openssl-devをインストールしておく。

sudo apt install libxml2-dev libssl-dev libcurl4-openssl-dev

RStudioでdevtoolsと必要パッケージを入れる。

install.package("devtools")

devtoolsの関連パッケージのインストールに失敗する場合は、エラーメッセージ

mv: cannot move '/usr/lib/R/library/00LOCK-ps/00new/xxxx' to
'/usr/lib/R/library/xxxx': Permission denied

が出ている場合は、

Sys.setenv(R_INSTALL_STAGED = FALSE)

を設定した上で、再度インストールする。
終わったら一応設定を戻しておく。

Sys.setenv(R_INSTALL_STAGED = TRUE)

または、Ubuntuに何かのパッケージが足りない場合は、エラーメッセージを探してsudo apt install package-nameする。

追記

加賀谷先生の解説も大変参考になります。ありがとうございます。



Rをソースからビルドしてみよう

多分そのままだとエラーが出たためsource.listを編集しておく。

sudo vi /etc/apt/sources.list

ファイルを開いた後に

:%s/# deb-src/deb-src/
:wq


ビルドに必要なものを入れておく。

sudo apt update
sudo apt build-dep r-base

バージョン番号を変数に代入しておき、ソースファイルをダウンロードして展開。

export R_VERSION=4.0.0
curl -O https://cran.rstudio.com/src/base/R-4/R-${R_VERSION}.tar.gz
tar -xzvf R-${R_VERSION}.tar.gz
cd R-${R_VERSION}

ビルドとインストール。

./configure \
    --prefix=/opt/R/${R_VERSION} \
    --enable-memory-profiling \
    --enable-R-shlib \
    --with-blas \
    --with-lapack
make
sudo make install

インストールの結果確認。

/opt/R/${R_VERSION}/bin/R --version

シンボリックリンクを貼っておく。

sudo ln -s /opt/R/${R_VERSION}/bin/R /usr/local/bin/R
sudo ln -s /opt/R/${R_VERSION}/bin/Rscript /usr/local/bin/Rscript


RStudio Serverを一旦止める。

sudo rstudio-server stop

設定ファイルを編集する。

sudo vi /etc/rstudio/rserver.conf

ファイルを開いた後にiで挿入モードにし

rsession-which-r=/usr/local/lib/R

を追記、escで挿入モードを抜けて、:wqで保存して閉じる。
RStudio Serverを起動してブラウザから開く。

sudo rstudio-server start


うまくいかない時はリセットしたりして対応した。

sudo rstudio-server force-suspend-all
sudo rstudio-server stop
rm -rf /home/user/.rstudio
sudo rstudio-server start

部分集団と集団全体での割合の差に対する推測

以下のような状況を想定する。

  • あるマーカーが陽性の時の真の有病率を$p_1$、あるマーカーが陰性の時の真の有病率を$p_2$とする
  • ある集団では、陽性の$n_1$人のうち病気は$X_1$人であり、陰性の$n_2$人のうち病気は$X_2$人であった
  • ある集団全体では$n=n_1+n_2$人中$X_1+X_2$人が病気であった

一般的には、陽性の有病率$p_1$と陰性の有病率$p_2$を比較したリスク差$p_2-p_1$を考える事が多いと思われるが、ここでは陽性の有病率$p_1$と全体の有病率$p$を比較したリスク差$p-p_1$の推測について考えてみたい。

リスク差$p-p_1$に対する推測

いわゆる二項分布の仮定$X_1 \sim \mathrm{Bin}(n_1, p_1)$および$X_2 \sim \mathrm{Bin}(n_2, p_2)$のもとでは、$j=1,2$について、$\mathrm{E}[X_j]=n_jp_j$、$\mathrm{Var}[X_j]=n_jp_j(1-p_j)$である。

各有病率$p_j$の推定量は$\hat{p}_j=\frac{X_j}{n_j}$である。$\mathrm{E}[\hat{p}_j]=\frac{n_jp_j}{n_j}=p_j$であり、これは不偏推定量である。

$p$の推定量は$\hat{p}=\frac{X_1+X_2}{n}$である。陽性割合を$w=n_1/n$とおけば、$\mathrm{E}[\hat{p}]=\frac{n_1p_1+n_2p_2}{n}=wp_1+(1-w)p_2$と書ける。

そして、リスク差$p-p_1$の点推定量は、
\[
\begin{split}
\hat{p}-\hat{p}_1
&=
\frac{X_1+X_2}{n}-\frac{X_1}{n_1}
\\&=
\frac{n_2}{n}\frac{X_2}{n_2}-\frac{n_2}{n}\frac{X_1}{n_1}
\\&=
(1-w)(\hat{p}_2-\hat{p}_1),
\end{split}
\]であり、期待値は$\mathrm{E}[\hat{p}-\hat{p}_1]=(1-w)(p_2-p_1)$となる。

また、$\hat{p}-\hat{p}_1$の分散は、
\[
\begin{split}
\mathrm{Var}[\hat{p}-\hat{p}_1]
&=
\mathrm{Var}[(1-w)(\hat{p}_2-\hat{p}_1)]
\\&=
(1-w)^2(\mathrm{Var}[\hat{p}_2]+\mathrm{Var}[\hat{p}_1]),
\end{split}
\]となり、$\hat{p}-\hat{p}_1$の分散推定量は、$\widehat{\mathrm{Var}}[\hat{p}-\hat{p}_1]=(1-w)^2(\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1])$である。ただし、$\mathrm{Var}[\hat{p}_j]=\frac{p_j(1-p_j)}{n_j}$、$\widehat{\mathrm{Var}}[\hat{p}_j]=\frac{\hat{p}_j(1-\hat{p}_j)}{n_j}$である。

正規近似できる事を想定すれば、$p-p_1=0$を帰無仮説とした検定統計量
\[
\begin{split}
Z
&=
\frac{(1-w)(\hat{p}_2-\hat{p}_1)}{\sqrt{(1-w)^2(\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1])}}
\\&=
\frac{\hat{p}_2-\hat{p}_1}{\sqrt{\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1]}},
\end{split}
\]が得られるが、これは$p_2-p_1=0$を帰無仮説としたナイーブな検定統計量に完全に一致する。

リスク差$p-p_1$の信頼区間は、
\[
(1-w)\left[(\hat{p}_2-\hat{p}_1) \pm \Phi^{-1}(1-\alpha/2) \sqrt{\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1]}\right]
\]と書け、リスク差$p_2-p_1$の信頼区間を$(1-w)$倍すれば求められる事が分かる。

まとめ

リスク差$p-p_1$の点推定量区間定量はリスク差$p_2-p_1$の点推定量区間定量から計算でき、$p-p_1=0$を帰無仮説としたナイーブな検定は$p_2-p_1=0$を帰無仮説としたナイーブな検定と同一の$P$値を与える事が分かる。

補足


篠崎先生、ありがとうございました。

層内に相関のあるK×2×2分割表の共通オッズ比の推測について

Liang KY. Odds ratio inference with dependent data. Biometrika 1985; 72(3): 678–682.


と言ったしたものの確認したことがなかったのでシミュレーションしてみました。
マルコフ連鎖よりもベータ二項分布の方が簡単そうだったので、積ベータ二項分布モデルについて検討しています(原典と微妙な違いがあるかもしれません)。

シミュレーション用コード

require("VGAM")
require("reshape2")
require("foreach")
require("doParallel")
mhsim <- function(k, phi, rho, n = 5, m = 5, r = 5000) {
  p2k0 <- 0.2
  cluster <- makeCluster(detectCores() - 1)
  registerDoParallel(cluster)
  res <- foreach(i = 1:r, .combine = "rbind") %dopar% {
    tab <- array(dim = c(2, 2, k))
    p2 <- p2k0
    for (l in 1:k) {
      p2 <- p2 + 0.6/k
      p1 <- (p2*phi) / (1 - p2 + p2*phi)
      n11 <- VGAM::rbetabinom(1, n, p2, rho)
      m11 <- VGAM::rbetabinom(1, m, p1, rho)
      tab[2, 1, l] <- n11
      tab[1, 1, l] <- m11
      tab[2, 2, l] <- n - n11
      tab[1, 2, l] <- m - m11
    }
    dat1 <- reshape2::melt(tab)
    dat1 <- data.frame(x=rep(dat1$Var1 - 1, dat1$value),
                       y=rep(dat1$Var2 - 1, dat1$value),
                       k=rep(dat1$Var3, dat1$value))
    mh <- mantelhaen.test(tab)
    cmle <- mantelhaen.test(tab, exact = TRUE)
    mle <- glm(y ~ factor(x) + factor(k),
                   family = binomial(link = "logit"),
                   data = dat1)
    data.frame(mh = mh$estimate, cmle = cmle$estimate,
               mle = exp(mle$coefficients[2]))
  }
  stopCluster(cluster)
  res
}
summary(mhsim(k =  25, phi = 5, rho = 0.2))
summary(mhsim(k =  50, phi = 5, rho = 0.2))
summary(mhsim(k = 100, phi = 5, rho = 0.2))
summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100))
summary(mhsim(k = 100, phi = 5, rho = 0.0))
summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100))

相関の影響についての結果

オッズ比$\phi=5$、相関のパラメータ$\rho=0.2$、各層内のサンプルサイズが$N_k=M_k=5$、$K=25, 50, 100$のときを検討します。
mh:Mantel–Haenszel推定量、cmle:条件付き最尤推定量、mle:無条件の最尤推定量です。
Mantel–Haenszel推定量は$K$が大きければ$\phi=5$に近づいていきます。
条件付き最尤推定量は$K$が大きくなっても$\phi=5$から離れたままです。
無条件の最尤推定量はもともと一致性がありません。

> summary(mhsim(k =  25, phi = 5, rho = 0.2))
       mh              cmle              mle         
 Min.   : 1.086   Min.   :  1.098   Min.   :  1.109  
 1st Qu.: 3.833   1st Qu.:  4.402   1st Qu.:  5.236  
 Median : 5.167   Median :  6.059   Median :  7.497  
 Mean   : 6.018   Mean   :  7.074   Mean   :  9.313  
 3rd Qu.: 7.125   3rd Qu.:  8.459   3rd Qu.: 11.038  
 Max.   :89.333   Max.   :115.649   Max.   :231.215  
> summary(mhsim(k =  50, phi = 5, rho = 0.2))
       mh              cmle             mle        
 Min.   : 2.042   Min.   : 2.148   Min.   : 2.340  
 1st Qu.: 4.130   1st Qu.: 4.786   1st Qu.: 5.764  
 Median : 5.064   Median : 5.968   Median : 7.395  
 Mean   : 5.410   Mean   : 6.381   Mean   : 8.096  
 3rd Qu.: 6.261   3rd Qu.: 7.511   3rd Qu.: 9.606  
 Max.   :35.167   Max.   :32.358   Max.   :54.375
> summary(mhsim(k = 100, phi = 5, rho = 0.2))
       mh              cmle             mle        
 Min.   : 2.458   Min.   : 2.770   Min.   : 3.109  
 1st Qu.: 4.365   1st Qu.: 5.101   1st Qu.: 6.186  
 Median : 5.041   Median : 5.961   Median : 7.392  
 Mean   : 5.218   Mean   : 6.167   Mean   : 7.730  
 3rd Qu.: 5.880   3rd Qu.: 7.012   3rd Qu.: 8.899  
 Max.   :14.608   Max.   :16.252   Max.   :24.016  
> summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100))
       mh             cmle            mle        
 Min.   :3.928   Min.   :4.636   Min.   : 5.555  
 1st Qu.:4.757   1st Qu.:5.628   1st Qu.: 6.930  
 Median :5.061   Median :6.040   Median : 7.504  
 Mean   :5.098   Mean   :6.042   Mean   : 7.514  
 3rd Qu.:5.433   3rd Qu.:6.390   3rd Qu.: 7.995  
 Max.   :6.756   Max.   :7.992   Max.   :10.339  

おまけ:無相関の積ベータ二項分布モデル

相関がない積ベータ二項分布モデルは積二項分布モデルに等しく、Mantel–Haenszel推定量も条件付き最尤推定量も一致性があります。

> summary(mhsim(k = 100, phi = 5, rho = 0.0))
       mh             cmle            mle        
 Min.   :2.945   Min.   :3.050   Min.   : 3.465  
 1st Qu.:4.498   1st Qu.:4.495   1st Qu.: 5.363  
 Median :5.005   Median :5.007   Median : 6.059  
 Mean   :5.090   Mean   :5.075   Mean   : 6.169  
 3rd Qu.:5.597   3rd Qu.:5.573   3rd Qu.: 6.842  
 Max.   :9.277   Max.   :9.877   Max.   :13.188  

おまけ:無条件のMLE

よく知られた結果ですが、無条件のMLEは1:1マッチングのときに$\phi^2$に収束します(Breslow, 1981)。

> summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100))
       mh             cmle            mle       
 Min.   :3.423   Min.   :3.423   Min.   :11.72  
 1st Qu.:4.311   1st Qu.:4.310   1st Qu.:18.58  
 Median :4.872   Median :4.872   Median :23.73  
 Mean   :5.013   Mean   :5.012   Mean   :25.99  
 3rd Qu.:5.457   3rd Qu.:5.457   3rd Qu.:29.78  
 Max.   :7.966   Max.   :7.965   Max.   :63.45  
  • Breslow N. Odds ratio estimators when the data are sparse. Biometrika 1981; 68(1): 73–84.

SAS/STAT 14.3

New Procedure

14.2に続いて因果推論系のプロシジャが追加されていますね。

  • CAUSALMED Procedure: Causal mediation analysis用のプロシジャですね (The CAUSALMED Procedure)。いくつかの方法で、直接効果・間接効果を分解した推定ができるようです。

Enhancements

コンパートメントモデルやcause-specific PHモデルなど実務で使いそうなものがサポートされています。

  • FREQ Procedure: COMMONRISKDIFFオプションでoverall risk differenceの推定値、信頼区間が出力されるようになった。METHOD=SCORE2オプションが追加され、リスク差やリスク比の正確な信頼区間の計算方法が追加された (Agresti and Min 2001)。Bowker’s symmetry testのexact版が追加された。(What’s New in SAS/STAT 14.3 FREQ Procedure)。
  • NLMIXED and MCMC Procedures: CMPTMODEL statementが追加され、コンパートメントモデルの解析が簡単に?できるようになった (CMPTMODEL Statement)。
  • PHREG Procedure: Cause-specific proportional hazards modelsのサポート (PHREG Procedure MODEL Statement)。
  • TTEST Procedure: ブートストラップ法に対応 (What’s New in SAS/STAT 14.3 TTEST Procedure)。

メモ

このページ に新情報が結構まとまっています。

SAS/STAT 14.2

New Procedures

最近は因果推論関係に力をいれているのか、procedure がふたつ追加されていました What's New in SAS/STAT 14.2 - New Procedures)。

  • CAUSALTRT Procedure: 二値 or 連続アウトカム、二値の治療変数の場合の average causal effect (ATT と ATE) を推定できる (The CAUSALTRT Procedure)。
  • PSMATCH Procedure: モデルを指定して簡単にマッチングが実施できる (The PSMATCH Procedure)。Exact マッチング、Mahalanobis distance マッチング、PS マッチングができる。

Enhancements

機能追加関連も、時間依存性ROC解析やNLMIXEDのマルチスレッド化など色々とあるようでした (What's New in SAS/STAT 14.2 - Highlights of Enhancements)。

  • PHREG Procedure: 時間依存性ROC解析に対応 (SUGI 2017 の paper に解説あり; Evaluating Predictive Accuracy of Survival Models with PROC PHREG)。Concordance Statistics (Harrell's C-statistic, Uno's C-statistic)、経時AUC・各時点のROC曲線の図示。
  • NLMIXED Procedure: 変量効果を指定している場合にマルチスレッド化して処理できるように (STAT14.1までは変量効果がある場合にNTHREADSを指定しても無視されていたらしいです)。

RcmdrPlugin.KMggplot2_0.2-4 is on CRAN

RcmdrPlugin.KMggplot2 v0.2-4 (an Rcmdr plug-in; a GUI for 'ggplot2') was released on Dec. 20, 2016.

NEWS

Changes in version 0.2-4 (2016-12-20)
  • ggplot2 2.2.0 was supported.
  • Added ggplot2's colour/fill scales (scale_hue, scale_grey).
  • New plot: Beeswarm plot in Box plot (Thanks to DailyFunky @EpiFunky).
  • New function: lastcom() function to list the last command.
  • Fixed a bug: Histogram with missing values (Thanks to Dr. José G. Conde).
  • Fixed a bug: Boxplot with an inappropriate jitter position (Thanks to Dr. M. Felix Freshwater).

Beeswarm plot

f:id:triadsou:20161220235144p:plain

Example: with box plot
f:id:triadsou:20161220235149p:plain

Example: with CI
f:id:triadsou:20161220235152p:plain

ggplot2's colour/fill scales

scale_*_hue
f:id:triadsou:20161220235156p:plain

scale_*_grey
f:id:triadsou:20161220235201p:plain

ggplot2-2.0.0 の拡張

ggproto と export

ggplot2 の v2.0.0 では OO の機構が ggproto というパッケージ内に含まれた独自のものに変更され、他のパッケージからの拡張が容易になったようです。
詳しくは、Extending ggplot2 に公式の解説があるので、これを熟読すると良いでしょう。
GeomXXX, geom_xxx, StatXXX, stat_xxx などが export されたので、かなり簡単に拡張が書けるようです。

独自の geom

私も自分のパッケージで独自の geom を作ってみました。
Kaplan-Meier plot という階段状の図を作図するコンポーネントで、信頼区間も階段状に帯の様に表示したかったので、geom_ribbon を継承し、前処理を少し加えたものを作成しました。

#' Step ribbon plots.
#'
#' \code{geom_stepribbon} is an extension of the \code{geom_ribbon}, and
#' is optimized for Kaplan-Meier plots with pointwise confidence intervals
#' or a confidence band.
#'
#' @section Aesthetics:
#' \Sexpr[results=rd,stage=build]{ggplot2:::rd_aesthetics("geom", "ribbon")}
#'
#' @seealso
#'   \code{\link[ggplot2:geom_ribbon]{geom_ribbon}} \code{geom_stepribbon}
#'   inherits from \code{geom_ribbon}.
#' @inheritParams ggplot2:::geom_ribbon
#' @param kmplot If \code{TRUE}, missing values are replaced by the previous
#' values. This option is needed to make Kaplan-Meier plots if the last
#' observation has event, in which case the upper and lower values of the
#' last observation are missing. This processing is optimized for results
#' from the survfit function.
#' @examples
#' huron <- data.frame(year = 1875:1972, level = as.vector(LakeHuron))
#' h <- ggplot(huron, aes(year))
#' h + geom_stepribbon(aes(ymin = level - 1, ymax = level + 1), fill = "grey70") +
#'     geom_step(aes(y = level))
#' h + geom_ribbon(aes(ymin = level - 1, ymax = level + 1), fill = "grey70") +
#'     geom_line(aes(y = level))
#' @rdname geom_stepribbon
#' @importFrom ggplot2 layer GeomRibbon
#' @export
geom_stepribbon <- function(
  mapping = NULL, data = NULL, stat = "identity", position = "identity",
  na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, kmplot = FALSE, ...) {
  layer(
    data = data,
    mapping = mapping,
    stat = stat,
    geom = GeomStepribbon,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      na.rm = na.rm,
      kmplot = kmplot,
      ...
    )
  )
}

#' @rdname geom_stepribbon
#' @format NULL
#' @usage NULL
#' @export
GeomStepribbon <- ggproto(
  "GeomStepribbon", GeomRibbon, 
  
  extra_params = c("na.rm", "kmplot"),
  
  draw_group = function(data, panel_scales, coord, na.rm = FALSE) {
    if (na.rm) data <- data[complete.cases(data[c("x", "ymin", "ymax")]), ]
    data <- rbind(data, data)
    data <- data[order(data$x), ]
    data$x <- c(data$x[2:nrow(data)], NA)
    data <- data[complete.cases(data["x"]), ]
    GeomRibbon$draw_group(data, panel_scales, coord, na.rm = FALSE)
  },
  
  setup_data = function(data, params) {
    if (params$kmplot) {
      data <- data[order(data$PANEL, data$group, data$x), ]
      tmpmin <- tmpmax <- NA
      for (i in 1:nrow(data)) {
        if (is.na(data$ymin[i])) {
          data$ymin[i] <- tmpmin
        }
        if (is.na(data$ymax[i])) {
          data$ymax[i] <- tmpmax
        }
        tmpmin <- data$ymin[i]
        tmpmax <- data$ymax[i]
      }
    }
    data
  }
)  
geom_stepribbon

継承しているため、ほぼ geom_ribbon のコピーですが、kmplot という引数を追加しています。
これは、survival パッケージの survfit という関数を内部で利用していて、この関数の出力を使う場合の例外処理の有無を指定できるようにしたかったので入れました。

GeomStepribbon

実際に処理している本体の部分です。

  • ggproto(`_class` = NULL, `_inherit` = NULL, ...) となっていて、ひとつ目の引数がクラス名、ふたつ目が継承するオブジェクトなので、GeomRibbon を継承しています、geom を一から作成する場合は Geom を継承すれば良いと思います。
  • extra_params は追加引数名の指定です。geom_xxx の方で追加引数を入れた場合はここで指定していないとエラーになります。
  • draw_group は描画処理の部分ですが、今回はここで前処理をして GeomRibbon$draw_group を呼んでいます。
  • setup_data はデータのセットアップ処理の部分です、今回はここに例外処理を入れています。

実際には、Geom 内部で順番に継承している他の処理と上述の処理が呼ばれていくのですが、今回は一部で事足りてしまいました。
今後は今回触れていない部分を勉強してみようと思っています。

出力例

出力は、geom_stepribbon または以下のようになりました。

require("ggplot2")
data(dataKm, package = "RcmdrPlugin.KMggplot2")

.df <- na.omit(data.frame(x = dataKm$time, y = dataKm$event, z = dataKm$trt))
.df <- .df[do.call(order, .df[, c("z", "x"), drop = FALSE]), , drop = FALSE]
.fit <- survival::survfit(
  survival::Surv(time = x, event = y, type = "right") ~ z, .df)
.fit <- data.frame(x = .fit$time, y = .fit$surv, nrisk = .fit$n.risk, 
  nevent = .fit$n.event, ncensor= .fit$n.censor, upper = .fit$upper,
  lower = .fit$lower)
.df <- .df[!duplicated(.df[,c("x", "z")]), ]
.df <- .fit <- data.frame(.fit, .df[, c("z"), drop = FALSE])
.df <- .fit <- rbind(unique(data.frame(x = 0, y = 1, nrisk = NA, nevent = NA,
  ncensor = NA, upper = 1, lower = 1, .df[, c("z"), drop = FALSE])), .fit)
.cens <- subset(.fit, ncensor == 1)

ggplot(data = .fit, aes(x = x, y = y, colour = z)) + 
  RcmdrPlugin.KMggplot2::geom_stepribbon(data = .fit,
    aes(x = x, ymin = lower, ymax = upper, fill = z), alpha = 0.25,
    colour = "transparent", show.legend = FALSE, kmplot = TRUE) +
  geom_step(size = 1.5) +
  geom_linerange(data = .cens, aes(x = x, ymin = y, ymax = y + 0.02),
    size = 1.5) + 
  scale_x_continuous(breaks = seq(0, 21, by = 7), limits = c(0, 21)) + 
  scale_y_continuous(limits = c(0, 1), expand = c(0.01, 0)) + 
  scale_colour_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  xlab("Time from entry") +
  ylab("Proportion of survival") +
  labs(colour = "trt") +
  theme_bw(base_size = 14, base_family = "sans") + 
  theme(legend.position = "right")

因みに、geom_ribbon では帯は折れ線状につながるので geom_line と同じ感じです。
今回の、geom_stepribbon はある意味で geom_step と合わせて使うために作っています。
あと、geom_stepribbon の引数で colour = "transparent" とすると帯の端が透明にならないため、次のバージョンアップで colour = NA にでも変えようかと思っています。

geom-.r ファイル

また、geom-.r というファイルを roxygen2 のコメントを含めて読むと、仕組みがよくわかりました。