Manetees and bombers

• conditioning: Dependence on state
• everything is conditional
  • on data
  • on model
• Interaction: Influence of predictor conditional on other predictor

x’s effect on y - simple

Sometimes, we just want to investigate one variable’s effect on an outcome the estimand = smoking on lung cancer

flowchart LR
    X[X] --> Y[Y]

\(y ~ x\)

In session 11, we learned…

We cannot say age have and effect on HD, as this effect is strong in younger age; but the effect dies off at older age.

The truth is complicated - we need a simple answer!

Often, we don’t care about every possible value of age. We might be interested in a specific age group, or we might be interested in the average effect of all age groups.

Marginalizing the problem

flowchart LR

    X((X)) --> Y[Y]

    subgraph Marginalized covariates
        C1[Age]
        C2[Income]
        C3[Education]
        C4[Region]
        C5[...]
    end

    C1 --> Y
    C2 --> Y
    C3 --> Y
    C4 --> Y
    C5 --> Y

    style X stroke:#d62728,stroke-width:4px

How might we do that?

Three good way of doing this:

1. predicting the effect at a specific age
2. predicting the average effect of all participants’ age in the current study
3. predicting the average effect of any individual’s, considering we want to generalize to the population

How?

age2use <- 50

age50_HD_true <- as.numeric(
    coef(model_interact)[1] +               # intercept
        coef(model_interact)[2]*age2use +   # age
        coef(model_interact)[3]*1 +         # HD = true
        coef(model_interact)[4]*1*age2use   # HD = true and age
    )

age50_HD_false <- as.numeric(
    coef(model_interact)[1] +               # intercept
        coef(model_interact)[2]*age2use +   # age
        coef(model_interact)[3]*0 +         # HD = false
        coef(model_interact)[4]*0*age2use   # HD = false
        )

Using {marginaleffects}

model_interact |>
    marginaleffects::predictions(
        newdata = datagrid(Age = 50),
        variables = "HD"
        )

 Age Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
  50      161       1.57 102.4   <0.001 Inf   158    164
  50      142       2.15  65.9   <0.001 Inf   137    146

Type: response

2. predictions() will compute predictions
3. newdata to make a new dataset with datagrid
4. variables takes the variable we examine

Hypothesis test on the effect of HD at age 50

 Age Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
  50    -19.4       2.66 -7.29   <0.001 41.5 -24.6  -14.2

Term: HD
Type: response
Comparison: yes, heart disease - no heart disease

\(E(Y∣HD = yes,Age=50)−E(Y∣HD = no,Age=50)\)

Uncertainty

    Age 
52.7862 

         [,1]
[1,] 14.57242

Suppose your model is

\(η=β0​+β1​Age+β2​HD+β3​(Age×HD)\)

and you’re interested in the difference in predictions at Age = 50:

\(g(β)=y^​(HD=1,Age=50)−y^​(HD=0,Age=50)\)

The delta method says:

\(Var(g(β^))≈∇g(β^)TVar(β^)∇g(β^)\)

The trick is finding the gradient \(∇g\)

The variance is: \(Var(g)=L⊤VL\)

Imagine…

We do this for all combinations of age and HD we have in our data set?

We’d get the average effect of all participants’ age on HD in the current study

Ready for a lot of homework?

{marginaleffects} to the rescue

# A tibble: 6 × 2
    Age HD                
  <dbl> <fct>             
1    63 no heart disease  
2    67 yes, heart disease
3    67 yes, heart disease
4    37 no heart disease  
5    41 no heart disease  
6    56 no heart disease  

It’s easy!

model_interact |>
  predictions(
    newdata = hd_data, # can be left alone
    variables = list(
      HD = c("no heart disease", "yes, heart disease")
    )
  )

 Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
      147       2.25 65.4   <0.001   Inf   143    152
      143       2.76 51.9   <0.001   Inf   138    149
      143       2.76 51.9   <0.001   Inf   138    149
      175       2.92 59.9   <0.001   Inf   169    181
      171       2.40 71.2   <0.001   Inf   166    175
--- 596 rows omitted. See ?print.marginaleffects ---
      144       2.93 49.0   <0.001   Inf   138    149
      135       2.89 46.8   <0.001   Inf   129    141
      139       1.65 84.3   <0.001   Inf   136    142
      139       1.65 84.3   <0.001   Inf   136    142
      146       4.22 34.7   <0.001 871.6   138    154
Type: response

Now, we average all HD = 1 and HD = 0 separately

                 HD Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 no heart disease        156       1.55 101.2   <0.001 Inf   153    159
 yes, heart disease      140       1.71  81.9   <0.001 Inf   137    143

Type: response

Counterfactuals

This is the counterfactual estimate, where we predict the maxHR for all participants, had we observed them with and without HD.

How does our sample look?

# A tibble: 2 × 2
  HD                     n
  <fct>              <int>
1 no heart disease     164
2 yes, heart disease   139

How would we like our target population to be?

We could consider

{marginaleffects} can predict the target population

                 HD Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
 no heart disease        151       1.93 78.2   <0.001 Inf   147    154
 yes, heart disease      138       1.79 77.0   <0.001 Inf   135    142

Type: response

Hypothesis testing on the target population

 Hypothesis Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
      b1=b2     12.5       2.63 4.76   <0.001 19.0  7.37   17.7

Type: response

So this model

• estimates a 12.5 difference in maxHR between HD and no-HD in the age group 50-70 years old
• the estimate has a SE of 2.63 and 95%CI[7.37, 17.7]
• the z-statistic is 4.76 and the p-value is < 0,001

so we reject the null that HD has no effect on maxHR given the average effect of an age group between 50-70.