# Events/Trials Syntax vs. Single-Trial Mixed Model - Beta-Binomial Model

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;


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;


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?