# Simulate a Weibull regression model

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

• Your code doesn't appear to match your model: in particular, why do you exponentiate the linear combination to create time?
– whuber
Commented Oct 11, 2022 at 17:13
• Sorry, it is actually log(Y) Commented Oct 11, 2022 at 17:15
• Please edit the question to specify which extreme value distribution parameterization you are using and the package from which you get revd(), as that function doesn't seem to be in base R. There's both a maximum and a minimum extreme value distribution; make sure you're using the correct one.
– EdM
Commented Oct 11, 2022 at 17:19
• I tried different extreme value functions or gumbel function. I have the same issue. I found out the issues with intercept only happen with extreme value distribution. If I specify W to be a normal or logistic, the intercept estimates are close to 1. Commented Oct 11, 2022 at 18:58

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.

• thanks. R functions have different definitions of extreme value distribution. As you said, they are for maximum extreme value distribution, not minimum extreme value distribution. I just can not understand why they can not standardize the definitions Commented Oct 12, 2022 at 1:47
• @Vincent it’s maybe even worse for Weibull. See this page, for example. The same words can be used to describe fundamentally different parameters, and you always have to check which parameterization is used. I still get confused every time I need to deal with Gumbel or Weibull.
– EdM
Commented Oct 12, 2022 at 2:59