www.phidot.org

I am working with nest survival data, collected from different plots at 2 different sites and in 2 different years.

I am interested in how survival differs between plot, site and year. The ‘plot’ variable is nested in site (and site nested in year) and I am wondering about the best way to deal with this using RMark.

I'd like to include plot, site and year as fixed effects in the RMark model. So far, I have a set of models, including Year+ Site + Plot, and also Year*Site*Plot etc.

My question: what can / should I do to deal with the nested fixed effects? Is it enough to make sure that year and site are included in the model, and also to include an interaction?

If I wasn't working with RMark, another option would be to try using 'aov' in R and to include an error term in the model, to allow for the plots within each site not being independent replicates (and also the sites within each year not being independent).

A random effect model wouldn't be useful for these data as I'm interested in testing all of the terms as fixed effects, and I only have 2 levels of both site and year (so not enough to include as random anyway).

Any help would be greatly appreciated.

Best wishes

El

It all depends on how you coded your plot factor variable. I'm assuming here that you used year,site and plot as group variables. If you coded plots with an identifier within site but they are the same across sites

For example

Site Plot
A 1
A 2
B 1
B 2

Then ~site+plot would fit 3 parameters: 1) an intercept for Site A, plot 1, 2) a site B value, and 3) a plot 2 value. This would treat plot 1 and 2 as the same across sites.

If you fitted ~site*plot it will fit 4 parameters: 1) an intercept for Site A, plot 1, 2) a site B value, and 3) a plot 2 value, 4) a site B-Plot 2 interaction value. This would treat plots differently across sites which is what you want because site is nested in plot.

Now if your data was setup as

Site Plot
A 1
A 2
B 3
B 4

Then ~plot will fit a model with 4 parameters and is the same as ~site*plot with plot defined as 1 or 2. If you used ~site*plot with plot defined as 1 to 4, it will be over-parameterized because it will try to fit 8 parameters. If instead you used ~-1+site:plot then it will fit 4 parameters.

Since year is crossed with site and plot then it can be treated as either additive or an interaction. The interaction with year could be over just site if the thinking was a year difference would differ across sites but the same within plot.

You can work out a lot of this for yourself by creating a dummy dataframe and using model.matrix to look at the design matrix. Even though RMark does the work of creating DMs for you, you need to understand DMs and what the formula is creating. Here are some examples of what I described above

Code: Select all

> df=data.frame(year=rep(c("1","2"),each=4),plot=rep(c("A","B"),each=2),site=rep(c("1","2"),4))
> df
  year plot site
1    1    A    1
2    1    A    2
3    1    B    1
4    1    B    2
5    2    A    1
6    2    A    2
7    2    B    1
8    2    B    2
> str(df)
'data.frame':   8 obs. of  3 variables:
 $ year: Factor w/ 2 levels "1","2": 1 1 1 1 2 2 2 2
 $ plot: Factor w/ 2 levels "A","B": 1 1 2 2 1 1 2 2
 $ site: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2

> model.matrix(~site+plot,df)
  (Intercept) site2 plotB
1           1     0     0
2           1     1     0
3           1     0     1
4           1     1     1
5           1     0     0
6           1     1     0
7           1     0     1
8           1     1     1
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$site
[1] "contr.treatment"

attr(,"contrasts")$plot
[1] "contr.treatment"

> model.matrix(~site*plot,df)
  (Intercept) site2 plotB site2:plotB
1           1     0     0           0
2           1     1     0           0
3           1     0     1           0
4           1     1     1           1
5           1     0     0           0
6           1     1     0           0
7           1     0     1           0
8           1     1     1           1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$site
[1] "contr.treatment"

attr(,"contrasts")$plot
[1] "contr.treatment"

> model.matrix(~year+site*plot,df)
  (Intercept) year2 site2 plotB site2:plotB
1           1     0     0     0           0
2           1     0     1     0           0
3           1     0     0     1           0
4           1     0     1     1           1
5           1     1     0     0           0
6           1     1     1     0           0
7           1     1     0     1           0
8           1     1     1     1           1
attr(,"assign")
[1] 0 1 2 3 4
attr(,"contrasts")
attr(,"contrasts")$year
[1] "contr.treatment"

attr(,"contrasts")$site
[1] "contr.treatment"

attr(,"contrasts")$plot
[1] "contr.treatment"

> model.matrix(~year*site+site*plot,df)
  (Intercept) year2 site2 plotB year2:site2 site2:plotB
1           1     0     0     0           0           0
2           1     0     1     0           0           0
3           1     0     0     1           0           0
4           1     0     1     1           0           1
5           1     1     0     0           0           0
6           1     1     1     0           1           0
7           1     1     0     1           0           0
8           1     1     1     1           1           1
attr(,"assign")
[1] 0 1 2 3 4 5
attr(,"contrasts")
attr(,"contrasts")$year
[1] "contr.treatment"

attr(,"contrasts")$site
[1] "contr.treatment"

attr(,"contrasts")$plot
[1] "contr.treatment"

>


Hi

Thank you very much for your response, it was really useful.

I wanted to try including an interaction between year and plot, as I think that survival in plots will behave differently between years.

I have coded the data as follows:

Site Plot
A 1
A 2
B 1
B 2



etc.

Then I’ve tried the model:

~year*site*plot

However, I’ve checked the design matrix and I think this may be including too many parameters. For example, the combination year*plot is included (for each plot), which will treat all Plot 1s as the same, whereas I want plots to be treated differently across sites.

Any further advice greatly appreciated.

You didn't really explain why you thought it was over-parameterized. I assume it was because some of the parameter had large Std errors or conf intervals looked strange, It isn't over-parameterized but it could be that for one of the combinations the parameter is at a boundary. As you explained it, you have 2 years, 2 sites and 2 plots so you have 8 combinations. At the most you should have 8 parameters to describe the mean values. You could have additional individual or group covariates if there is contrast within the 8 groups (eg sex).

Now your confusion may be due to interpretation, so what I suggest is that you fit ~-1+Year:Site:Plot which will fit 8 separate parameters and will identify any parameter that is at a boundary if that is the problem.

--jeff

When I first emailed, I simplified the design of the data slightly.

The data are actually as follows:

Year 1

Site 1, plot 1,2,3

Site 2, plot 1,2

Year 2

Site 1, plot 1,2,3,4,

Site 2, plot 1,2

This is the error message I am getting:

Note: only 11 parameters counted of 12 specified parameters

Looking at the standard error of the beta estimates gives a clue to what is going on (very high for S:fPlot4 and S:fYear2:fPlot4)

The problem seems to be because there is no plot 4 in year 1, Site 1. The RMark workshop notes suggest re-levelling the data, to change the intercept.

I changed the order of the data, so that Site 1 and year 2 were in the intercept, and the standard errors have now gone down, but I'm still getting the error message:

Note: only 11 parameters counted of 14 specified parameters

Also, I'm still not sure that it's right that e.g. in this case:

S:fYear1:fPlot2

All plot 2s are being treated the same, even though they are in different sites.

Thanks for any help, I just want to make sure that I am specifying the model properly.

Yes this is different because numbers of plots are unequal. What was described previously had the same number of plots per site and year. If you use ~-1+year:site:plot it will fit a separate estimate for each plot in each site in each year assuming each is a factor variable.

Hi Jeff

This has worked, thanks for all your help. Very much appreciated.

El