risk(f, g, delta = 0, log = TRUE, h = NULL, adaptive = FALSE, res = 50,
WIN = NULL, tolerate = FALSE, plotit = TRUE, comment = TRUE)"bivden" representing the `case' density estimate, or an object of type data.frame, list, f, but for the controls. Whatever the type, the class of g must match that of f.delta requests no additive constant (default).TRUE with the alternative being the raw density ratio.f and g are already "bivden" objects. An optional numeric vector of length 1 OR 2, giving the global bandwidth(s) for internal estimation of the case and control densities if adaptive = TRUE, f and g are already "bivden" objects. A boolean value specifying whether or not to employ adaptive smoothing for internally estimating the densities. A value of FALSE (default) elects use of ff and g are already "bivden" objects. A numeric value giving the desired resolution (of one side) of the evaluation grid. Higher values increase resolution at the expense of computational efficiency. Defauf and g are already "bivden" objects OR objects of class ppp (in which case the study region is set to the value of the resident window comf and g are already "bivden" objects. A boolean value specifying whether or not to calculate a corresponding asymptotic p-value surface (for tolerance contours) for the estimated relative risk function. IfTRUE (default), a heatplot of the estimated relative risk function is produced. If tolerate = TRUE, asymptotic tolerance contours are automatically added to the plot at a significance level of 5%.f and g are already "bivden" objects. Boolean. Whether or not to print function progress (including starting and ending date-times) during execution. Defaults to TRUE."rrs". This is a marked list with the following components:res*res matrix (where res is the grid resolution as specified in the calls to bivariate.density for calculation of f and g) giving the values of the risk surface over the evaluation grid. Values corresponding to grid coordinates outside the study region are assigned NA"bivden" used as the numerator or `case' density estimate"bivden" used as the denominator or `control' density estimate"bivden" (based on the pooled data) calculated internally if f and g were raw data arguments, NA otherwisetolerate = TRUE and f and g were raw data arguments, NA otherwiserisk, as opposed to previously computed objects of class "bivden", the running time of this function will be greater. This is particularly the case if the user has also selected tolerate = TRUE. In the same fashion as bivariate.density and tolerance, setting comment = TRUE can keep the user appraised of the function progress during run-time.bivariate.density. In geographical epidemiology, the argument f represents the spatial distribution of the disease cases, and g the at-risk (control) population.
The option to supply the raw case and control data is available. If this is done, the function runs bivariate.density internally, abstracting certain decisions about the density estimation away from the user. If the user sets adaptive = TRUE (and h remains at NULL), the smoothing parameters are calculated as per the approach taken in Davies and Hazelton (2010): a common global bandwidth using the pooled data from OS, and differing pilot bandwidths using CV.sm on the case and control data separately. The scaling parameter gamma is common for the case and control density estimates, set as the gamma component of the pooled estimate. If a fixed relative risk is desired (adaptive = FALSE) and no specific bandwidths are given via the argument h, the case and control densities share a common bandwidth computed from the pooled data using OS. In supplying raw data to risk, the user must also specify an evaluation grid resolution (defaulting to 50 by 50) and the study region WIN (unless f and g are objects of class ppp, in which case the resident window component overrides WIN). All other arguments are set to their defaults as in bivariate.density.
If more flexibility is required for estimation of the case and control densities, the user must supply `pre-calculated' objects of class "bivden" (from bivariate.density) as the f and g arguments. This drastically reduces the running time of a call to risk (as the density estimation step is already complete). However, the option of internally computing the asymptotic p-value surfaces (via the argument tolerate) is unavailable in this case; the user must run the tolerance function separately if tolerance contours are desired.
The relative risk function is defined here as the ratio of the `case' density to the `control' (Bithell, 1990; 1991). Using kernel density estimation to model these densities (Diggle, 1985), we obtain a workable estimate thereof. This function defines the risk function r in the following fashion:
r= (f + delta*max(g))/(g + delta*max(g))
Note the (optional) additive constants defined by delta times the maximum of each of the densities in the numerator and denominator respectively (see Bowman and Azzalini, 1997).
The log-risk function rho, given by rho = log[r], is argued to be preferable in practice as it imparts a sense of symmetry in the way the case and control densities are treated (Kelsall and Diggle, 1995a;b). The option of log-transforming the returned risk function is therefore selected by default.data(PBC)
PBC.casedata <- PBC$data[PBC$data[,3]==1,1:2]
PBC.controldata <- PBC$data[PBC$data[,3]==0,1:2]
pbc.h <- OS(PBC$data[,1:2],
nstar = sqrt(nrow(PBC.casedata) * nrow(PBC.controldata)))
pbc.pool <- bivariate.density(data = PBC$data[,1:2], pilotH = pbc.h,
adaptive = FALSE, WIN = PBC$owin)
pbc.case <- bivariate.density(data = PBC.casedata, pilotH = pbc.h,
adaptive = FALSE, WIN = PBC$owin)
pbc.con <- bivariate.density(data = PBC.controldata, pilotH = pbc.h,
adaptive = FALSE, WIN = PBC$owin)
pbc.rrs <- risk(f = pbc.case, g = pbc.con, plotit = FALSE)
summary(pbc.rrs)Run the code above in your browser using DataLab