Hello everyone,
I have a question about the necessary scripts to make a recruitment analysis in avian studies.
What script could you recommend me?
Would it be correct?
Proc Logistic data=Birds descending;
Class treatment Sex;
model Presence_2011 (Event='1')= treatment Sex treatment*Sex weight;
run;
Where for example Variables = Presence in 2011, 2012, 2013, sex, treatment, weight of chick ...
Would it be necessary to add more commands?
Thank you very much for your time
Jaime
So you will have to do separate analyses for each year, and rely on visual inspection to compare between years and globally. It appears that presence_global = max(presence_2011, presence_2012, presence_2013), so an analysis that puts it in with the year values may be biased. However, you may wish to compare any of the given years to the global value. To do that, you could put your data in long format, and have a CLASS variable that reflected year/global status.
data longbirds;
set birds;
presence=presence_2011; indicator='2011 '; output;
presence=presence_2012, indicator='2012 '; output;
presence=presence_2013; indicator='2013 '; output;
presence=presence_global; indicator='Global';output;
keep bird treatment sex weight indicator presence;
run;
Proc logistic data=longBirds descending;
Class treatment Sex indicator;
model presence (Event='1')= treatment Sex treatment*Sex indicator indicator*treatment indicator*Sex indicator*treatment*sex weight;
lsmeans treatment Sex treatment*Sex indicator indicator*treatment indicator*Sex indicator*treatment*sex /ilink oddsratio;
<possibly include ODDSRATIO or LSMESTIMATE statements to get comparisons to where you need to have them>
run;
Steve Denham
This would be a good first model. You will probably want to add an ODDSRATIO statement, an LSMEANS statement and an LSMESTIMATE statement, to get estimates and test for differences between them..
My question is how do you wish to deal with 'year' in this data? Do you wish to compare between years, or are estimates for each year separately sufficient?
Steve Denham
The years are presence (1) and absence (0) of birds. I really want to know if treatment had effects in the bird recruitment (for each year and in global). For example, I would like say that treatment 1 increased the recruitment in 2011 relative to control.
It is part of my database:
Bird | Treatment | Sex | Presence_2011 | Presence_2012 | Presence_2013 | Presence_Global | weight |
353701 | 0 | 2 | 0 | 0 | 0 | 0 | 56 |
353702 | 0 | 2 | 1 | 1 | 1 | 1 | 39 |
353703 | 0 | 2 | 1 | 0 | 0 | 1 | 54 |
353704 | 0 | 2 | 0 | 0 | 1 | 1 | 63 |
353705 | 1 | 1 | 0 | 0 | 1 | 1 | 60 |
353706 | 1 | 1 | 1 | 1 | 0 | 1 | 45 |
353707 | 1 | 1 | 0 | 1 | 0 | 1 | 52 |
353708 | 1 | 1 | 0 | 0 | 0 | 0 | 49 |
353709 | 1 | 1 | 0 | 0 | 1 | 1 | 38 |
Thank you very much
All the best;
Jaime Muriel
So you will have to do separate analyses for each year, and rely on visual inspection to compare between years and globally. It appears that presence_global = max(presence_2011, presence_2012, presence_2013), so an analysis that puts it in with the year values may be biased. However, you may wish to compare any of the given years to the global value. To do that, you could put your data in long format, and have a CLASS variable that reflected year/global status.
data longbirds;
set birds;
presence=presence_2011; indicator='2011 '; output;
presence=presence_2012, indicator='2012 '; output;
presence=presence_2013; indicator='2013 '; output;
presence=presence_global; indicator='Global';output;
keep bird treatment sex weight indicator presence;
run;
Proc logistic data=longBirds descending;
Class treatment Sex indicator;
model presence (Event='1')= treatment Sex treatment*Sex indicator indicator*treatment indicator*Sex indicator*treatment*sex weight;
lsmeans treatment Sex treatment*Sex indicator indicator*treatment indicator*Sex indicator*treatment*sex /ilink oddsratio;
<possibly include ODDSRATIO or LSMESTIMATE statements to get comparisons to where you need to have them>
run;
Steve Denham
Thank you very much for your useful assistance!!
All the best,
Jaime
Your scripts were very useful. However, I think something changes the sign that directs the results when comparing the results obtained with the original data... Could it be possible? Should I add or remove a word in the script?
Thank you very much again for your help,
Best wishes,
Jaime
First I would try changing to (Event='0') in the model statement, to see if that changes the sign. If you are looking at the LSMEANS, remember that the estimate is on the logit scale, and so might well be negative. Look at the mean value, making sure you used the ILINK option. It should be on the original scale as a proportion.
Steve Denham
Message was edited by: Steve Denham
If I change the 1 by 0 in the model statement, the result is the same. Instead of considering the recruitment, calculates the probability of not recruiting. But I mean I found the opposite effect of treatment if I make the graphics from raw data, that if I perform them with SAS. Hence I ask myself if a new word was needed or just left over a word in the script.
Thank you very much,
Best wishes,
Jaime
Be very careful about your interpretation, because of the interaction. For instance, if there is a significant interaction, looking at main effects is really futile, as the effect of the first variable differs at various levels of the second variable. A single plot is unlikely to capture this--you would need plots of the first variable at all levels of the second. Can you share your output, so that we might check this assumption?
Steve Denham
Proc logistic data=Vips descending;
Class treat Sex event year;
model presence (Event='1')= treat Sex treat*Sex year treat*weight event*treat weight*event weight*treat*event sex*event sex*event*treat weight event;
lsmeans treat Sex treat*Sex year treat*Sex/ilink oddsratio;
run;
Sistema SAS |
Procedimiento LOGISTIC
Conjunto de datos | WORK.VIPS |
---|---|
Variable de respuesta | presence |
Número de niveles de respuesta | 2 |
Modelo | logit binario |
Técnica de optimización | Puntuación de Fisher |
Número de observaciones leí | 3330 |
---|---|
Número de observaciones usa | 3213 |
1 | 1 | 463 |
---|---|---|
2 | 0 | 2750 |
La probabilidad modelada es presence =1. |
Note: | 117 observations were deleted due to missing values for the response or explanatory variables. |
Treat | 0 | 1 | |
---|---|---|---|
1 | -1 | ||
Sex | 1 | 1 | |
2 | -1 | ||
event | 1 | 1 | 0 |
2 | 0 | 1 | |
3 | -1 | -1 | |
year | 2012 | 1 | 0 |
2013 | 0 | 1 | |
2014 | -1 | -1 |
Criterio de convergencia (GCONV=1E-8) satisfecho. |
AIC | 2651.701 | 2517.243 |
---|---|---|
SC | 2657.776 | 2638.742 |
-2 LOG L | 2649.701 | 2477.243 |
Ratio de verosim | 172.4584 | 19 | <.0001 |
---|---|---|---|
Puntuación | 160.0191 | 19 | <.0001 |
Wald | 148.2303 | 19 | <.0001 |
Treat | 1 | 13.5012 | 0.0002 |
---|---|---|---|
Sex | 1 | 5.2251 | 0.0223 |
Treat*Sex | 1 | 0.3069 | 0.5796 |
year | 2 | 85.8229 | <.0001 |
weight*Treat | 1 | 13.6287 | 0.0002 |
Treat*event | 2 | 18.1990 | 0.0001 |
weight*event | 2 | 7.5996 | 0.0224 |
weight*treat*event | 2 | 16.6071 | 0.0002 |
Sex*event | 2 | 1.4899 | 0.4748 |
Treat*Sex*weight | 2 | 9.9884 | 0.0068 |
weight | 1 | 16.2378 | <.0001 |
event | 2 | 7.5174 | 0.0233 |
Intercept | 1 | -4.3347 | 0.6141 | 49.8221 | <.0001 | |||
---|---|---|---|---|---|---|---|---|
Treat | 0 | 1 | -2.2545 | 0.6136 | 13.5012 | 0.0002 | ||
Sex | 1 | 1 | 0.1430 | 0.0626 | 5.2251 | 0.0223 | ||
Treat*Sex | 0 | 1 | 1 | -0.0346 | 0.0625 | 0.3069 | 0.5796 | |
year | 2012 | 1 | 0.6182 | 0.0699 | 78.1565 | <.0001 | ||
year | 2013 | 1 | -0.00655 | 0.0759 | 0.0075 | 0.9312 | ||
weight*Treat | 0 | 1 | 0.0350 | 0.00948 | 13.6287 | 0.0002 | ||
treat*event | 0 | 1 | 1 | 0.1580 | 0.7670 | 0.0424 | 0.8368 | |
treat*event | 0 | 2 | 1 | -3.5230 | 1.0052 | 12.2831 | 0.0005 | |
weight*event | 1 | 1 | -0.0236 | 0.0116 | 4.1426 | 0.0418 | ||
weight*event | 2 | 1 | -0.00398 | 0.0156 | 0.0647 | 0.7993 | ||
weight*treat*event | 0 | 1 | 1 | -0.00515 | 0.0116 | 0.1972 | 0.6570 | |
weight*treat*event | 0 | 2 | 1 | 0.0547 | 0.0156 | 12.2368 | 0.0005 | |
Sex*event | 1 | 1 | 1 | -0.0194 | 0.0746 | 0.0678 | 0.7945 | |
Sex*event | 1 | 2 | 1 | 0.1158 | 0.1027 | 1.2701 | 0.2597 | |
Trt1*Sex*event | 0 | 1 | 1 | 1 | -0.0274 | 0.0746 | 0.1350 | 0.7133 |
Trt1*Sex*event | 0 | 1 | 2 | 1 | -0.2429 | 0.1027 | 5.5896 | 0.0181 |
weight | 1 | 0.0382 | 0.00948 | 16.2378 | <.0001 | |||
event | 1 | 1 | 1.5374 | 0.7670 | 4.0181 | 0.0450 | ||
event | 2 | 1 | 0.2637 | 1.0050 | 0.0689 | 0.7930 |
year 2012 vs 2014 | 3.421 | 2.615 | 4.475 |
---|---|---|---|
year 2013 vs 2014 | 1.831 | 1.375 | 2.439 |
Concordancia de porcentaje | 68.3 | D de Somers | 0.372 |
---|---|---|---|
Discordancia de porcentaje | 31.1 | Gamma | 0.374 |
Porcentaje ligado | 0.6 | Tau-a | 0.092 |
Pares | 1273250 | c | 0.686 |
On the other hand,
Would it be correct also use this other model?
proc genmod data=Vips descending;
class year treat sex;
model presence = treat Sex year....... / dist = binomial link = logit;
run;
Thank you very much,
Jaime
GENMOD would be a good choice, in my opinion. Make sure that all of the interactions are included, as it appears that many of them are significant, such that odds ratios are not constant at all levels of various CLASS variables. It almost certainly explains why the raw data plots have differing signs from the solution vector.
A key critical element is that there is an interaction between the continuous covariate 'weight' and the treatment by sex levels. This means that separate slopes need to be fit, and that the test for treatment by sex is a test of the intercept values in the current approach. Eventual comparison is going to involve use of the AT= option in either the LSMEANS or LSMESTIMATE statement.
Steve Denham
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.