Correct coding of glimmix for both g and r-side random effects

[BirdMan]

Hi there,

I hope someone might be able to help me out with some coding for proc glimmix. I'm not as familiar with coding r-side random effects and have been banging my head against this problem for too long so I hope that this forum might be able to offer some suggestions. I'll try to give a detailed description of the study design first. This is a study dealing with bird counts in Canada. Our first step was to randomly select 9 41-km square "landscapes" that fell into three "treatments" (the landscapes were dominated by different vegetation types) and each treatment was replicated three times (so there were three 41-km square landscapes that were treatment 1, three landscapes that were treatment 2, and three that were treatment 3), for a total of 9 different 41-km square landscapes. Within these landscapes we surveyed for birds in accessible parcels of land using a point count method (i.e., stand at one spot and listen for birds). This study was done over 2 years and each point count location was surveyed once per year. There were approximately 170 point count locations. Many of our species were either present or absent (binary dependent variable) at each point count and we are interested in exmaining how bird presence or absence varies in response to various vegetation parameters (vegetation height, etc.), as well as some other parameters of interest (landscape metrics, etc.). But given the repeated measures I am having some difficulty coding the glimmix script. This is what I've come up with so far and I'd be interested in getting some feedback. I hope my descriptions of the study design helps.

proc glimmix data=XX ;

class Year Land_Type LandID;

model occ(event='1')=vht /dist=binary link=logit solution;

random landid(Land_type);

random Year/Subject=PCSite type=ar(1) residual;

run;

A collaborator suggested that we incldue the individual "landscape" nested within treatment effect (the landid is simply an indentification number given to each individual landscape irrespective of treatment) and land_type is the treatment type (3 possible treatments), but I am stuck on how to code the repeated measures of the same site between years. Any suggestions would be extremely helpful! Please if I did not present enough information let me know.

Thanks for any help/suggestions!

Ryan

[SteveDenham]

If I understand everything correctly, you have two observations for each PCSite--one in year1 and one in year2.  If this is the case, you might as well specify type=chol (an unstructured Cholesky parameterization) for the covariance structure.  There are only three R-side parameters to be estimated. 

So I would only slightly change the R side specification to:

random year/type=chol subject=pcsite residual;

and make sure that pcsite is included in the class statement.

Now back to the G side. I get the impression that landid is unique, and so a given sequence of landid's is nested within each landscape, as some sort of sub-sample.  How does landid differ from pcsite?  This will motivate what I say next.  Are the point count sites clustered within the landid's?  If so, then maybe something like:

random intercept land_type/subject=landid;

would be appropriate.  This should give two variance components: one arising from landid's within land_types and one arising from between land_types.

If the pcsites are NOT clustered within the landid's, then we need to consider just what the experimental unit is, and maybe start to think about spatial location covariance structures.

Steve Denham

[BirdMan]

Hi Steve,

Thanks so much for your quick and thoughtful response. Hopefully I can clear up your second question related to the g-side random effects. You are correct in that the point count sites are clustered within the landid's and that landid is unique.

So would I remove the: random landid(Land_type) statement and replace with your suggestion of random intercept land_type/subject=landid (or would your suggestion get added to my existing code?).

Thanks so much for the help!

Ryan

[SteveDenham · Accepted solution]

Replace.  Otherwise, it will try to estimate the between land_types twice, and set one to zero, or some other odd thing.

Final suggested code:

proc glimmix data=XX ;

class Year Land_Type LandID PCSite;

model occ(event='1')=vht /dist=binary link=logit solution;

random intercept Land_type/subject=LandID;

random Year/Subject=PCSite type=chol residual;

run;

Hope this works...

Steve Denham