#' @title Formulas for local contributions to EJ metrics
#' @description Formulas for pop counts, pop demog pct, risk, cases, RR,
#' excess risk, excess cases (like EJ index sum),
#' and how those are different under 4 alt scenarios like E of D is set to that of dref.
#' @param P Population count per place (default is 1)
#' @param D Demographic group per place as fraction of P (default is 1)
#' @param E Environmental indicator (e.g., risk, exposure, or any local indicator, but if it is individual risk then number of cases can be calculated rather than just people-points)
#' @param na.rm Default is TRUE.
#' @param n Default is 5.
#' @param ... NOT USED. Would pass additional parameters to other functions
#' @export
ejRRadded <- function(E,D=1,P=1, na.rm=TRUE, n=5, ...) {
if (1==0) {
E=bg$cancer
D=bg$pctlowinc
P=bg$pop
na.rm=TRUE
}
###################################################################################################################### #
# **** This treats P as the full universe or denominator to use for D & d *****
# ******* this does not handle zero as denominator in zone ****
# **** need to handle NA in E here. .. **************
###################################################################################################################### #
######## #
# People (population counts)
######## #
# pop here
count <- P
# pop.overall
count.overall <- sumP <- sum(P, na.rm = na.rm)
# Dcount here
count.D <- PD <- P * D
PD[is.na(PD)] <- 0
count.D[is.na(count.D)] <- 0
# Dcount.overall
count.D.overall <- sumPD <- sum(PD)
# dcount here
count.dref <- Pd <- P - count.D # should be same as pct.dref * P, etc.
Pd[is.na(Pd)] <- 0
count.dref[is.na(count.dref)] <- 0
# dcount.overall
count.dref.overall <- sumPd <- sum(Pd)
######## #
# percentage of population
# (again, P must be proper universe)
######## #
# Dpct here
pct.D <- D
# Dpct.overall (******* this does not handle zero as denominator in zone ****)
pct.D.overall <- sumPD / sumP
# dpct here
pct.dref <- 1 - pct.D
# dpct.overall
pct.dref.overall <- sumPd / sumP
######## #
# cases # or "wtdE"
######## #
# cases here everybody **** need to handle NA in E here **************
cases <- PE <- P * E
PE[is.na(PE)] <- 0
cases[is.na(cases)] <- 0
# cases overall everybody
cases.overall <- sumPE <- sum(PE)
#(sum(cases)/1e6)/70 # if lifetime cancer risk per million is E
#cases.overall/1e6/70 # if lifetime cancer risk per million is E
# dref cases here
cases.dref <- PdE <- Pd * E
PdE[is.na(PdE)] <- 0
cases.dref[is.na(cases.dref)] <- 0
# dref cases overall
cases.dref.overall <- sumPdE <- sum(PdE)
# D cases here
cases.D <- PDE <- PD * E
cases.D[is.na(cases.D)] <- 0
PDE[is.na(PDE)] <- 0
cases.D[is.na(cases.D)] <- 0
# D cases overall
cases.D.overall <- sumPDE <- sum(PDE)
####### #
# E (individual risk)
######## #
risk <- E
risk.D <- E
risk.dref <- E
risk.overall <- sumPE / sumP # P-wtd.mean of E
risk.D.overall <- sumPDE / sumPD # PD-wtd.mean of E
risk.dref.overall <- sumPdE / sumPd # Pd-wtd.mean of E # one goal for EJ is for D overall to have this E
risk.elsewhere <- (sumPE - PE) / (sumP - P)
risk.D.elsewhere <- (sumPDE - PDE) / (sumPD - PD)
risk.dref.elsewhere <- (sumPdE - PdE) / (sumPd - Pd) # one goal for EJ is for D here to have this E
# riskodds <- risk / (1-risk) ???
# riskodds.overall <- ???
####### #
# Excess risk (has implicit goal)
######## #
risk.excess.overall <- risk.D.overall - risk.dref.overall
####### #
# Excess cases (has implicit goal)
######## #
cases.excess.overall <- risk.excess.overall * count.D.overall
####### #
# RR
######## #
#RR.overall <- (sumPDE / sumPD) / (sumPdE / sumPd)
RR.overall <- (risk.D.overall) / (risk.dref.overall)
# OR
#OR.overall <-
############################################### #
# ALTERNATIVE SCENARIOS, AND THEN CONTRIBUTION BECAUSE NOT AT THAT SCENARIO:
############################################### #
# if.here.was.zero means if E was zero and also P was zero !!! i.e. If nobody lived here at all! Not just if E = 0
egoal.v <- 0 # can be used to provide scenarios where local E is cut to zero
####### #
# cases in other scenarios
####### #
# THESE WORK EVEN IF egoal.v IS A VECTOR
cases.D.here.if.here.was.egoal <- count.D * egoal.v # = PD * egoal.v
#
cases.D.overall.if.here.was.egoal.vector <- cases.D.overall - cases.D + (egoal.v * count.D) # = sumPDE - PDE + (egoal.v * PD)
cases.dref.overall.if.here.was.egoal.vector <- cases.dref.overall - cases.dref + (egoal.v * count.dref)
cases.D.overall.if.here.was.risk.dref.overall <- cases.D.overall - cases.D + (risk.dref.overall * count.D)
cases.dref.overall.if.here.was.risk.dref.overall <- cases.dref.overall - cases.dref + (risk.dref.overall * count.dref)
cases.D.overall.if.here.was.risk.dref.elsewhere <- cases.D.overall - cases.D + (risk.dref.elsewhere * count.D)
cases.dref.overall.if.here.was.risk.dref.elsewhere <- cases.dref.overall - cases.dref + (risk.dref.elsewhere * count.dref)
cases.D.overall.if.here.was.zero <- cases.D.overall - cases.D + (0 * count.D) # = sumPDE - PDE
cases.dref.overall.if.here.was.zero <- cases.dref.overall - cases.dref + (0 * count.dref) # = sumPdE - PdE
####### #
# risk in other scenarios
####### #
risk.D.overall.if.here.was.egoal.vector <- (cases.D.overall.if.here.was.egoal.vector / (count.D.overall) )
risk.dref.overall.if.here.was.egoal.vector <- (cases.dref.overall.if.here.was.egoal.vector / (count.dref.overall) )
risk.D.overall.if.here.was.risk.dref.overall <- (cases.D.overall.if.here.was.risk.dref.overall / (count.D.overall) )
risk.dref.overall.if.here.was.risk.dref.overall <- (cases.dref.overall.if.here.was.risk.dref.overall / (count.dref.overall) )
risk.D.overall.if.here.was.risk.dref.elsewhere <- (cases.D.overall.if.here.was.risk.dref.elsewhere / (count.D.overall) )
risk.dref.overall.if.here.was.risk.dref.elsewhere <- (cases.dref.overall.if.here.was.risk.dref.elsewhere / (count.dref.overall) )
risk.D.overall.if.here.was.zero <- (cases.D.overall.if.here.was.zero / (count.D.overall - count.D) ) # i.e. if local pop was zero
risk.dref.overall.if.here.was.zero <- (cases.dref.overall.if.here.was.zero / (count.dref.overall - count.dref ) ) # i.e. if local pop was zero
####### #
# RR in other scenarios *** these use different versions of variables: but can use risk. variables above now that they exist ***
####### #
RR.if.here.was.egoal <- (
risk.D.overall.if.here.was.egoal.vector /
risk.dref.overall.if.here.was.egoal.vector
)
RR.if.here.was.risk.dref.overall <- (
risk.D.overall.if.here.was.risk.dref.overall /
risk.dref.overall.if.here.was.risk.dref.overall)
RR.if.here.was.risk.dref.elsewhere <- (
risk.D.overall.if.here.was.risk.dref.elsewhere /
risk.dref.overall.if.here.was.risk.dref.elsewhere
)
RR.if.here.was.zero <- (
risk.D.overall.if.here.was.zero /
risk.dref.overall.if.here.was.zero
)
RR.if.here.was.zero2 <- (
(sumPDE - PDE) / (sumPD - PD)) /
((sumPdE - PdE) / (sumPd - Pd))
####### #
# excess risk in other scenarios
####### #
#risk.excess.overall <- risk.D.overall - risk.dref.overall
risk.excess.if.here.was.egoal.vector <- risk.D.overall.if.here.was.egoal.vector - risk.dref.overall
risk.excess.if.here.was.risk.dref.overall <- risk.D.overall.if.here.was.risk.dref.overall - risk.dref.overall
risk.excess.if.here.was.risk.dref.elsewhere <- risk.D.overall.if.here.was.risk.dref.elsewhere - risk.dref.overall
risk.excess.if.here.was.zero <- risk.D.overall.if.here.was.zero - risk.dref.overall
####### #
# excess cases in other scenarios # **** need to check if these are correct and vectorize right: ****
####### #
cases.excess.if.here.was.egoal.vector <- risk.excess.if.here.was.egoal.vector * count.D.overall
cases.excess.if.here.was.risk.dref.overall <- risk.excess.if.here.was.risk.dref.overall * count.D.overall
cases.excess.if.here.was.risk.dref.elsewhere <- risk.excess.if.here.was.risk.dref.elsewhere * count.D.overall
cases.excess.if.here.was.zero <- risk.excess.if.here.was.zero * count.D.overall
############################################### #
# print some results for now
############################################### #
printstuff <- function(xname, x0,x1,x2,x3,x4){
cat('\n',
xname, '\n',
x0, ' now \n')
print(
head(
cbind(
if.egoal=x1,
if.here.was.risk.dref.overall=x2,
if.here.was.risk.dref.elsewhere=x3,
if.here.was.zero=x4
), n
)
)
print( summary(cbind(x1,x2,x3,x4)) )
}
xname='cases excess overall'
x0=cases.excess.overall
x1=cases.excess.if.here.was.egoal.vector
x2=cases.excess.if.here.was.risk.dref.overall
x3=cases.excess.if.here.was.risk.dref.elsewhere
x4=cases.excess.if.here.was.zero
printstuff(xname,x0,x1,x2,x3,x4)
xname='risk excess overall'
x0=risk.excess.overall
x1=risk.excess.if.here.was.egoal.vector ##### ok
x2=risk.excess.if.here.was.risk.dref.overall
x3=risk.excess.if.here.was.risk.dref.elsewhere
x4=risk.excess.if.here.was.zero
printstuff(xname,x0,x1,x2,x3,x4)
xname='RR now'
x0=RR.overall
x1=RR.if.here.was.egoal
x2=RR.if.here.was.risk.dref.overall
x3=RR.if.here.was.risk.dref.elsewhere
x4=RR.if.here.was.zero
printstuff(xname,x0,x1,x2,x3,x4)
xname='Risk D overall now'
x0=risk.D.overall
x1=risk.D.overall.if.here.was.egoal.vector
x2=risk.D.overall.if.here.was.risk.dref.overall
x3=risk.D.overall.if.here.was.risk.dref.elsewhere
x4=risk.D.overall.if.here.was.zero
printstuff(xname,x0,x1,x2,x3,x4)
xname='cases D overall now'
x0=cases.D.overall
x1=cases.D.overall.if.here.was.egoal.vector ## #
x2=cases.D.overall.if.here.was.risk.dref.overall ## #
x3=cases.D.overall.if.here.was.risk.dref.elsewhere ## #
x4=cases.D.overall.if.here.was.zero ## #
printstuff(xname,x0,x1,x2,x3,x4)
############## #
printstuff2 <- function() {
cat('\n')
cat('count.overall', 'count.dref.overall', 'count.D.overall: ', '\n')
cat(count.overall, count.dref.overall, count.D.overall, '\n','\n')
cat('pct.dref.overall', 'pct.D.overall:' , '\n')
cat(pct.dref.overall, pct.D.overall, '\n','\n')
cat('count.overall', 'count.dref.overall', 'count.D.overall:', '\n')
cat(count.overall, count.dref.overall, count.D.overall, '\n','\n')
cat('risk.overall', 'risk.D.overall', 'risk.dref.overall:', '\n')
cat(risk.overall, risk.D.overall, risk.dref.overall, '\n','\n')
}
printstuff2()
return()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.