Fitting a Nested Random Effects Logistic Regression Model using GENLINMIXED

[Ryan]

All,

 

I was interested in determining if one could fit a nested random effects logistic regression model by using two RANDOM statements within the GENLINMIXED procedure. I generated data from a model which included nested random effects along with two fixed effects predictors, one of which is at the "person" level while the other is at the "clinic" level. Data simulation code and the GENLINMIXED code are below. The answer is in the affirmative. Moreover, the parameters were actually recovered better than expected. Still, simulation studies suggest that pseudo-likelihood estimation methods tend to yield biased estimates. I hope that one day soon the GENLINMIXED procedure is capable of employing integral approximation methods. Anyway, I thought I'd share the code. Hope this is of interest to others.

 

Best,

 

Ryan

--

 

*Random Effects Logistic Regression Simulation.
 set seed 65923454.
 
new file.
 inp pro.
 
compute clinic_ID = -99.
compute #clinic_rand_eff = -99.
compute clinic_SES = -99.
compute person_ID = -99.
compute #person_rand_eff = -99.
compute gender = -99.
compute observation = -99.
 
leave clinic_ID to observation.
 
loop clinic_ID = 1 to 250.
   compute #clinic_rand_eff = sqrt(0.30)*rv.normal(0,1).
   compute clinic_SES = rv.normal(0,1).

loop person_ID= 1 to 50.
   compute #person_rand_eff = sqrt(0.50)*rv.normal(0,1).
   compute gender = rv.bernoulli(.50).
 
  loop observation = 1 to 10.
   compute y = rv.bernoulli(1 / (1+exp(-(-3.0 + 2.0*(gender=0) + 1.5*clinic_SES + #clinic_rand_eff + #person_rand_eff)))).

     end case.
    end loop.
   end loop.
  end loop.
 end file.
 end inp pro.
 exe.

 

VARIABLE LEVEL  clinic_ID(NOMINAL).
VARIABLE LEVEL  person_ID(NOMINAL).
EXECUTE.

 

GENLINMIXED
  /DATA_STRUCTURE SUBJECTS=person_ID*clinic_ID
  /FIELDS TARGET=y TRIALS=NONE OFFSET=NONE
  /TARGET_OPTIONS REFERENCE=0 DISTRIBUTION=BINOMIAL LINK=LOGIT
  /FIXED  EFFECTS=gender clinic_SES USE_INTERCEPT=TRUE
  /RANDOM USE_INTERCEPT=TRUE SUBJECTS=clinic_ID COVARIANCE_TYPE=VARIANCE_COMPONENTS
  /RANDOM USE_INTERCEPT=TRUE SUBJECTS=person_ID*clinic_ID COVARIANCE_TYPE=VARIANCE_COMPONENTS.

 

 

[Ryan]

To be complete, here are the parameter estimates from GENLINMIXED:

 

1. clinic_SES = 1.341

2. gender = 1.851

3. clinic_ID VC = .240

4. person_ID VC = .438

5. fixed intercept = -2.767

 

...and here are the parameter estimates from SAS GLIMMIX using the Laplace estimation method:

 

1. clinic_SES = 1.432

2. gender = 1.976

3. clinic_ID VC = .259

4. person_ID VC = .516

5. fixed intercept = -2.954

 

As a reminder, the simulation parameters were:

 

1. clinic_SES = 1.50

2. gender = 2.00

3. clinic_ID VC = .300

4. person_ID VC = .500

5. fixed intercept = -3.00

 

Ryan
 

On Wed, Apr 10, 2013 at 10:56 PM, R B <[hidden email]> wrote:

All,

 

I was interested in determining if one could fit a nested random effects logistic regression model by using two RANDOM statements within the GENLINMIXED procedure. I generated data from a model which included nested random effects along with two fixed effects predictors, one of which is at the "person" level while the other is at the "clinic" level. Data simulation code and the GENLINMIXED code are below. The answer is in the affirmative. Moreover, the parameters were actually recovered better than expected. Still, simulation studies suggest that pseudo-likelihood estimation methods tend to yield biased estimates. I hope that one day soon the GENLINMIXED procedure is capable of employing integral approximation methods. Anyway, I thought I'd share the code. Hope this is of interest to others.

 

Best,

 

Ryan

--

 

*Random Effects Logistic Regression Simulation.
 set seed 65923454.
 
new file.
 inp pro.
 
compute clinic_ID = -99.
compute #clinic_rand_eff = -99.
compute clinic_SES = -99.
compute person_ID = -99.
compute #person_rand_eff = -99.
compute gender = -99.
compute observation = -99.
 
leave clinic_ID to observation.
 
loop clinic_ID = 1 to 250.
   compute #clinic_rand_eff = sqrt(0.30)*rv.normal(0,1).
   compute clinic_SES = rv.normal(0,1).

loop person_ID= 1 to 50.
   compute #person_rand_eff = sqrt(0.50)*rv.normal(0,1).
   compute gender = rv.bernoulli(.50).
 
  loop observation = 1 to 10.
   compute y = rv.bernoulli(1 / (1+exp(-(-3.0 + 2.0*(gender=0) + 1.5*clinic_SES + #clinic_rand_eff + #person_rand_eff)))).

     end case.
    end loop.
   end loop.
  end loop.
 end file.
 end inp pro.
 exe.

 

VARIABLE LEVEL  clinic_ID(NOMINAL).
VARIABLE LEVEL  person_ID(NOMINAL).
EXECUTE.

 

GENLINMIXED
  /DATA_STRUCTURE SUBJECTS=person_ID*clinic_ID
  /FIELDS TARGET=y TRIALS=NONE OFFSET=NONE
  /TARGET_OPTIONS REFERENCE=0 DISTRIBUTION=BINOMIAL LINK=LOGIT
  /FIXED  EFFECTS=gender clinic_SES USE_INTERCEPT=TRUE
  /RANDOM USE_INTERCEPT=TRUE SUBJECTS=clinic_ID COVARIANCE_TYPE=VARIANCE_COMPONENTS
  /RANDOM USE_INTERCEPT=TRUE SUBJECTS=person_ID*clinic_ID COVARIANCE_TYPE=VARIANCE_COMPONENTS.