Paradoxul lui Simpson într -o regresie logistică

Exemplu sintetic: încercare de pierdere în greutate

Aceasta este o declarație a problemei în care ne-am aștepta ca coeficienții unei regresii logistice să fie non-negative (cu excepția interceptării).

Luați în considerare un proces care testează eficacitatea unui regim alimentar specific (să zicem post intermitent 16/8, pe care îl vom numi ifasting) și un regim de exerciții specifice (o plimbare rapidă de 30 de minute în fiecare zi, pe care doar o vom suna exercise) Obiectivul („succes”) este de a pierde cel puțin cinci kilograme până la sfârșitul perioadei de încercare. Am creat trei grupuri de tratament, după cum urmează:

• 200 de subiecți încearcă să exercite singuri
• 300 de subiecți încearcă IFasting singur
• 300 de subiecți încearcă IFASTING plus exercițiul

Înainte de proces, toți subiecții au condus stiluri de viață destul de sedentare și nu au fost dieta în niciun mod formal.

Pentru acești subiecți, s -ar putea aștepta în mod rezonabil ca nici exercițiul și nici ifasting -ul Mai puțin succes pentru a pierde în greutate decât a nu face nimic. De asemenea, s -ar aștepta în mod rezonabil că IFASTING plus exercițiul nu ar trebui să facă mai rău decât să facă unul singur. Prin urmare, modelarea rezultatelor unui astfel de experiment ca regresie logistică ar trebui să conducă la un model în care coeficienții şi Sunt ambele non-negative, deoarece orice tratament ar trebui să crească (sau cel puțin, să nu scadă) șansele ca subiectul să piardă în greutate.

Să arătăm un exemplu în care așteptările noastre nu sunt îndeplinite. Cel mai simplu mod de a face acest lucru este de a genera un set de date care are paradoxul lui Simpson ascuns în el.

În primul rând, să încărcăm pachetele de care avem nevoie.

library(poorman) # or dplyr
library(ggplot2)
library(kableExtra)
library(WVPlots)

Iată o funcție care va genera un subset specific de date, după cum este necesar.

# ifasting: 1 if this group fasted, else 0
# exercise: 1 if this group exercised, else 0
# total: total number of subjects in this group
# successes:  number of subjects who successfully lost weight
# label: label for the group.
generate_samples = function(ifasting, exercise, total, successes, label) {
  failures = total-successes
  data.frame(ifasting = ifasting,
             exercise = exercise,
             success = c(rep(1, successes), rep(0, failures)),
             label=label)
}

Efecte ale populației ascunse, nemodelate

Să presupunem acum că, necunoscut cercetătorilor, există două sub -populații în rândul subiecților. Aceste două sub -populații răspund diferit la exercițiile fizice și la regimurile de post intermitente. Deci, vom genera datele noastre sintetice de către populație.

popA = data.table::rbindlist(list(
  generate_samples(ifasting=0, exercise=1, total=100, successes=2, "Population A"),
  generate_samples(ifasting=1, exercise=0, total=200, successes=160, "Population A"),
  generate_samples(ifasting=1, exercise=1, total=100, successes=90, "Population A")
))

#| code-fold: true
#| code-summary: "Show the code"
popA$treatment = with(popA, ifelse(ifasting, ifelse(exercise, 'both', 'ifast alone'),
                                   ifelse(exercise, 'exercise alone', 'none')))
popA$treatment = factor(popA$treatment,
                            levels = c('none', 'exercise alone', 'ifast alone', 'both'))
popA$outcome = factor(with(popA, ifelse(success==1, "success", "failure")),
                      levels = c('success', 'failure'))

pvals = c(success="darkblue", failure="gray")

ggplot(popA, aes(x=treatment, color=outcome, fill=outcome)) + 
  geom_bar() +
  scale_fill_manual(values=pvals) + 
  scale_color_manual(values=pvals) +
  ggtitle("Outcomes for population A") + 
  theme(legend.position="bottom")

df1 = popA |>
  group_by(treatment) |>
  summarize(success_rate=mean(success)) |>
  ungroup() 

df2 = popA |>
  summarize(success_rate = mean(success)) |>
  mutate(treatment = "overall") 

rbind(df1, df2) |>
  knitr::kable(digits=3, caption = "Success rates, Population A")

popB = data.table::rbindlist(list(
  generate_samples(ifasting=0, exercise=1, total=100, successes=1, "Population B"),
  generate_samples(ifasting=1, exercise=0, total=100, successes=40, "Population B"),
  generate_samples(ifasting=1, exercise=1, total=200, successes=100, "Population B")
))


#| code-fold: true
#| code-summary: "Show the code"
popB$treatment = with(popB, ifelse(ifasting, ifelse(exercise, 'both', 'ifast alone'),
                                   ifelse(exercise, 'exercise alone', 'none')))
popB$treatment = factor(popB$treatment,
                            levels = c('none', 'exercise alone', 'ifast alone', 'both'))
popB$outcome = factor(with(popB, ifelse(success==1, "success", "failure")),
                      levels = c('success', 'failure'))



pvals = c(success="darkblue", failure="gray")

ggplot(popB, aes(x=treatment, color=outcome, fill=outcome)) + 
  geom_bar() +
  scale_fill_manual(values=pvals) + 
  scale_color_manual(values=pvals) +
  ggtitle("Outcomes for population B") + 
  theme(legend.position="bottom")

df1 = popB |>
  group_by(treatment) |>
  summarize(success_rate=mean(success)) |>
  ungroup() 

df2 = popB |>
  summarize(success_rate = mean(success)) |>
  mutate(treatment = "overall") 

rbind(df1, df2) |>
  knitr::kable(digits=3, caption = "Success rates, Population B")

bothpops = rbind(popA, popB)
palette = c('Population A' = '#1b9e77', 'Population B' = '#d95f02')
ggplot(bothpops, aes(x=treatment, color=label, fill=label)) + 
  geom_bar() + 
  scale_fill_manual(values=palette) + 
  scale_color_manual(values=palette) +
  ggtitle("Treatment Group Sizes, with population labels") + 
  theme(legend.position="bottom")

Modelare

Acum să ne potrivim unui model de regresie logistică pentru a încerca să deducem efectele diferitelor tratamente asupra pierderii în greutate. O vom face mai întâi pe întreaga populație, deoarece aceasta a fost sarcina inițială.

tab_coeff = function(model, caption) {
  coeff = summary(model)$coefficients(, c(1, 4)) |>
    as.data.frame()

  colnames(coeff) = c('Estimate', 'pval')

  # using cell_spec below breaks the digits setting
  # (because of course it does) so round the numbers first.
  coeff = coeff |>
    mutate(Estimate = as.numeric(formatC(Estimate, format="f", digits=3)),
           pval = as.numeric(format(pval, format="g", digits=3)))

  coeff = coeff |>
    mutate(Estimate = cell_spec(Estimate, color=ifelse(Estimate < 0, "red", "black")),
           pval = cell_spec(pval, color=ifelse(pval < 0.05, "darkblue", "darkgray")))

  knitr::kable(coeff, caption=caption)

}

bothpops = rbind(popA, popB)
mAll = glm(success ~ ifasting + exercise, data=bothpops, family=binomial)

tab_coeff(mAll, "Model coefficients, whole population")

Postul intermitent are un coeficient pozitiv, ceea ce înseamnă că postul intermitent este corelat pozitiv cu succesul pierderii în greutate. Dar exercițiul are un negativ coeficientul, implicând exercițiul este corelat negativ cu pierderea în greutate și că a face ambele împreună va fi Mai puțin reușit decât postul intermitent singur!

Și, într -adevăr, dacă ne uităm la rezumatele brute, vom vedea că datele suportă aceste inferențe.

df1 = bothpops |>
  group_by(treatment) |>
  summarize(success_rate=mean(success)) |>
  ungroup() 

df2 = bothpops |>
  summarize(success_rate = mean(success)) |>
  mutate(treatment = "overall") 

rbind(df1, df2) |>
  knitr::kable(digits=3, caption = "Success rates, entire population")

Acesta este un exemplu al modului în care paradoxul lui Simpson s -ar putea manifesta într -un model de regresie logistică și se datorează variabilei de confuzie nemodificate, tipul de populație. Acest lucru, plus un anumit ghinion în dimensiunile relative ale grupurilor de tratament în ceea ce privește tipul de populație, duce la rezultatele de mai sus, contra intuitive.

Rețineți că am raportat valori p și, în acest caz, coeficientul pentru exerciții fizice este nesemnificativ (la ), implicând că exercițiile fizice nu pot avea niciun efect notabil asupra pierderii în greutate. Cu toate acestea, este posibil să vedem în continuare coeficienți cu semne contra intuitive în populații arbitrar mari, care pot părea apoi semnificative.

mA = glm(success ~ ifasting + exercise, data=popA, family=binomial)

tab_coeff(mA, "Model coefficients, Population A")

mB = glm(success ~ ifasting + exercise, data=popB, family=binomial)

tab_coeff(mB, "Model coefficients, Population B")

mAllplus = glm(success ~ ifasting + exercise + label, data=bothpops, family=binomial)

#| code-fold: true
#| code-summary: "Show the code"
tab_coeff(mAllplus, "Model coefficients, whole population with population labels")

Concluzie

Paradoxul lui Simpson este un caz în care inferența modelului pare să contrazică cunoștințele de domeniu. De obicei, este doar un simptom al unei combinații de prejudecăți variabile omise, studii dezechilibrate sau specificații cauzale greșite (sau poate că nu există o astfel de specificație). Dacă aveți grijă să căutați acest efect, acesta poate oferi, de fapt, indicii pentru o analiză mai bună.

Nina Zumel este un om de știință de date cu sediul în San Francisco, cu peste 20 de ani de experiență în învățarea automată, statistici și analize. Ea este co-fondatorul firmei de consultanță în domeniul științei de date Win-Vector LLC și (cu John Mount) coautorul științei datelor practice cu R, aflat acum în a doua ediție.