BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
aska_ujita
Obsidian | Level 7

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?

 

 

1 ACCEPTED SOLUTION

Accepted Solutions
Rajesh3
Obsidian | Level 7
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.

View solution in original post

8 REPLIES 8
Rajesh3
Obsidian | Level 7
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
Obsidian | Level 7

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
Obsidian | Level 7
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
Obsidian | Level 7

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
Obsidian | Level 7

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
Obsidian | Level 7

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
Obsidian | Level 7

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
Obsidian | Level 7

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

sas-innovate-2024.png

Available on demand!

Missed SAS Innovate Las Vegas? Watch all the action for free! View the keynotes, general sessions and 22 breakouts on demand.

 

Register now!

Mastering the WHERE Clause in PROC SQL

SAS' Charu Shankar shares her PROC SQL expertise by showing you how to master the WHERE clause using real winter weather data.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 8 replies
  • 890 views
  • 0 likes
  • 2 in conversation