Triad sou.

層内に相関のある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 のコメントを含めて読むと、仕組みがよくわかりました。

RcmdrPlugin.KMggplot2_0.2-3 now on CRAN

New version (v0.2-3) of 'RcmdrPlugin.KMggplot2' (an Rcmdr plug-in; a GUI for 'ggplot2') released.

NEWS

Changes in version 0.2-3 (2015-12-30)
  • New geom_stepribbon().
  • Pointwise confidence intervals of Kaplan-Meier plots with band (Thanks to Dr. M. Felix Freshwater).
  • Notched box plots (Thanks to Dr. M. Felix Freshwater).
  • Stratified (colored) box plots (Thanks to Dr. M. Felix Freshwater).
Changes in version 0.2-2 (2015-12-22)
  • New ggplot2's theme (theme_linedraw, theme_light).
  • New ggthemes's themes (theme_base, theme_par).
  • Fixed a bug was caused by new ggplot2's theme.
  • Fixed a bug related to ggthemes 3.0.0.
  • Fixed a bug related to windows fonts.
  • Fixed a bug related to file saving.
  • Added a vignette for dataset requirements to make Kaplan-Meier plot.
  • Added a vignette for extrafont.

geom_stepribbon()

The geom_stepribbon is an extension of the geom_ribbon, and is optimized for Kaplan-Meier plots with pointwise confidence intervals or a confidence band.

Usage
geom_stepribbon(mapping = NULL, data = NULL, stat = "identity",
  position = "identity", na.rm = FALSE, show.legend = NA,
  inherit.aes = TRUE, kmplot = FALSE, ...)

The additional argument

  • kmplot
    • If TRUE, missing values are replaced by the previous values. This option is needed to make Kaplan-Meier plots if the last observation has an 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.

Other arguments are the same as the geom_ribbon.

Example
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")

Pointwise confidence intervals of Kaplan-Meier plots with band


Notched box plots and stratified (colored) box plots


RcmdrPlugin.KMggplot2_0.2-1 is on CRAN now

I posted a new version of the "RcmdrPlugin.KMggplot2" package, which is an Rcmdr plug-in for a "ggplot2" GUI front-end.

NEWS

Changes in version 0.2-1 (2015-12-14)
  • Kaplan-Meier plot: fixed a bug was caused by a tie data handling. (Thanks to Dr. M. Felix Freshwater).
  • Slovenian translation by Matjaž Jeran (Thanks to Matjaž Jeran and their friends).
  • New ggplot2's theme (theme_dark).
  • New ggthemes's themes (theme_calc, theme_fivethirtyeight, theme_gdocs, theme_hc, theme_pander).
  • Preview button (Thanks to Prof. Erich Neuwirth).
  • GUI refinement (Thanks to Prof. Erich Neuwirth).
  • Updates for release of the next version ggplot2.
Preview button

The new version has a preview mode.
You can make plots step-by-step.

Slovenian translation

themes
  • theme_dark (ggplot2)


  • theme_calc (ggthemes)


  • theme_fivethirtyeight (ggthemes)


  • theme_gdocs (ggthemes)


  • theme_hc (ggthemes)


  • theme_pander (ggthemes)