Simulate a Weibull regression model

Question

I am trying to simulate a AFT Weibull model in R using a log-linear model

\begin{align} log(Y)=\beta_0 +\beta_1 X1+\beta_2 X2+c W \end{align} where $W ~$ a extreme value distribution. below is my simulation code:

install.packages("extRemes")
library(extRemes)

n<-500000    
x1<-rnorm(n,0,1)    
x2<-rnorm(n,0,1)    
error<-revd(n, loc = 0, scale = 1)    
b0<- 1    
b1<-0.5    
b2<--1.2    
c<-1    
time<-exp(b0+b1*x1+b2*x2+c*error)    
status<-rep(1,n)    
survreg(Surv(time, status==1) ~ x1+x2,dist="exponential")

####output 

Coefficients:    
(Intercept)          x1          x2 

  2.9390956   0.4931199  -1.2061369 

Scale fixed at 1 
#####

I can not understand why intercept is 2.94, not close to 1 ?

Thanks for the help from EDM. R Gumbel and extreme value functions are all for maximum extreme value distribution

\begin{align} f(x)=e^{(-x-e^{-x})} \end{align}

What we really need is a random variable from the minimum extreme value distribution with density function like this \begin{align} f(x)=e^{(x-e^{x})} \end{align}

Once I generate the error from a maximum extreme value distribution, I need to set it to -1*error.

n<-500000    
x1<-rnorm(n,0,1)    
x2<-rnorm(n,0,1)    
error<-revd(n, loc = 0, scale = 1)    
b0<- 1    
b1<-0.5    
b2<--1.2    
c<-1    
time<-exp(b0+b1*x1+b2*x2-c*error)    
status<-rep(1,n)    
survreg(Surv(time, status==1) ~ x1+x2,dist="exponential")


Coefficients:
(Intercept)          x1          x2 
   1.000095    0.500679   -1.200210 

Scale fixed at 1 

Answer 1

The problem is that this formulation of a Weibull model requires use of the minimum extreme value distribution for $W$. This page explains the difficulty of assuming that a "Gumbel distribution" or "extreme value distribution" is what you need without further specification.

For a location of 0, a scale of 1, and a shape of 0, the code for the revd() function you used, or for the similar rgev() function in the R evd package, is equivalent to -log(rexp(n)), where rexp(n) represents n random samples from a standard exponential distribution. That returns random samples from the standard maximum extreme value distribution, the form used by Wikipedia:

This article uses the Gumbel distribution to model the distribution of the maximum value. To model the minimum value, use the negative of the original values.

As exp(-x) = 1/exp(x), you can use -log(1/rexp(n)) = log(rexp(n)) (or just multiply all the values you got from the maximum extreme value distribution by $-1$) to get samples from the standard minimum extreme value distribution.

A quick illustration without covariates, intercept-only model, starting with a sample from the standard maximum extreme value distribution as you used:

set.seed(202)
maxEV1000 <- -log(rexp(1000)) ## what you used, for maximum EV
survreg(Surv(exp(maxEV1000))~1,dist="exponential")
# Call:
# survreg(formula = Surv(exp(maxEV1000)) ~ 1, dist = "exponential")
# 
# Coefficients:
# (Intercept) 
#    2.039815 
# 
# Scale fixed at 1 
# 
# Loglik(model)= -3039.8   Loglik(intercept only)= -3039.8
# n= 1000 

That gives a positive intercept similar to what you found for the excess intercept over what you had intended in your simulation. If you instead work with samples from the standard minimum extreme value distribution:

minEV1000 <- -maxEV1000      ## what you need, for minimum EV
survreg(Surv(exp(minEV1000))~1,dist="exponential")
# Call:
# survreg(formula = Surv(exp(minEV1000)) ~ 1, dist = "exponential")
# 
# Coefficients:
# (Intercept) 
# -0.01854772 
# 
# Scale fixed at 1 
# 
# Loglik(model)= -981.5   Loglik(intercept only)= -981.5
# n= 1000 

you get what you expect, an intercept near 0. If you don't specify an exponential model and allow a general Weibull fit, you get the correct scale (1.005) from this latter data sample, versus 1.653 when you try to fit the former data sample.