topic Re: Avian recruitment analysis in Statistical Procedures

Avian recruitment analysis

SpanishPhD — Fri, 16 Jan 2015 16:02:29 GMT

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

Re: Avian recruitment analysis

SteveDenham — Fri, 16 Jan 2015 16:07:33 GMT

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

Re: Avian recruitment analysis

SpanishPhD — Fri, 16 Jan 2015 16:32:06 GMT

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

Re: Avian recruitment analysis

SteveDenham — Fri, 16 Jan 2015 18:32:00 GMT

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;

run;

Steve Denham

Re: Avian recruitment analysis

SpanishPhD — Mon, 19 Jan 2015 10:28:58 GMT

Thank you very much for your useful assistance!!

All the best,

Jaime

Re: Avian recruitment analysis

SpanishPhD — Mon, 19 Jan 2015 16:49:59 GMT

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

Re: Avian recruitment analysis

SteveDenham — Mon, 19 Jan 2015 20:11:28 GMT

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

Re: Avian recruitment analysis

SpanishPhD — Tue, 20 Jan 2015 12:26:34 GMT

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

Re: Avian recruitment analysis

SteveDenham — Tue, 20 Jan 2015 19:02:48 GMT

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

Re: Avian recruitment analysis

SpanishPhD — Wed, 21 Jan 2015 16:03:35 GMT

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
weighttreatevent	2	16.6071	0.0002
Sex*event	2	1.4899	0.4748
TreatSexweight	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
weighttreatevent	0	1		1	-0.00515	0.0116	0.1972	0.6570
weighttreatevent	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
Trt1Sexevent	0	1	1	1	-0.0274	0.0746	0.1350	0.7133
Trt1Sexevent	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

Re: Avian recruitment analysis

SpanishPhD — Wed, 21 Jan 2015 16:16:03 GMT

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

Re: Avian recruitment analysis

SteveDenham — Fri, 23 Jan 2015 15:24:59 GMT

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