Quando considerare un numero come un conteggio

Aug 24 2020

Sto testando le differenze di gruppo nel numero di giorni in cui i partecipanti hanno utilizzato un farmaco nei 28 giorni precedenti. Questi sono i dati, ma ho difficoltà a decidere quale approccio utilizzare: regressione gaussiana standard o regressione binomiale aggregata. Ho fatto domande simili in precedenza su CV (ad esempio qui ) ma sono ancora un po 'insicuro.

Ho fornito il codice R per la replicabilità, ma ovviamente chiunque voglia intervenire - utente R o altro - è più che benvenuto.

df <- data.frame(group = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                 baseline = as.integer(c(28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 12, 28, 28, 28, 28, 28, 28, 24, 28, 28, 28, 28, 28, 28, 28, 28, 20, 28, 28, 24, 24, 28, 28, 28, 28, 28, 28, 28, 24, 28, 28, 28, 28, 28, 16, 28)),
                 outcome = as.integer(c(28, 0, 28, 0, 0, NA, NA, 16, 28, 10, 12, 0, 28, 12, 0, 0, 28, 8, 0, 28, 28, 0, 4, NA, NA, 0, NA, 28, NA, 20, 1, 3, 28, 26, NA, 0, 20, 16, 16, 0, NA, 3, 0, 1, 20, 0)),
                 coverage = 28)

groupè il trattamento ricevuto dai partecipanti; baselineil numero di giorni che hanno utilizzato nei 28 giorni precedenti l'inizio dello studio; outcomeil numero di giorni utilizzati durante lo studio di 28 giorni ( coverageè il numero di giorni nello studio).

Ecco le statistiche riassuntive:

library(tidyverse)

df %>%
  group_by(group) %>%
    drop_na(outcome) %>%
      summarise(mean = mean(outcome, na.rm = T),
                sd = sd(outcome, na.rm = T),
                median = median(outcome, na.rm = T),
                firstQuartile = quantile(outcome, probs = 0.25, na.rm = T),
                thirdQuartile = quantile(outcome, probs = 0.75, na.rm = T),
                tally = n()) 

# output
# group  mean    sd median firstQuartile thirdQuartile tally
# <dbl> <dbl> <dbl>  <int>         <dbl>         <dbl> <int>
#     0  10.7  11.3      3             0            20    17
#     1  12.3  12.3     10             0            28    21

E la distribuzione dei risultati in ogni gruppo

ggplot(df, aes(x = outcome, group = group)) + geom_histogram() + facet_wrap(~group) + scale_x_continuous(breaks = seq(0,28,7))

Come è tipico per i dati sull'uso di sostanze, i risultati sono distribuiti in modo abbastanza bimodale.

Quando analizzo il risultato, i giorni di regresso utilizzati, trattati come una variabile continua, sui giorni di trattamento groupe baselineutilizzati ...

summary(contMod <- lm(formula = outcome ~ group + baseline, 
                      data = df, 
                      na.action = na.exclude))

# output
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  17.7783    16.0047   1.111    0.274
# group         1.7020     3.9248   0.434    0.667
# baseline     -0.2660     0.5919  -0.449    0.656

Il modello non indica differenze significative tra i gruppi nella media dei giorni utilizzati durante il controllo per i giorni di riferimento utilizzati. Tuttavia, i residui del modello sono molto non normali:

hist(resid(contMod))

La trama quantile-quantile racconta la stessa storia

plot(contMod,2)

Quindi a me sembra che la regressione gaussiana standard potrebbe non essere appropriata per modellare questi dati.

Dato che i dati sono fondamentalmente conteggi interi delle occorrenze di un evento binario (utilizzato il giorno x vs non utilizzato il giorno x ) entro un numero noto di "prove" (28 giorni). Mi chiedevo se una regressione binomiale aggregata potesse essere un modo più appropriato per modellare i dati?

summary(contMod <- glm(formula = cbind(outcome, coverage-outcome) ~ group + baseline, 
                       data = df, 
                       family = binomial,
                       na.action = na.exclude))

# output
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)  
# (Intercept)  0.54711    0.50908   1.075   0.2825  
# group        0.25221    0.12634   1.996   0.0459 *
# baseline    -0.03866    0.01886  -2.050   0.0403 *

Ora la differenza di gruppo è significativa quando si controlla la linea di base.

Una differenza così drammatica nei risultati di due diversi modelli dello stesso è abbastanza sorprendente per me. Ovviamente sapevo che era possibile, ma non l'avevo mai incontrato prima in natura.

Quindi ho diverse domande per gli utenti intelligenti di CV

1. La regressione binomiale aggregata è un modo migliore per modellare questi dati data l'estrema non normalità sia del risultato che dei residui del modello? E se è opportuno l'ho fatto correttamente? E, anche se l'ho fatto correttamente, c'è un altro modo ancora migliore (non parametrico per esempio? Il test di Kruskal-Wallis ha kruskal.test(outcome ~ group, data = df)prodotto risultati simili al gaussiano,$\chi^2 = 0.07, p = 0.80$, ma non controlla la linea di base).

2. Come si interpreta l'output della regressione logistica aggregata? Se il risultato fosse un processo di Bernoulli, userei una semplice regressione logistica binaria e interpretare i risultati sarebbe semplice, esponenziale il coefficiente di gruppo e che rappresenta quanto sono maggiori le probabilità di utilizzo nel singolo giorno in questione nel 1gruppo rispetto al 0gruppo. Ma sono meno sicuro di come riportare il risultato dal binomio aggregato.

Aiuto molto apprezzato come sempre.

Risposte

2 NickCox Aug 25 2020 at 16:14

Stai facendo una domanda sui metodi qui, ma preferirei iniziare una risposta dai tuoi dati e da ciò che vuoi sapere.

Ecco una versione dei tuoi dati che può essere utile alle persone che non usano abitualmente R; le linee di apertura e chiusura sono specifiche per Stata, ma gli utenti di altri software dovrebbero essere in grado di adattarsi a seconda delle necessità. I periodi sono il codice generico di Stata per i mancati numeri e corrispondono a NA in R.

clear
input byte(id group baseline outcome coverage)
 1 1 28 28 28
 2 1 28  0 28
 3 1 28 28 28
 4 1 28  0 28
 5 1 28  0 28
 6 1 28  . 28
 7 1 28  . 28
 8 1 28 16 28
 9 1 28 28 28
10 1 28 10 28
11 1 12 12 28
12 1 28  0 28
13 1 28 28 28
14 1 28 12 28
15 1 28  0 28
16 1 28  0 28
17 1 28 28 28
18 1 24  8 28
19 1 28  0 28
20 1 28 28 28
21 1 28 28 28
22 1 28  0 28
23 1 28  4 28
24 1 28  . 28
25 0 28  . 28
26 0 28  0 28
27 0 20  . 28
28 0 28 28 28
29 0 28  . 28
30 0 24 20 28
31 0 24  1 28
32 0 28  3 28
33 0 28 28 28
34 0 28 26 28
35 0 28  . 28
36 0 28  0 28
37 0 28 20 28
38 0 28 16 28
39 0 24 16 28
40 0 28  0 28
41 0 28  . 28
42 0 28  3 28
43 0 28  0 28
44 0 28  1 28
45 0 16 20 28
46 0 28  0 28
end

Il nocciolo del problema è il confronto della media outcomeper due valori di group. Una distrazione è che baselinevaria e sembra essere più semplice almeno all'inizio ignorare i casi che non sono per 28 giorni baseline. Non è ovvio per me che l'aggiunta baselinecome predittore sia il modo migliore per regolare la variazione baseline; un'alternativa è scalare outcomefino a una frazione di baseline, ma anche questo potrebbe creare confusione.

Il confronto delle medie può naturalmente essere ripreso come un problema di regressione. Di seguito è riportato un grafico con la linea di regressione sovrapposta per la regressione di outcomeon groupper baseline28 giorni. Con questa semplificazione, la linea collega solo i due mezzi di gruppo.

Sto ruotando i tuoi istogrammi e trattando i dati per quello che sono, punti dati in un problema di regressione che confronta i mezzi. L'impilamento di risultati identici è solo una convenzione grafica e non influisce sui risultati di regressione.

Anche se ti riferisci alla "regressione gaussiana", la condizione ideale degli errori gaussiani o normali è l'aspetto meno importante della regressione lineare. Il recente testo di Gelman e amici

https://www.cambridge.org/core/books/regression-and-other-stories

sconsiglia anche i normali grafici quantili dei residui come perdita di tempo. Non andrei così lontano, ma è un punto di vista serio.

Uno sguardo al grafico e ai risultati della regressione indica una differenza di 2,9 giorni; la mia ipotesi è che una differenza di tale entità sarebbe clinicamente o scientificamente interessante, ma i risultati della regressione mostrano che il campione è troppo piccolo per confermarlo come significativo a livelli convenzionali. Se si è preoccupati che tale indicazione dipenda eccessivamente dall'assunzione implicita di errori normali, alcuni bootstrap di quei risultati di regressione implicano che una differenza di zero è ben all'interno di qualsiasi intervallo di confidenza per la differenza delle medie. Ritirarsi a Kruskal-Wallis mi sembra un consiglio di disperazione; perché usare la tecnologia degli anni '50 quando la tecnologia degli anni '70 (bootstrap) è disponibile e ti consente di concentrarti sulla differenza di mezzi che è di primo interesse?

In generale, è davvero una buona idea essere sensibili al fatto che i tuoi dati vengano conteggiati o misurati; pensare alle loro distribuzioni condizionali; e per notare se un risultato è necessariamente limitato. In questo caso particolare, questi semplici risultati di regressione implicano che poco importa ciò che si assume o è ciò che si presume o è l'ideale per i metodi utilizzati. La differenza tra le medie sembra interessante ma non è convenzionalmente significativa e tale indicazione è robusta a qualsiasi cosa tu faccia a titolo di analisi.

Tuttavia, se provo a far corrispondere la tua regressione binomiale, ma concentrandomi su baselineuguale a 28, trovo che sia sufficiente invertire la differenza in convenzionalmente significativa. All'inizio non ho capito perché ci sia una così grande differenza nell'indicazione.

Ma dovremmo preoccuparci di ciò che si presume sulle distribuzioni. Noto che i binomi non possono essere a forma di U. All'inizio dubitavo che fosse quello il problema, ma quel pensiero era viscerale, non logico. Se ripeti l'analisi con errori standard robusti (Eicker-Huber-White), il significato evapora.

In breve, applicando la regressione binomiale piuttosto che una semplice regressione stai sostituendo un'ipotesi distributiva che non morde - anche se sembra del tutto sbagliata - con un'ipotesi distributiva che morde! Questa è la mia diagnosi.

FWIW, l'uso dei giorni qui come numero intero è in parte naturale (le persone seguono ritmi giornalieri, a volte in modo rigido e talvolta in modo lasco) e in parte una convenzione (in linea di principio i dati potrebbero essere basati anche su ore del giorno, producendo giorni frazionari) .

Infine, il confronto dei mezzi non è l'unico gioco in città. Noto che nel gruppo 0 solo 2 su 13 ma nel gruppo 1 7 persone su 19 hanno riferito i 28 giorni completi. Queste differenze hanno naturalmente influenzato i mezzi, ma anche i dettagli possono essere importanti.

Nitty-gritty segue come output Stata. Le persone R si aspettano che siamo abbastanza intelligenti da decodificare l'output R se siamo abbastanza perversi da non usarlo (non usarlo di routine, nel mio caso) e restituisco il complimento. Il minimalismo dell'output R è ammirevole, tranne per il fatto che non mostrare la dimensione del campione utilizzata anche nel riepilogo predefinito mi lascia perplesso.