Avian recruitment analysis

[SpanishPhD]

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

[SteveDenham]

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

[SpanishPhD]

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

[SteveDenham · Accepted solution]

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

[SpanishPhD]

Thank you very much for your useful assistance!!

All the best,

Jaime

[SpanishPhD]

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

[SteveDenham]

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

[SpanishPhD]

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

[SteveDenham]

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

[SpanishPhD]

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

[SpanishPhD]

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

[SteveDenham]

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