survreg 関数のサンドイッチ分散

R の survival package には、パラメトリックな生存時間解析を行うための survreg 関数があります。 survreg 関数には、サンドイッチ分散 (or ロバスト分散) の推定値を計算するための robust オプションを指定することができます。 しかし、クラスタ変数がない場合に robust = TRUE を指定してしまうと、エラーが表示されてしまいます。

  require(survival)
  data(leukemia)

  leukemia$id <- factor(1:nrow(leukemia))

  result <- survreg(
    Surv(time, status) ~ x,
    data = leukemia,
    dist = "weibull", robust = TRUE
  )
  

  Error in rowsum.default(residuals.survreg(fit, "dfbeta")) : 
    argument "group" is missing, with no default
  

エラーログとソースコードから、rowsum 関数の引数が一つしかないためにエラーが表示されており、また、ヘルプファイルには、

Use robust 'sandwich' standard errors, based on independence of individuals if there is no cluster() term in the formula, based on independence of clusters if there is.

と記載されているため、survreg 関数のバグであると思われます。

クラスタ変数がない場合にサンドイッチ分散を求めたい場合には、個人を単位とするクラスタ変数を作成し、cluster の項を付け加えると良いでしょう。

  require(survival)
  data(leukemia)

  leukemia$id <- factor(1:nrow(leukemia))

  result <- survreg(
    Surv(time, status) ~ x + cluster(id),
    data = leukemia,
    dist = "weibull", robust = TRUE
  )
  summary(result)
  

  Call:
  survreg(formula = Surv(time, status) ~ x, data = leukemia, dist = "weibull")
                  Value Std. Error     z        p
  (Intercept)     4.109      0.300 13.70 9.89e-43
  xNonmaintained -0.929      0.383 -2.43 1.51e-02
  Log(scale)     -0.235      0.178 -1.32 1.88e-01

  Scale= 0.791 

  Weibull distribution
  Loglik(model)= -80.5   Loglik(intercept only)= -83.2
          Chisq= 5.31 on 1 degrees of freedom, p= 0.021 
  Number of Newton-Raphson Iterations: 5 
  n= 23
  

また、coxph で robust = TRUE を指定した場合、"一応" サンドイッチ分散が出力されますが、廃止予定のオプションですので cluster を利用するように警告が表示されます。

  require(survival)
  data(leukemia)

  result2 <- coxph(
    Surv(time, status) ~ x,
    data = leukemia,
    robust = TRUE
  )
  summary(result2)
  

  Warning message:
  In coxph(Surv(time, status) ~ x, data = leukemia, robust = TRUE) :
    The robust option is depricated

  Call:
  coxph(formula = Surv(time, status) ~ x, data = leukemia, robust = TRUE)

    n= 23, number of events= 18 

                    coef exp(coef) se(coef) robust se     z Pr(>|z|)  
  xNonmaintained 0.9155    2.4981   0.5119    0.4645 1.971   0.0487 *
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

                  exp(coef) exp(-coef) lower .95 upper .95
  xNonmaintained     2.498     0.4003     1.005     6.209

  Concordance= 0.619  (se = 0.073 )
  Rsquare= 0.137   (max possible= 0.976 )
  Likelihood ratio test= 3.38  on 1 df,   p=0.06581
  Wald test            = 3.88  on 1 df,   p=0.04873
  Score (logrank) test = 3.42  on 1 df,   p=0.06454,   Robust = 4.69  p=0.03031

    (Note: the likelihood ratio and score tests assume independence of
        observations within a cluster, the Wald and robust score tests do not).
      

履歴

• 2012/05/07 公開