使用rpy2从Python调用带se.fit参数的survreg predict函数遇错求助

Got it, let's break down why you're seeing that "unavailable keyword" error for se.fit=True when using predict() on your survreg model, and how to fix it.

The Root Cause #

The predict.survreg() function from R's survival package does not support the se.fit parameter when type="quantile". If you check the official R docs for predict.survreg, you’ll see that se.fit is only valid when type="link" (the linear predictor) or type="response" (the mean of the response distribution, the default). When requesting quantile predictions, this parameter isn’t recognized.

Solutions #

Depending on what you need, here are two clear paths forward:

1. You don’t need standard errors for the quantiles

Simply remove the se.fit=True argument from your call. Your code will work perfectly like this:

# Assuming you're using rpy2 to interface with R
predict(fit_hs, type="quantile", p=0.5)

This will return the median (p=0.5) predictions without any errors.

2. You need standard errors for the quantiles

To get standard errors for quantile predictions, you’ll calculate them manually using the linear predictor’s standard errors and the delta method. Here’s a step-by-step implementation:

First, fetch the linear predictor values and their standard errors (using type="link" which supports se.fit):

# Get linear predictions and their SEs
link_results = predict(fit_hs, type="link", se_fit=True)
linear_pred = link_results[0]
se_link = link_results[1]

Next, extract the scale parameter from your survreg model (this varies slightly by distribution):

# Pull the scale parameter (sigma) from the model object
sigma = fit_hs.rx2("scale")[0]

Finally, calculate the quantile and its standard error based on your model’s distribution:

• Log-Normal models: The median (p=0.5) is exp(linear_pred) (since the 0.5 quantile of the normal distribution is 0). The standard error uses the delta method:

  import rpy2.robjects as ro
  median_pred = ro.r.exp(linear_pred)
  median_se = median_pred * se_link
  

• Weibull models: The median is exp(linear_pred + ro.r.qgumbel(0.5)) (where qgumbel(0.5) = -log(-log(0.5))). The standard error is exp(linear_pred + qgumbel_p) * se_link (the derivative of the quantile with respect to the linear predictor matches the quantile itself for this parameterization).

Full Example Code #

Here’s a complete snippet using rpy2 for a log-normal survreg model:

from rpy2.robjects import r, pandas2ri
from rpy2.robjects.packages import importr

# Initialize rpy2 and load survival package
pandas2ri.activate()
survival = importr("survival")

# Assume fit_hs is your pre-trained survreg model
# Fetch linear predictions and their SEs
link_output = survival.predict(fit_hs, type="link", se_fit=True)
linear_pred = link_output[0]
se_link = link_output[1]

# Calculate median prediction and its standard error
median_pred = r.exp(linear_pred)
median_se = median_pred * se_link

# Print results
print("Median Predictions:", median_pred)
print("Median Standard Errors:", median_se)

内容的提问来源于stack exchange，提问作者user3302302