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