This function organizes input and output for the analysis of attributable risk (for categorical variables). The analysis data, dframe, can be either a data frame or a simple features (sf) object. If an sf object is used, coordinates are extracted from the geometry column in the object, arguments xcoord and ycoord are assigned values "xcoord" and "ycoord", respectively, and the geometry column is dropped from the object.

attrisk_analysis(
  dframe,
  vars_response,
  vars_stressor,
  response_levels = NULL,
  stressor_levels = NULL,
  subpops = NULL,
  siteID = NULL,
  weight = "weight",
  xcoord = NULL,
  ycoord = NULL,
  stratumID = NULL,
  clusterID = NULL,
  weight1 = NULL,
  xcoord1 = NULL,
  ycoord1 = NULL,
  sizeweight = FALSE,
  sweight = NULL,
  sweight1 = NULL,
  fpc = NULL,
  popsize = NULL,
  vartype = "Local",
  conf = 95,
  All_Sites = FALSE
)

Arguments

dframe

Data to be analyzed (analysis data). A data frame or sf object containing survey design variables, response variables, stressor variables, and subpopulation (domain) variables.

vars_response

Vector composed of character values that identify the names of response variables in dframe. Each response variable must have two category values (levels), where one level is associated with poor condition and the other level is associated with good condition.

vars_stressor

Vector composed of character values that identify the names of stressor variables in dframe. Each stressor variable must have two category values (levels), where one level is associated with poor condition and the other level is associated with good condition.

response_levels

List providing the category values (levels) for each element in the vars_response argument. Each element in the list must contain two values, where the first value identifies poor condition, and the second value identifies good condition. This argument must be named and must be the same length as argument vars_response. Names for this argument must match the values in the vars_response argument. If this argument equals NULL, then a named list is created that contains the values "Poor" and "Good" for the first and second levels, respectively, of each element in the vars_response argument and that uses values in the vars_response argument as names for the list. If response_levels is provided without names, then the names of response_levels are set to vars_response. The default value is NULL.

stressor_levels

List providing the category values (levels) for each element in the vars_stressor argument. Each element in the list must contain two values, where the first value identifies poor condition, and the second value identifies good condition. This argument must be named and must be the same length as argument vars_stressor. Names for this argument must match the values in the vars_stressor argument. If this argument equals NULL, then a named list is created that contains the values "Poor" and "Good" for the first and second levels, respectively, of each element in the vars_stressor argument and that uses values in the vars_stressor argument as names for the list. If stressor_levels is provided without names, then the names of stressor_levels are set to vars_stressor. The default value is NULL.

subpops

Vector composed of character values that identify the names of subpopulation (domain) variables in dframe. If a value is not provided, the value "All_Sites" is assigned to the subpops argument and a factor variable named "All_Sites" that takes the value "All Sites" is added to dframe. The default value is NULL.

siteID

Character value providing the name of the site ID variable in dframe. For a two-stage sample, the site ID variable identifies stage two site IDs. The default value is NULL, which assumes that each row in dframe represents a unique site.

weight

Character value providing the name of the design weight variable in dframe. For a two-stage sample, the weight variable identifies stage two weights. The default value is "weight".

xcoord

Character value providing name of the x-coordinate variable in dframe. For a two-stage sample, the x-coordinate variable identifies stage two x-coordinates. Note that x-coordinates are required for calculation of the local mean variance estimator. If dframe is an sf object, this argument is not required (as the geometry column in dframe is used to find the x-coordinate). The default value is NULL.

ycoord

Character value providing name of the y-coordinate variable in dframe. For a two-stage sample, the y-coordinate variable identifies stage two y-coordinates. Note that y-coordinates are required for calculation of the local mean variance estimator. If dframe is an sf object, this argument is not required (as the geometry column in dframe is used to find the t-coordinate). The default value is NULL.

stratumID

Character value providing the name of the stratum ID variable in dframe. The default value is NULL.

clusterID

Character value providing the name of the cluster (stage one) ID variable in dframe. Note that cluster IDs are required for a two-stage sample. The default value is NULL.

weight1

Character value providing the name of the stage one weight variable in dframe. The default value is NULL.

xcoord1

Character value providing the name of the stage one x-coordinate variable in dframe. Note that x coordinates are required for calculation of the local mean variance estimator. The default value is NULL.

ycoord1

Character value providing the name of the stage one y-coordinate variable in dframe. Note that y-coordinates are required for calculation of the local mean variance estimator. The default value is NULL.

sizeweight

Logical value that indicates whether size weights should be used during estimation, where TRUE uses size weights and FALSE does not use size weights. To employ size weights for a single-stage sample, a value must be supplied for argument weight. To employ size weights for a two-stage sample, values must be supplied for arguments weight and weight1. The default value is FALSE.

sweight

Character value providing the name of the size weight variable in dframe. For a two-stage sample, the size weight variable identifies stage two size weights. The default value is NULL.

sweight1

Character value providing the name of the stage one size weight variable in dframe. The default value is NULL.

fpc

Object that specifies values required for calculation of the finite population correction factor used during variance estimation. The object must match the survey design in terms of stratification and whether the design is single-stage or two-stage. For an unstratified design, the object is a vector. The vector is composed of a single numeric value for a single-stage design. For a two-stage unstratified design, the object is a named vector containing one more than the number of clusters in the sample, where the first item in the vector specifies the number of clusters in the population and each subsequent item specifies the number of stage two units for the cluster. The name for the first item in the vector is arbitrary. Subsequent names in the vector identify clusters and must match the cluster IDs. For a stratified design, the object is a named list of vectors, where names must match the strata IDs. For each stratum, the format of the vector is identical to the format described for unstratified single-stage and two-stage designs. Note that the finite population correction factor is not used with the local mean variance estimator.

Example fpc for a single-stage unstratified survey design:

fpc <- 15000

Example fpc for a single-stage stratified survey design:

fpc <- list( Stratum_1 = 9000, Stratum_2 = 6000)

Example fpc for a two-stage unstratified survey design:

fpc <- c( Ncluster = 150, Cluster_1 = 150, Cluster_2 = 75, Cluster_3 = 75, Cluster_4 = 125, Cluster_5 = 75)

Example fpc for a two-stage stratified survey design:

fpc <- list( Stratum_1 = c( Ncluster_1 = 100, Cluster_1 = 125, Cluster_2 = 100, Cluster_3 = 100, Cluster_4 = 125, Cluster_5 = 50), Stratum_2 = c( Ncluster_2 = 50, Cluster_1 = 75, Cluster_2 = 150, Cluster_3 = 75, Cluster_4 = 75, Cluster_5 = 125))

popsize

Object that provides values for the population argument of the calibrate or postStratify functions in the survey package. If a value is provided for popsize, then either the calibrate or postStratify function is used to modify the survey design object that is required by functions in the survey package. Whether to use the calibrate or postStratify function is dictated by the format of popsize, which is discussed below. Post-stratification adjusts the sampling and replicate weights so that the joint distribution of a set of post-stratifying variables matches the known population joint distribution. Calibration, generalized raking, or GREG estimators generalize post-stratification and raking by calibrating a sample to the marginal totals of variables in a linear regression model. For the calibrate function, the object is a named list, where the names identify factor variables in dframe. Each element of the list is a named vector containing the population total for each level of the associated factor variable. For the postStratify function, the object is either a data frame, table, or xtabs object that provides the population total for all combinations of selected factor variables in the dframe data frame. If a data frame is used for popsize, the variable containing population totals must be the last variable in the data frame. If a table is used for popsize, the table must have named dimnames where the names identify factor variables in the dframe data frame. If the popsize argument is equal to NULL, then neither calibration nor post-stratification is performed. The default value is NULL.

Example popsize for calibration:

popsize <- list( Ecoregion = c( East = 750, Central = 500, West = 250), Type = c( Streams = 1150, Rivers = 350))

Example popsize for post-stratification using a data frame:

popsize <- data.frame( Ecoregion = rep(c("East", "Central", "West"), rep(2, 3)), Type = rep(c("Streams", "Rivers"), 3), Total = c(575, 175, 400, 100, 175, 75))

Example popsize for post-stratification using a table:

popsize <- with(MySurveyFrame, table(Ecoregion, Type))

Example popsize for post-stratification using an xtabs object:

popsize <- xtabs(~Ecoregion + Type, data = MySurveyFrame)

vartype

Character value providing the choice of the variance estimator, where "Local" indicates the local mean estimator and "SRS" indicates the simple random sampling estimator. The default value is "Local".

conf

Numeric value providing the Gaussian-based confidence level. The default value is 95.

All_Sites

A logical variable used when subpops is not NULL. If All_Sites is TRUE, then alongside the subpopulation output, output for all sites (ignoring subpopulations) is returned for each variable in vars. If All_Sites is FALSE, then alongside the subpopulation output, output for all sites (ignoring subpopulations) is not returned for each variable in vars. The default is FALSE.

Value

The analysis results. A data frame of population estimates for all combinations of subpopulations, categories within each subpopulation, response variables, and categories within each response variable. Estimates are provided for proportion and size of the population plus standard error, margin of error, and confidence interval estimates. The data frame contains the following variables:

Type

subpopulation (domain) name

Subpopulation

subpopulation name within a domain

Response

response variable

Stressor

stressor variable

nResp

sample size

Estimate

attributable risk estimate

StdError_log

attributable risk standard error (on the log scale)

MarginofError_log

attributable risk margin of error (on the log scale)

LCBxxPct

xx% (default 95%) lower confidence bound

UCBxxPct

xx% (default 95%) upper confidence bound

WeightTotal

sum of design weights

Count_RespPoor_StressPoor

number of observations in the poor response and poor stressor group

Count_RespPoor_StressGood

number of observations in the poor response and good stressor group

Count_RespGood_StressPoor

number of observations in the good response and poor stressor group

Count_RespGood_StressGood

number of observations in the good response and good stressor group

Prop_RespPoor_StressPoor

weighted proportion of observations in the poor response and poor stressor group

Prop_RespPoor_StressGood

weighted proportion of observations in the poor response and good stressor group

Prop_RespGood_StressPoor

weighted proportion of observations in the good response and poor stressor group

Prop_RespGood_StressGood

weighted proportion of observations in the good response and good stressor group

Details

Attributable risk measures the proportional reduction in the extent of poor condition of a response variable that presumably would result from eliminating a stressor variable, where the response and stressor variables are classified as either good (i.e., reference condition) or poor (i.e., different from reference condition). Attributable risk is defined as one minus the ratio of two probabilities. The numerator of the ratio is the conditional probability that the response variable is in poor condition given that the stressor variable is in good condition. The denominator of the ratio is the probability that the response variable is in poor condition. Attributable risk values close to zero indicate that removing the stressor variable will have little or no impact on the probability that the response variable is in poor condition. Attributable risk values close to one indicate that removing the stressor variable will result in extensive reduction of the probability that the response variable is in poor condition.

References

Sickle, J. V., & Paulsen, S. G. (2008). Assessing the attributable risks, relative risks, and regional extents of aquatic stressors. Journal of the North American Benthological Society, 27(4), 920-931.

See also

relrisk_analysis

for relative risk analysis

diffrisk_analysis

for risk difference analysis

Author

Tom Kincaid Kincaid.Tom@epa.gov

Examples

dframe <- data.frame(
  siteID = paste0("Site", 1:100),
  wgt = runif(100, 10, 100),
  xcoord = runif(100),
  ycoord = runif(100),
  stratum = rep(c("Stratum1", "Stratum2"), 50),
  RespVar1 = sample(c("Poor", "Good"), 100, replace = TRUE),
  RespVar2 = sample(c("Poor", "Good"), 100, replace = TRUE),
  StressVar = sample(c("Poor", "Good"), 100, replace = TRUE),
  All_Sites = rep("All Sites", 100),
  Resource_Class = rep(c("Agr", "Forest"), c(55, 45))
)
myresponse <- c("RespVar1", "RespVar2")
mystressor <- c("StressVar")
mysubpops <- c("All_Sites", "Resource_Class")
attrisk_analysis(dframe,
  vars_response = myresponse,
  vars_stressor = mystressor, subpops = mysubpops, siteID = "siteID",
  weight = "wgt", xcoord = "xcoord", ycoord = "ycoord",
  stratumID = "stratum"
)
#>             Type Subpopulation Response  Stressor nResp     Estimate
#> 1      All_Sites     All Sites RespVar1 StressVar   100 -0.146822490
#> 2      All_Sites     All Sites RespVar2 StressVar   100 -0.029705172
#> 3 Resource_Class           Agr RespVar1 StressVar    55 -0.220942488
#> 4 Resource_Class        Forest RespVar1 StressVar    45 -0.072784416
#> 5 Resource_Class           Agr RespVar2 StressVar    55 -0.044451332
#> 6 Resource_Class        Forest RespVar2 StressVar    45 -0.006427495
#>   StdError_log MarginofError_log  LCB95Pct  UCB95Pct WeightTotal
#> 1    0.3746480         0.7342965 -1.389996 0.4497054    5571.035
#> 2    0.3825463         0.7497770 -1.179400 0.5134933    5571.035
#> 3    0.5034165         0.9866781 -2.274945 0.5448167    2944.924
#> 4    0.5775314         1.1319407 -2.327423 0.6541268    2626.110
#> 5    0.5252535         1.0294780 -1.924050 0.6269289    2944.924
#> 6    0.5317098         1.0421320 -1.853479 0.6450311    2626.110
#>   Count_RespPoor_StressPoor Count_RespPoor_StressGood Count_RespGood_StressPoor
#> 1                        23                        34                        25
#> 2                        24                        25                        24
#> 3                        11                        20                        14
#> 4                        12                        14                        11
#> 5                        14                        18                        11
#> 6                        10                         7                        13
#>   Count_RespGood_StressGood Prop_RespPoor_StressPoor Prop_RespPoor_StressGood
#> 1                        18                0.2184026                0.3360798
#> 2                        27                0.2182399                0.2605832
#> 3                        10                0.1873002                0.3371157
#> 4                         8                0.2532809                0.3349180
#> 5                        12                0.2600568                0.3177389
#> 6                        15                0.1713464                0.1964886
#>   Prop_RespGood_StressPoor Prop_RespGood_StressGood
#> 1                0.2530810                0.1924366
#> 2                0.2532438                0.2679332
#> 3                0.2861882                0.1893958
#> 4                0.2159546                0.1958464
#> 5                0.2134317                0.2087726
#> 6                0.2978892                0.3342759