What model to use for analyzing age frequency data? issues with linear model in R

There are standard methods of analyzing fish-catch data in terms of the distribution of ages over calendar years, evaluating baseline mortality and the influence of factors like environmental conditions and depletion due to commercial fishing effort. This 2001 document* from the Food and Agriculture Organization of the United Nations (FAO) describes what's called "Virtual Population Analysis" in the fisheries world.

This cohort analysis models each year's population of fish of a particular age as a function of the population of fish one year younger from the previous year, thus going back to the calendar year of age-0 for each birth cohort. The models can be simple exponential decay models.

From that perspective, the linear model proposed in the question, modeling the number of age-0 fish as a function of numbers of each of the older-age fish to estimate some "correlation," has the causality backwards. You need to be modeling each year's catch of a certain age of fish as a function of the previous year's catch of 1-year younger fish. Also, as these are small count data, an ordinary least square analysis with lm() isn't appropriate.

A simple way to proceed is Poisson regression of counts versus age, with its log link to represent the usually assumed exponential decay of numbers of fish with age. The fish counts are modeled as a function of fish age (numeric), with birth cohort and catch year as factors. This way all of the data on a birth cohort are used, rather than relying on differences from the age-0 count alone, and correlations from year to year within each cohort are accounted for.

Re-format the data into long form, with one row for each count value and columns for the count, age, birth cohort, and year of catch. You can then examine all the data simply with xyplots from the lattice package:

xyplot(count~age|birthCohort,data=longDF)

which shows that there is some useful information available on birth cohorts going back to 2006 even though collection data don't start until 2008. In general, counts of each birth cohort decrease with age, documenting the expected "temporal correlations."

To account for differences in sizes of birth cohorts and the catch effort among calendar years, include those as random effects. As there seem to have been differences in efforts made to examine different ages in different catch years, allow the apparent influence of age on numbers to vary among catch years, represented by a random slope for age within catch year. Restrict to birth cohorts after 2005, based on the above visual inspection of the data.

The function call (using the lme4 package in R) is:

glmer(count~ age + (1|birthCohort) + (age|CatchYr),data=longDF,subset=as.numeric(as.character(birthCohort))>2005,family=poisson)

The log link, modeling the exponential decay in the coefficient for age, is the default for the Poisson family when called this way. A quick check with the DHARMa package suggests that the Poisson fit is fairly good, supporting this exponential decay model (not shown). The age coefficient, from the summary() of the model, is then the time constant for the exponential decay:

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.87544    0.23658   16.38  < 2e-16 ***
age         -0.54679    0.09781   -5.59 2.27e-08 ***

The intercept is the log of the estimated overall number of age-0 counts per birth cohort; the variability among birth cohorts is captured by the fairly large variance of its random effect, 0.48 (standard deviation, 0.69). The random slope and intercept associated with catch year both seem to be important (not shown).

So yes, there is a strong relationship between catch number and fish age within each birth cohort: an exponential decay with time constant approximately -0.55 per year of age.

---

From the initial version of this question, it seemed that the problem was missing data. After some back-and-forth among several individuals, it's now clear that there are no missing data, but rather true 0 catch values of fish at high ages overall, with perhaps some additional data problems at early years of the study. Much of the original answer has thus now been deleted. Please look at the editing history of both the question and this answer if you wish to make sense of some of the comments here.

*Lassen, H and Medley, P. Virtual Population Assessment--a practical manual for stock assessment (2001). FAO Fisheries Technical Paper 400.

Some notes on how this document applies here:

Much of the document is about using proxies of fish length, weight, etc for the actual age. As the current data evidently have correct ages based on otolith analysis those approximations and estimates (and associated efforts to estimate age distributions from bulk-scale catches) aren't needed here.

On the other hand, much of what's in the document is based on large scale data for which linear regression models would be expected to work well. The current data aren't, they are small-count data for which ordinary least squares analysis with lm() isn't appropriate. Count-based analysis with Poisson or related generalized linear models, noted but not emphasized in the document, is needed here.

---

Note on data re-formatting. It's important to develop some facility in moving from wide-format data as in your age.matrix (often an easy form for data capture from spreadsheets) to the long format that often is more useful for regression analysis. As an example, here's how I did it in this case, starting with a data-frame version of your age.matrix that I called age.df.

First, give more informative names to clarify the distinction between the year of the catch and the ages of the fish caught in each year, in a way that will simplify determining the birth year of the fish of a certain age caught in any year:

names(age.df) <- c("CatchYr", paste("age",as.character(0:9),sep="."))

Although there may be more intuitive functions for reshaping, I just used the standard reshape() function in R to create a long-form data frame with one count per row and associated annotations, longDF.

longDF <- reshape(age.df,direction="long",idvar="CatchYr",v.names="count",timevar="age",varying = paste("age",as.character(0:9),sep="."),sep=".")

The syntax for that function is tricky, and (as always) it took me a couple of tries to get it right. The direction specifies the direction of the output data frame. The idvar says which column of the starting data frame to use as the basis for identifying rows in the output, here CatchYr. v.names is what to call the column with the single values that are pulled out from the wide-form input into separate rows, in this case the "count" of fish of each age for each CatchYr. varying specifies the names of the columns in the wide format that will be parsed into corresponding identifiers in the long format, here the various age columns in the wide-form age.df. The results of that parsing are placed in an output column with name specified by timevar, here "age".

When I inspected the data frame I saw that the output "age" values ran from 1 to 10, so I subtracted 1 from all values to put them into the desired 0 to 9 range.

> longDF[,"age"] <- longDF[,"age"] - 1

Then I set up a new column to represent the birth cohort, the year in which the fish of a specific age caught in a specific year would have been born. That just required subtracting the age values from the CatchYr values, taking care about whether values were to this point specified as numeric or character variables.

> longDF[,"birthCohort"] <- as.character(as.numeric(longDF[,"CatchYr"])-longDF[,"age"])

Then I transformed the CatchYr and birthCohort values (currently character variables) into factors:

> longDF$CatchYr <- factor(longDF$CatchYr)
> longDF$birthCohort <- factor(longDF$birthCohort)

The summary of the resulting data frame:

> summary(longDF)
    CatchYr        age          count         birthCohort
 2008   :10   Min.   :0.0   Min.   :  0.00   2008   :10  
 2009   :10   1st Qu.:2.0   1st Qu.:  1.00   2009   :10  
 2010   :10   Median :4.5   Median :  3.00   2007   : 9  
 2011   :10   Mean   :4.5   Mean   : 14.67   2010   : 9  
 2012   :10   3rd Qu.:7.0   3rd Qu.: 16.00   2006   : 8  
 2013   :10   Max.   :9.0   Max.   :170.00   2011   : 8  
 (Other):50                                  (Other):56  

That was the data frame used for the mixed model.