Bounding analysis applied to attributed risk

Minh Ha-Duong, 2002, 2003

This notebook demonstrates a bounding analysis method to attribute the lung cancer risk to various environmental pollutants.

This code is freely re-useable according to the GNU General Public License.
This notebook should be distributed along with the  BeliefMatrices.m  package because it will not run without it.

The corresponding theory is explained in:
Granger-Casman NSF research grant number SES-0216897 and SRB-9521914
Ha-Duong, Casman and Granger, (2003) Bounding Analysis applied to lung cancer risk (submitted to ISIPTA´03, Lugano)
Ha-Duong, Casman and Granger, (2003) Bounding analysis as a companion for Risk Analysis: Theory and methods (for submission to Risk Analysis)

The Frame of reference

Ω = {smoke, radon, asbestos, 
          unknown};
n = Length[Ω];

Subsets of Ω are indexed in the canonical order (canonical only in this theory and for Mathematica only up to version 4.1 included):
using binary notation from 0 to (2^n) - 1, with 0 corresponding to the empty set and  (2^n) - 1 to Ω.
We load the standard package to index subsets this way, with  RankSubset and NthSubset   functions. 
It also defines the  KSubsets  function that returns all subsets with  k  elements.

Needs["DiscreteMath`Combinatorica`"]

Attributable risk: Upper bound is plausibility, lower bound is belief

Next we load the package implemening belief function theory for small n. Small means n less than 10, because the package implements Philippe Smets (2001) Matrix calculus for belief functions, which uses square matrices of size  (2^n) - 1.

Needs["BeliefMatrices`"]

Here is the central idea, see Ha-Duong and ... (2002) Discussion paper on bounding analysis:

The central variable is m , representing the distribution of exposure among all recorded cases.
For example,  m[smoke] = 0.50  means that within all lung cancer deaths in the reference year,
fifty percent of cases were exposed to  smoke  and no other factor.
We adopt m  as the lower bound on attributable risk, and 
plausibility for the bound that means the UPPER bound.
They are corrected by the background rate, which is NOT m[0], see below.

bound[attributableRisk[cause_Symbol], 
        backgroundRate_] := Plausibility[m, n][[
          RankBinarySubset[Ω, {cause}]]] (1 - backgroundRate)

    
bound[attributableRisk[cause1_  &&  cause2_], 
        backgroundRate_] := Commonality[m, n][[
          RankBinarySubset[Ω, {cause1, 
              cause2}]]] (1 - 
            backgroundRate)

    
bound[attributableRisk[cause1_  ||  cause2_], 
        backgroundRate_] := Plausibility[m, n][[
          RankBinarySubset[Ω, {cause1, 
              cause2}]]] (1 - 
            backgroundRate)

    
lowerBound[attributableRisk[cause_Symbol], backgroundRate_] := 
      m[RankBinarySubset[Ω, {cause}]] (1 - 
            backgroundRate)

    
lowerBound[attributableRisk[cause1_  &&  cause2_], 
        backgroundRate_] := 
      m[RankBinarySubset[Ω, {cause1, cause2}]] (1 - 
            backgroundRate)

    
lowerBound[attributableRisk[cause1_  ||  cause2_], 
        backgroundRate_] := Belief[m, n][[
          RankBinarySubset[Ω, {cause1, 
              cause2}]]] (1 - 
            backgroundRate)

Coherence constraints

See Walley (1991) page 198.
Note on the data structure: information about lower and upper probabilities shall be given as a list of pairs:
  boundsList =  {{lowerProba1, upperProba1 }, {lowerProba2, 
          upperProba2 },  ... }

sensibleBoundQ[{lowerBound_, upperBound_}] := 
      0 ≤ lowerBound ≤ upperBound ≤ 
        1

avoidSureLossQ[
        boundsList_] := (Plus @@ 
              boundsList)[[1]] ≤ 
        1 ≤ (Plus @@ 
              boundsList)[[2]]

coherentQ[boundsList_, 
        i_] := (Plus @@ 
                boundsList)[[1]] - 
          boundsList[[i, 1]] + 
          boundsList[[i, 
            2]] ≤ 
        1 ≤  (Plus @@ 
                boundsList)[[2]] - 
          boundsList[[i, 2]] + 
          boundsList[[i, 1]]
coherentQ[boundsList_] := 
      Table[coherentQ[boundsList, i], {i, 1, Length[boundsList]}]

coherenceConstraints[boundsList_] := 
      Flatten[{sensibleBoundQ /@ 
            boundsList, 
avoidSureLossQ[
            boundsList], 
coherentQ[
            boundsList]}]

Maximum unspecificity principle

Constraint that no death is related to  num  factors or more.
This is a simplification. In a further study we will introduce the concept of interactions/independance/irrelevance to replace this constraint.

noMultipleSynergy[
        num_] := (m[RankBinarySubset[Ω, #]] == 
              0) &  /@ (Join  @@  
            Table[KSubsets[Ω, k], {k, num, n}])

In addition we specify that m is positive and sums up to unity.

constraintSet[knowledge__] := 
      Join[Flatten[{knowledge}], 
Table[
          m[k] ≥ 0,  {k, 0, 
            (2^n) - 1}], 
{Sum[
              m[k],  {k, 0, (2^n) - 1}]  ==  
            1 }
]

Lets simplify the constraints.

simpleConstraints[set_] := 
      Cases[set, m[_] == _?NumberQ]
presolve[
        set_] := 
Join[simpleConstraints[set], 
          set /. (simpleConstraints[
                  set] /. {Equal ->  Rule })] //. {a___, True, 
            b___} -> {a, b}

We retain the most unspecific m compatible with available knowledge (given as a set of constraints).
Unicity is NOT garanteed.
There is a solution to the maximization program if and only if the coherence constraints can be satisfied.

mostUnspecificExposure[
          knowledge__] := ConstrainedMax[ 
            Unspecificity[m, n], 
presolve[
              constraintSet[knowledge]], 
Table[
              m[k],  {k, 0, 
                (2^n) - 1}]
][[2]];

Relative risk

The relative risk of smoking says how many times more likely you are to get cancer when you smoke compared to non-smokers.
The formula is   ar ==  ((p (rr - 1))/(p (rr - 1) + 1))  
where ar is the attributable risk, rr the relative risk and p the probability of exposure,
see for example Beir VI appendix C p. 229, or Breslow & Day Measures of Disease p. 74

attributableFromRelativeRisk[cause_, rr_] = 
        1 - (1/(exposureProbability[cause] (rr - 1)  + 1));

relativeFromAttributableRisk[cause_, ar_] = 
      rr /. Solve[ar == attributableFromRelativeRisk[cause, rr], 
            rr][[1]]

It is the same to bound ar or rr, since they both increase together.

Here is how to say that the relative risk of smoke is between 3 and 5,
so that the inequality stays linear in our control variables m.

bound[attributableRisk[smoke]] ≤ 
      attributableFromRelativeRisk[smoke, 5]

lowerBound[attributableRisk[smoke], br] >= 
      attributableFromRelativeRisk[smoke, 3]

Pretty printing

toOneDigit[x_] := Round[10  x]/10. 
toTwoDigits[x_] := Round[100  x]/100. 
toThreeDigits[x_] := Round[1000  x]/1000. 
Attributes[toOneDigit] = Listable;
Attributes[toTwoDigits] = Listable;
Attributes[toThreeDigits] = Listable;

dataToDisplay := {
toThreeDigits[bound[
                  attributableRisk[#], 
                  backgroundRate] & /@ Ω], 
toThreeDigits[lowerBound[attributableRisk[#], 
                  backgroundRate] & /@ Ω], 
exposureProbability[#] &  /@ Ω, 
toTwoDigits[
                relativeFromAttributableRisk[#, 
                  bound[attributableRisk[#], 
                    backgroundRate]]] & /@ Ω, 
toTwoDigits[
                relativeFromAttributableRisk[#, 
                  lowerBound[attributableRisk[#], 
                    backgroundRate]]] & /@ Ω
 (*,PignisticProbability[m, 
              n]*) 
};

\(resultSummary[
        rule_] := 
ColumnForm[{TableForm[
dataToDisplay /. 
              rule, 
TableHeadings -> {{\*""\!\(attributable\&_\) %"", \*""\!\(attributable\+_\) %"", "Exposure %", \*""\!\(\(relative\\ risk\)\&_\)"", \*""\!\(\(relative\\ risk\)\+_\)""\  (*\(,\)\("betP"\)*) }, Ω}
], 
{"Non attributed", \ Chop[m[0] /. rule]}
}]\)

resultDetail[
        rule_] := 
rule /. Rule[m[k_], 
            rr_]->{k, UnrankBinarySubset[k, Ω], 
            rr}

statTable[
        rule_] := {Rest[BinarySubsets[Ω]], 
          Table[m[k] (1 - backgroundRate), {k, 1, (2^n) - 1}], 
          Belief[m, n] (1 - backgroundRate), 
          Commonality[m, n] (1 - backgroundRate), 
          Plausibility[m, 
              n] (1 - backgroundRate)}/.rule

Numeric real-world values

The total population is 270 millions and the lung cancer deaths are 152 thousands.
We assume statistical independance of exposure probabilities

backgroundFraction[backgroundRate_] := 
      backgroundRate (1 - exposureProbability[smoke]) (1 - 
            exposureProbability[radon]) (1 - 
            exposureProbability[asbestos]) 270/
          0.152

    
exposureProbability[smoke] = 0.45;
exposureProbability[radon] = 0.5;
exposureProbability[asbestos] = 0.05; 
exposureProbability[unknown] = 0.05;

Our information is given as upper and lower bound on the attributable risk.

commonKnowledge[backgroundRate_, 
        boundsList_] := 
{m[0] ==  
          backgroundFraction[
            backgroundRate], 

bound[
            attributableRisk[smoke], backgroundRate] == 
          boundsList[[1, 
            2]], 
lowerBound[
            attributableRisk[smoke], backgroundRate] == 
          boundsList[[1, 
            1]], 
 
bound[attributableRisk[radon], backgroundRate] ==  
          boundsList[[2, 
            2]], 
lowerBound[
            attributableRisk[radon], backgroundRate] == 
          boundsList[[2, 
            1]], 

bound[attributableRisk[asbestos], 
            backgroundRate] == 
          boundsList[[3, 
            2]], 
lowerBound[
            attributableRisk[asbestos], backgroundRate] == 
          boundsList[[3, 
            1]] , 

noMultipleSynergy[
          3]
}

expert1Bounds = {{0.7, 0.95}, {0.02, 0.21}, {0.01, 
            0.05}};
expert1BGRate  = 3. /100000;

Coherence check

The bounds avoid sure loss if the first five lines can be True, and coherent if all constraints can be satisfied  (neglecting the background rate).

constSet = 
        coherenceConstraints[
          Append[expert1Bounds, {loX, upX}]];
constSet // ColumnForm

Here are the bounds on X so that coherence is maintained (neglecting the background rate).

ConstrainedMin[upX - loX, 
        constSet, {upX, 
          loX}][[2]]
ConstrainedMax[upX - loX, 
        constSet, {upX, 
          loX}][[2]]

Results: Bounds on X

Here are most unspecific bounds on X so that coherence is maintained AND 3-way synergies are disallowed AND background rate is non-zero.

\(\(\(result1 = 
      Append[mostUnspecificExposure[
          commonKnowledge[expert1BGRate, expert1Bounds]], 
        backgroundRate -> expert1BGRate]\)\(
\)
    \)\)
resultSummary[result1]
resultDetail[result1] // ColumnForm

    
\(TableForm[Transpose[statTable[result1]], 
      TableHeadings -> {Range[
            15], {"s", \*""(1-(Subscript[r, 0])) m (\!\(and\+_\))"", \*""(1-(Subscript[r, 0])) bel (\!\(or\+_\))"", \*""(1-(Subscript[r, 0]))q (\!\(and\&_\))"", \*""(1-(Subscript[r, 0]))pl (\!\(or\&_\))""}}]\)