A transformation can be used to reparametrize parameters internally
(i.e., such that the likelihoods are evaluated at
transform(parameters), where parameters are
unconstrained), and should be specified when fitting linear excess
relative risk and linear-exponential models to ensure nonnegative
odds/risk/hazard. The included function transform1 applies
an exponential transformation to the desired parameters, see below. When
supplying a function to transform, this should be a
function of the full parameter vector, returning a full (transformed)
parameter vector. In particular, the full parameter vector contains
parameters in the following order (see ?ameras for the
model definitions): \(\alpha_0, \mathbf
\alpha, \beta_1, \beta_2, \mathbf \beta_{m1}, \mathbf \beta_{m2},
\sigma\), where \(\mathbf
\alpha\), \(\mathbf\beta_{m1}\)
and \(\mathbf \beta_{m2}\) can be
vectors, with lengths matching X and M,
respectively. \(\sigma\) is only
included for the linear model (Gaussian family), and no intercept is
included for the Cox and conditional logistic models. For the
multinomial model, the full parameter vector is the concatenation of
\(Z-1\) parameter vectors in the order
as given above, where \(Z\) is the
number of outcome categories. When no transformation is specified and
the linear ERR model is used, transform1 is used for ERR
parameters \(\beta_1\) and \(\beta_2\) by default, with lower limits
\(-1/max(D)\) for linear dose-response
and \((0,-1/max(D^2))\) for
linear-quadratic dose-response, respectively (see below). For the
linear-exponential model, a lower limit of 0 is used for \(\beta_1\), and no transformation is used
for \(\beta_2\). If effect modifiers
M are specified, no transformation is used for those
parameters. When negative RRs are obtained during optimization, an error
will be generated and a different transformation or bounds should be
used. All output is returned in the original parametrization given in
?ameras. The Jacobian of the transformation
(transform.jacobian) is required when using a
transformation with methods other than BMA. For transform1,
the Jacobian is given by transform1.jacobian.
The included function transform1 applies the exponential
transformation \(f(\theta_i)=\exp(\theta_i)+LB_i\) to one or
multiple components of parameter vector \(\mathbf \theta\), where \(LB_i\) are lower limits that can be
different for each component. In particular, a vector of indices of
parameters to be transformed and a vector of corresponding lower bounds
LB can be supplied to arguments index.t and
lowlimit, respectively, resulting in transformed parameters
\(f(\theta_i)=\exp(\theta_i)+\text{LB}_i\).
In particular, transform1 and
transform1.jacobian are defined as follows:
transform1 <- function(params, index.t=1:length(params), lowlimit=rep(0,length(index.t)),
boundcheck=FALSE, boundtol=1e-3, ...){
if(length(index.t)!=length(lowlimit))
stop("Length mismatch between index.t and lowlimit")
if(any(!(index.t %in% 1:length(params))))
stop("Incorrect indices for transformation specified")
params[index.t] <- exp(params[index.t]) + lowlimit
if(boundcheck){
if(any(params[index.t]-lowlimit < boundtol))
warning(paste0("WARNING: one or multiple parameter estimates within ", boundtol, " of
lower bounds. Try different bounds or starting values."))
}
return(params)
}
transform1.jacobian <- function(params, index.t=1:length(params), ...){
if(any(!(index.t %in% 1:length(params))))
stop("Incorrect indices for transformation specified")
grad <- rep(1, length(params))
grad[index.t] <- exp(params[index.t])
if(length(params)>1){
return(diag(grad))
} else{
return(matrix(grad))
}
}If you wish to supply your own transformation, it is helpful to start
from the definition of transform1. It is also important to
keep the following in mind:
params
and ..., where params is the full parameter
vector. Parameters are in a specific order (see above), in case of doubt
it is always possible to run ameras with the default
settings to verify correct the order from the result.amerasboundcheck is an argument to the function,
it should be a logical. When transforming parameters after the
optimization, the transformation is called with
boundcheck=TRUE. If boundcheck is not an
argument of the function, this is ignored.See the definition of transform1 above for an example of
how to apply the transformation only to specific parameters, and how to
use extra arguments.
As an example, suppose instead of the exponential transformation from
transform1, for the parameters \(\beta_1\) and \(\beta_2\) we wish to use the sigmoid
transformation \(f: \mathbb{R} \rightarrow
(a,b)\) given by \[f_{a_i,b_i}(\theta_i)= a_i + (b_i-a_i)
\frac{1}{1+\exp(-\theta_i)}.\] Then, using
transform1 as a starting point, we can define the
transformation as follows (note that since \(\mathbf{a}\) and \(\mathbf{b}\) act as bounds, we use an
updated bound check):
transform.sigmoid <- function(params, index.t=1:length(params), a=rep(0,length(index.t)),
b=rep(1,length(index.t)), boundcheck=FALSE, boundtol=1e-3, ...){
if(length(index.t)!=length(a) | length(index.t) != length(b))
stop("Length mismatch between index.t, a, and b")
if(any(!(index.t %in% 1:length(params))))
stop("Incorrect indices for transformation specified")
params[index.t] <- a + (b-a) * 1/(1+exp(-1*params[index.t]))
if(boundcheck){
if(any( (params[index.t]-a < boundtol) | (b-params[index.t] < boundtol)))
warning(paste0("WARNING: one or multiple parameter estimates within ", boundtol,
" of bounds. Try different bounds or starting values."))
}
return(params)
}Next, noting that \(\mathrm{d}f_{a_i,b_i}/d\theta_i =
(b_i-a_i)\exp(-\theta_i)/\{1+\exp(-\theta_i) \}^2\), we can
define the Jacobian as follows, using transform1.jacobian
as a starting point:
transform.sigmoid.jacobian <- function(params, index.t=1:length(params),
a=rep(0,length(index.t)), b=rep(1,length(index.t)), ...){
if(length(index.t)!=length(a) | length(index.t) != length(b))
stop("Length mismatch between index.t, a, and b")
if(any(!(index.t %in% 1:length(params))))
stop("Incorrect indices for transformation specified")
grad <- rep(1, length(params))
grad[index.t] <- (b-a)*exp(-1*params[index.t])/(1+exp(-1*params[index.t]))^2
if(length(params)>1){
return(diag(grad))
} else{
return(matrix(grad))
}
}Now let us try this transformation on the example data using
regression calibration for a linear-quadratic ERR model. Note that all
parameters are returned after transformation, and so there should be no
difference between using transform1 and
transform.sigmoid. First, we fit the model using the
sigmoid transformation:
data(data, package="ameras")
dosevars <- paste0("V", 1:10)
fit.ameras.sigmoid <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data,
family="binomial", deg=2, doseRRmod = "ERR", methods="RC",
transform=transform.sigmoid, transform.jacobian=transform.sigmoid.jacobian,
index.t=4:5)
#> Fitting RC
#> Obtaining profile likelihood CI for dose
#> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
#> WARNING: Lower bound for dose is < 0 and may not exist if rescaling the
#> variable does not help
#> Obtaining profile likelihood CI for dose_squared
summary(fit.ameras.sigmoid)
#> Call:
#> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars,
#> X = c("X1", "X2"), methods = "RC", deg = 2, doseRRmod = "ERR",
#> transform = transform.sigmoid, transform.jacobian = transform.sigmoid.jacobian,
#> index.t = 4:5)
#>
#> Total run time: 3 seconds
#>
#> Runtime in seconds by method:
#>
#> Method Runtime
#> RC 3.0
#>
#> Summary of coefficients by method:
#>
#> Method Term Estimate SE CI.lowerbound CI.upperbound
#> RC (Intercept) -0.87358 0.09758 NA NA
#> RC X1 0.44587 0.07672 NA NA
#> RC X2 -0.33552 0.09610 NA NA
#> RC dose 0.04875 0.21281 0.0000 0.5130
#> RC dose_squared 0.28764 0.08100 0.1328 0.4121Next with default settings, using transform1:
fit.ameras.transform1 <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data,
family="binomial", deg=2, doseRRmod = "ERR", methods="RC")
#> Fitting RC
#> Obtaining profile likelihood CI for dose
#> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
#> WARNING: Lower bound for dose is < 0 and may not exist if rescaling the
#> variable does not help
#> Obtaining profile likelihood CI for dose_squared
summary(fit.ameras.transform1)
#> Call:
#> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars,
#> X = c("X1", "X2"), methods = "RC", deg = 2, doseRRmod = "ERR")
#>
#> Total run time: 3.3 seconds
#>
#> Runtime in seconds by method:
#>
#> Method Runtime
#> RC 3.3
#>
#> Summary of coefficients by method:
#>
#> Method Term Estimate SE CI.lowerbound CI.upperbound
#> RC (Intercept) -0.87359 0.09759 NA NA
#> RC X1 0.44587 0.07672 NA NA
#> RC X2 -0.33552 0.09610 NA NA
#> RC dose 0.04878 0.21283 0.0000 0.5115
#> RC dose_squared 0.28763 0.08100 0.1325 0.4108As expected, there is no difference between the results.