Binomial analysis

[aska_ujita]

Hello, I have two treatments (treated and control) with 20 animals each treatment.

And I have a reactivity score (0- no reactive and 1-reactive) from 10 days (day 1, 3, 5, 7, ..., 17, 19).

I would like to compare A and B by day just for escore 1.

So the objective is see if the treated group had less reactivity than group control on days 1 till 19.

They told me to use proc genmod or glimmix.

I have measure repeated in time: I have 10 information (1 information by day collected) per animal.

Can somebody help me, please?

[Rajesh3 · Accepted solution]

proc genmod data=fake descending
order=data;
class treatment animal;
model reactivity=treatment day/ d=bin link=logit;
repeated subject=animal / type=cs;
estimate "treatment" treatment 1 -1 / exp;
estimate "day" day / exp;
run;

The above code fits a GEE model with the binomial distribution.

Thank you,

Rajesh.

[aska_ujita]

Hi Rajesh, Thank you for your fast reply! Very kind of you, thank you very much for helping me.

So I tried to do it and it works for few, but I got errors in others.

I put it like this:

proc genmod data=C descending order=data;
class group cow dl po;
model react=group dl po/ d=bin link=logit;
repeated subject=cow/ type=cs;
estimate "group" group1 -1 / exp;
estimate "dl" dl1 -1/ exp;
run;

group is control and treat

dl is the day of lactation (1, 3, 5, 7, 9, 11, 13, 15, 17 and 19)

po is parturition order

react is the reactivity scored as 0 (without reactivity) and 1 (with reactivity).

In my table, I will just compare the score 1 (with reactivity). As example below:

DL	Reactivity
Group Treated	Group Control	P value
01	50	80
03	30	65
05	20	35

I would like to have the P value.

Sometimes the difference between groups is large and the P value isn't significative.

And sometimes have errors like this: 

NOTE: PROC GENMOD is modeling the probability that react='1'.
WARNING: The negative of the Hessian is not positive definite. The convergence is questionable.
WARNING: The procedure is continuing but the validity of the model fit is questionable.
WARNING: The specified model did not converge.
WARNING: Negative of Hessian not positive definite.
WARNING: The generalized Hessian matrix is not positive definite. Iteration will be terminated.
ERROR: Error in parameter estimate covariance computation.
ERROR: Error in estimation routine.

I'm bit new in this... and a desesperate student.

Please, if you can help me, I am really glad.

Thank you again.

Aska.

[Rajesh3]

proc genmod data=C descending order=data;
class group cow po;
model react=group dl po/ d=bin link=logit;
repeated subject=cow/ type=cs;
estimate "group" group1 -1 / exp;
estimate "dl" dl/ exp;
run;

The 'dl' variable need not be mentioned in the class statement. If you do not mention it in the class statement, it will be treated as a continuous variable. Also, if you do not need to control for 'po' variable, please remove it from the model. Please check if the model converges with the above code.

Based on the error that you are getting, you have specified a more complex model than your data can support.

Hope this helps!

Thank you,

Rajesh.

[aska_ujita]

Thank you for you reply and help!

It works now, thank you for you explanation.

But, abusing again your goodwill and knowledge, now the results havent the information of estimate, SE, Z and Pr>Z values for each DL...

Thank you again!

Aska.

[Rajesh3]

It will not give estimates for each time point separately. Now the interpretation becomes "change in reactivity for a day increase in time". If you want to interpret it for 2 days increase multiply the estimate for 'dl' by 2 and exponentiate it to get the appropriate change. Else, replace the relevant estimate statement with the below code.

estimate 'dl'   dl 2/exp;

Thanks,

Rajesh.

[aska_ujita]

Hello Rajesh, thank you.

This is the last thing, I promise. When I put: estimate "dl" dl 1 -1/ exp;

What this "1 -1" means?

Thank you!

Aska.

[Rajesh3]

If you included 'dl' in the class statement, the below SAS code will apply.

estimate "dl" dl 1 -1/ exp;

It means for a change in day from 1 to 3/1 to 5 (this will change based on the reference chosen but the concept is this), the change (increase/decrease) in the odds of having reactivity=1.

If you did not include 'dl' in the class statement, the estimate statement should be:

estimate "dl" dl 1/ exp;

It means, on average, for a unit increase in day, the change (increase/decrease) in the odds of having reactivity=1.

Thank you,

Rajesh.

[aska_ujita]

Thank you Rajesh again.

So that mean that is comparing day with day by genmod procedure.

And probably I have to change my PROC because what I would like to do is compare the difference between groups for each day and not compare between days for presence of reactivity.

Thank you for helping me!!! It really helped me!!

Aska