0
$\begingroup$

I'm a touch confused about modeling an events/trials outcome when the Bernoulli trials are not independent, such as when a series of Bernoulli trials are observed from the same person. I cannot find a direct answer as most texts speak to events/trials as shown in this post. Further, this post is related, but does not directly suggest accounting for non-independence of trials within each person.

Below is a small simulation to show my issue. I simulate 10 independent Bernoulli trials for each of 1,000 people (think, whether or not a student attended school on a given day over a 10-day period). The number of events is simulated to follow a beta-binomial distribution. I also have a person-level treatment variable (Tx; constant within a person).

%LET t = 10; /*N trials*/

DATA EventTrial;
    CALL STREAMINIT(12);
    DO PersonID = 1 TO 1000;
        Trials = &t.;
        Tx = RAND("BERNOULLI",0.50);
        Hold = 1; /*Placeholder for PROC TRANSPOSE to be used later*/
        lp = 10 + 1*Tx; /*note no random intercept here*/
        prob = RAND("BETA",lp,2);
        Events = RAND("BINOMIAL",prob,Trials);
        OUTPUT;
    END;
RUN;

The beta-binomial model can be estimated in SAS via PROC FMM or PROC NLMIXED (see here; for NLMIXED see Dale McLerran's additions in the comment section). I'll use NLMIXED because I can add random effects (later).

TITLE "NLMIXED: Beta-Binomial Model";
PROC NLMIXED DATA=EventTrial TECH=QUANEW FCONV=1E-50 GCONV=1E-50;
    PARMS b0 1.6 b1 0.1 log_a 2.5;
    x=Events;
    Ntrials=Trials;
    eta = b0 + b1*Tx;
    pi = 1 / (1 + 1/exp(eta));
    a = exp(log_a);
    phi = a/pi;
    mu = nTrials*pi;
    ll = lgamma(nTrials+1) - lgamma(x+1) - lgamma(nTrials - x + 1)
         + lgamma(phi) - lgamma(nTrials+phi) + lgamma(x+a)
         + lgamma(nTrials - x + phi - a) - lgamma(a)
         - lgamma(phi-a);
    MODEL ll ~ GENERAL(ll);
    ESTIMATE "Intercept" b0 DF=998;
    ESTIMATE "Tx" b1 DF=998;
    *ESTIMATE "a" a;
    *ESTIMATE "b" a*(1-(1 / (1 + 1/exp(b0 + b1)))) / (1 / (1 + 1/exp(b0 + b1)));
    ESTIMATE "Scale" a / (1 / (1 + 1/exp(b0 + b1)));
RUN;

Beta-Binomial Model

The model above (incorrectly?) assumes that all Bernoulli trials were independent, whether from the same person or not. If, however, I transpose the data to the single-trial level, I can add random effects into the NLMIXED code; specifically, a random intercept to account for the nesting of single-trial outcomes within the same person.

DATA Wide; SET EventTrial;
    BY PersonID;
    ARRAY T{&t.} T1-T&t.;
    DO i = 1 TO &t.;
        IF (i <= Events) THEN T{i} = 1; ELSE T{i} = 0;
    END;
    DROP i;
RUN;

PROC TRANSPOSE DATA=Wide OUT=Long(RENAME=(COL1=Event Hold=Trial) DROP=_NAME_);
    BY PersonID Tx Hold;
    VAR T1-T&t.;
RUN;

TITLE "NLMIXED: Mixed-Effects Beta-Binomial Model";
PROC NLMIXED DATA=Long TECH=QUANEW FCONV=1E-50 GCONV=1E-50;
    PARMS b0 1.6 b1 0.1 log_a 2.5 u0var 0.60;
    x=Event;
    Ntrials=Trial;
    eta = (b0 + u0i) + b1*Tx;
    pi = 1 / (1 + 1/exp(eta));
    a = exp(log_a);
    phi = a/pi;
    mu = nTrials*pi;
    ll = lgamma(nTrials+1) - lgamma(x+1) - lgamma(nTrials - x + 1)
         + lgamma(phi) - lgamma(nTrials+phi) + lgamma(x+a)
         + lgamma(nTrials - x + phi - a) - lgamma(a)
         - lgamma(phi-a);
    MODEL ll ~ GENERAL(ll);
    RANDOM u0i ~ NORMAL(0,u0var) SUBJECT=PersonID;
    ESTIMATE "Intercept" b0 DF=998;
    ESTIMATE "Tx" b1 DF=998;
    *ESTIMATE "a" a;
    *ESTIMATE "b" a*(1-(1 / (1 + 1/exp(b0 + b1)))) / (1 / (1 + 1/exp(b0 + b1)));
    ESTIMATE "Random Intercept Variance" u0var;
    ESTIMATE "Scale" a / (1 / (1 + 1/exp(b0 + b1)));
RUN;

enter image description here

At first glance, including random intercept does not seem to make a big difference; though, that is an empirical question that I could test by simulating out to 1,000 reps or so. It is my understanding that the model-scale estimates remain comparable between these two models because the unit-specific-ness resulting from logit link and random effects only applies to the data-scale (i.e., probability) estimates.

Anyways, my overarching question is whether non-independence when using events/trials syntax something that I should be worrying about?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.