Programming the statistical procedures from SAS

Avian recruitment analysis

Accepted Solution Solved
Reply
Occasional Contributor
Posts: 18
Accepted Solution

Avian recruitment analysis

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


Accepted Solutions
Solution
‎01-16-2015 01:32 PM
Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

View solution in original post


All Replies
Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

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:

BirdTreatmentSexPresence_2011Presence_2012Presence_2013Presence_Globalweight
35370102000056
35370202111139
35370302100154
35370402001163
35370511001160
35370611110145
35370711010152
35370811000049
35370911001138

Thank you very much

All the best;

Jaime Muriel

Solution
‎01-16-2015 01:32 PM
Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

Thank you very much for your useful assistance!!

All the best,

Jaime

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

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

Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

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

Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

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 datosWORK.VIPS
Variable de respuestapresence
Número de niveles de respuesta2
Modelologit binario
Técnica de optimizaciónPuntuación de Fisher

   
Número de observaciones leí3330
Número de observaciones usa3213

    
11463
202750


La probabilidad modelada es presence =1.

Note:117 observations were deleted due to missing values for the response or explanatory variables.

     
Treat01
1-1
Sex11
2-1
event110
201
3-1-1
year201210
201301
2014-1-1

  
Criterio de convergencia (GCONV=1E-8) satisfecho.

    
AIC2651.7012517.243
SC2657.7762638.742
-2 LOG L2649.7012477.243

     
Ratio de verosim172.458419<.0001
Puntuación160.019119<.0001
Wald148.230319<.0001

      
Treat113.50120.0002
Sex15.22510.0223
Treat*Sex10.30690.5796
year285.8229<.0001
weight*Treat113.62870.0002
Treat*event218.19900.0001
weight*event27.59960.0224
weight*treat*event216.60710.0002
Sex*event21.48990.4748
Treat*Sex*weight29.98840.0068
weight116.2378<.0001
event27.51740.0233

           
Intercept 1-4.33470.614149.8221<.0001
Treat0 1-2.25450.613613.50120.0002
Sex1 10.14300.06265.22510.0223
Treat*Sex01 1-0.03460.06250.30690.5796
year2012 10.61820.069978.1565<.0001
year2013 1-0.006550.07590.00750.9312
weight*Treat0 10.03500.0094813.62870.0002
treat*event01 10.15800.76700.04240.8368
treat*event02 1-3.52301.005212.28310.0005
weight*event1 1-0.02360.01164.14260.0418
weight*event2 1-0.003980.01560.06470.7993
weight*treat*event01 1-0.005150.01160.19720.6570
weight*treat*event02 10.05470.015612.23680.0005
Sex*event11 1-0.01940.07460.06780.7945
Sex*event12 10.11580.10271.27010.2597
Trt1*Sex*event0111-0.02740.07460.13500.7133
Trt1*Sex*event0121-0.24290.10275.58960.0181
weight 10.03820.0094816.2378<.0001
event1 11.53740.76704.01810.0450
event2 10.26371.00500.06890.7930

     
year 2012 vs 20143.4212.6154.475
year 2013 vs 20141.8311.3752.439

     
Concordancia de porcentaje68.3D de Somers0.372
Discordancia de porcentaje31.1Gamma0.374
Porcentaje ligado0.6Tau-a0.092
Pares1273250c0.686

Occasional Contributor
Posts: 18

Re: Avian recruitment analysis

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

Respected Advisor
Posts: 2,655

Re: Avian recruitment analysis

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

🔒 This topic is solved and locked.

Need further help from the community? Please ask a new question.

Discussion stats
  • 11 replies
  • 500 views
  • 0 likes
  • 2 in conversation