[R] Linear relative rate / excess relative risk models

[Wollschlaeger, Daniel]

A while ago, I inquired about fitting excess relative risk models in R. This is a follow-up about what I ended up doing in case the question pops up again.

While I was not successful in using standard tools, switching to Bayesian modeling using rstan (mc-stan.org/rstan.html) worked better. The results closely match those from Epicure.

Using the data here: http://dwoll.de/err/dat.txt
The stan model fit below replicates the results from Epicure here: http://dwoll.de/err/epicure.log

Of course I am still interested in learning about other options or approaches.

Daniel

##-------
## rstan code for fitting an excess relative risk model with linear dose-response
## events = personYears * exp(beta0 + beta1*age + beta2*sex) * (1 +beta3*dose)

dat_code <- '
data {
    int<lower=0> N;
    vector[N] pyears;
    vector[N] age;
    vector[N] sex;
    vector[N] dose;
    int<lower=0> event[N];
}
transformed data {
    vector[N] log_pyears;
    log_pyears <- log(pyears);
}
parameters {
    vector[2] beta;
    real beta0;              # baseline intercept param
    real<lower=0> betaD;     # dose param -> non-negative
}
model {
    # beta0 unspecified -> uniform prior (improper)
    beta  ~ normal(0, 100);  # flat normal prior for params
    betaD ~ cauchy(0, 2.5);  # ok even if not truncated, cf. stan reference
    event ~ poisson_log(log_pyears + beta0 + beta[1]*age + beta[2]*sex +
                        log(1 + betaD*dose));
}
'

library(rstan)
stan_dat <- with(dat,
                 list(pyears=pyears,
                      age=age,
                      sex=sex,
                      N=length(pyears),
                      event=event,
                      dose=dose))

stanFit <- stan(model_code=dat_code, data=stan_dat,
                thin=5, iter=10000, chains=2)

traceplot(stanFit)
stanFit

-----Original Message-----
From: Wollschlaeger, Daniel 
Sent: Thursday, January 9, 2014 10:44 AM
To: David Winsemius
Cc: r-help at r-project.org
Subject: RE: AW: [R] Linear relative rate / excess relative risk models

Thanks for your suggestions! Here are links to simulated data and the Epicure syntax + reference fit:
http://dwoll.de/err/dat.txt
http://dwoll.de/err/epicure.log

The model tested in Epicure is

lambda = exp(alpha0 + alpha1*agePyr)*(1 + beta0*dosePyr*exp(beta1*agePyr))

with counts in variable event and offset pyears.

Many thanks, D

> -----Original Message-----
> From: David Winsemius [mailto:dwinsemius at comcast.net]
> Sent: Thursday, January 09, 2014 4:33 AM
> To: Wollschlaeger, Daniel
> Cc: r-help at r-project.org
> Subject: Re: AW: [R] Linear relative rate / excess relative risk models
> 
> 
> On Jan 8, 2014, at 3:22 PM, Wollschlaeger, Daniel wrote:
> 
> > If I understand you correctly, that is exactly the approach taken by
> > Atkinson & Therneau: They get the baseline rates from published rate
> > tables from the general population, multiply them by the appropriate
> > person-time from their data to get expected counts, and use this as
> > offset.
> >
> > Unfortunately, I won't have comparable baseline rate tables. And
> > while I could fit a separate model only to the unexposed group for
> > expected counts, I'd prefer to fit both factors (lambda0 and 1+ERR)
> > simultaneously - as it is typically done in the existing literature.
> 
> If you would describe your data situation more completely (ideally
> with a reproducible example) you might get a better answer. It's also
> considered polite on this mailing list to include the email chain, so
> appending original question:
> 
> --
> David
> 
> >
> > Best, Daniel
> >
> > ________________________________________
> > Von: David Winsemius [dwinsemius at comcast.net]
> > Gesendet: Mittwoch, 8. Januar 2014 19:06
> > An: Wollschlaeger, Daniel
> > Cc: r-help at r-project.org
> > Betreff: Re: [R] Linear relative rate / excess relative risk models
> >
> > I would fit a Poisson model to the dose-response data with offsets
> > for the baseline expecteds.
> 
> David Winsemius, MD
> Alameda, CA, USA
> 
> ============================
> My question is how I can fit linear relative rate models (= excess
> relative risk models, ERR) using R. In radiation epidemiology, ERR
> models are used to analyze dose-response relationships for event rate
> data and have the following form [1]:
> 
> lambda = lambda0(z, alpha) * (1 + ERR(x, beta))
> 
> * lambda is the event rate
> * lambda0 is the baseline rate function for non-exposed persons and
> depends on covariates z with parameters alpha
> * ERR is the excess relative risk function for exposed persons and
> depends on covariates x (among them dose) with parameters beta
> * lambda/lambda0 = 1 + ERR is the relative rate function
> 
> Often, the covariates z are a subset of the covariates x (like sex and
> age). lambda is assumed to be log-linear in lambda0, and ERR typically
> has a linear (or lin-quadratic) dose term as well as a log-linear
> modifying term with other covariates:
> 
> lambda0 = exp(alpha0 + alpha1*z1 + alpha2*z2 + ...)
> ERR = beta0*dose * exp(beta1*x1 + beta2*x2 + ...)
> 
> The data is often grouped in form of life tables with the observed
> event counts and person-years (pyr) for each cell that results from
> categorizing and cross-classifying the covariates. The counts are
> assumed to have a Poisson-distribution with mean mu = lambda*pyr, and
> the usual Poisson-likelihood is used. The interest is less in lambda0,
> but in inference on the dose coefficient beta0 and on the modifier
> coefficients beta.
> 
> In the literature, the specialized Epicure program is almost
> exclusively used. Last year, a similar question on R-sig-Epi [2] did
> not lead to a successful solution (I contacted the author). Atkinson &
> Therneau in [3] discuss excess risk models but get lambda0 separately
> from external data instead of fitting lambda0 as a log-linear term.
> Some R packages sound promising to me (eg., gnm, timereg) but I
> currently don't see how to correctly specify the model.
> 
> Any help on how to approach ERR models in R is highly appreciated!
> With many thanks and best regards