\documentclass{article}
\usepackage[latin1]{inputenc}
%\usepackage[brazil,brazilian]{babel}
\title{Towards a joint model for abundance and age compositions}
\author{Ernesto Jardim \& Paulo Justiniano Ribeiro Jr}
\date{\today}
\begin{document}
\maketitle
\section{Introduction}
The model for analysis proposed in Jardim \& Ribeiro in the paper submitted to CJFAS
relies on a fundamental assumption of independence between the total abundance (henceforth $I$)
and the proportions (henceforth $P$) of age groups.
We argue this is not a fundamental issue for the analysis and the method proposed in the paper
remains valid for the target of the analysis within reasonable deviations of such assumption.
However, this has risen questions by the referes and a simulation study is suggested in order to
further investigate such assertion.
In this document we discuss an way to generate simulations for total abundance and
age compositions in order to run a simulation study to assess the effect of the interaction between
I and P.
We seek for mechanism for generating different scenarios ranging from independence,
weak dependence and strong dependence between $I$ and $P$.
Although the focus here is in how to generate the simulations the ideas could be
used as an attempt to build a joint model.
We propose a mechanism (model) for generating age compositions and abundance at age,
possibly related with the total abundance.
The intuition for the model is as follows: the basic ``ingredient'' is a Gaussian random field
which describes the spatial variation determined by the environmental conditions.
The total abundance is then related with such field.
Considering the proportions for each age group could be dependent of
total abundance, each one of them is generated controled
by three parameters: one related with the overall level of the proportion,
another controlling the strength of the relation with the total abundance and
a third allowing for variation not related to the total abundance.
An alternative model could be proposed by generating the abundance for each age class
and them adding them to produce the total abundance.
Although apparently more appealing this takes a fundamental and different view.
This model would not consider (at least in its basic format)
the fact the age groups could be in some sort of competition and therefore interacting.
Ideias for such models are given, for instance, by "common component models" where the
abundances at age are genetrated including a common spatial component.
Alternative models options such as this one could be proposed and investigated in other pieces of work.
Comparisons and intuition for the models could be further explored.
\section{The proposed model I (mechanism)}
Denote the spatial locations by $x$.
Define $R(x)$ a standardised Gaussian random field, i.e. with zero mean, unit variance and
correlation function $\rho(u; \phi)$ where $u$ denotes the distance between a pair of spatial
locations and $\phi$ is a parameter controlling the rate of decay of the correlation as function
of the distance. In particular we assume $\rho(u) = \exp(-u/\phi)$.
The total abundance is modelled as log-normal Gaussian random field
given by $Y(x) = \exp\{\mu + \sigma R(x) + \epsilon \}$
where $\epsilon$ is a vector of
i.i.d. Gaussian random variables with variance $\tau^2$.
Consider now $k = 1, 2, \dots, K$ denotes the age classes.
For each age class consider $S(x) = \beta_{0k} + \beta_{1k} R(x) + \epsilon_k$
where $\epsilon_k$ is a vector of
i.i.d. Gaussian random variables with variance $\tau_k^2$.
Compute $q_k(x) = \exp\{S_k(x)\}/(1+\exp\{S_k(x)\})$
and the proportions at each location by $p_k(x) = q_k(x)/\sum_k(q_k(x))$
The abundance at age are given simply by $Y_k(x) = p_k(x) \dot Y(x)$.
The proposed model has $3 K + 4$ parameters. A more parsimonious version could consider
a common parameter $\tau_k \forall k$ ($5 + 2 K$ parameters)
or even $\tau_k = tau \forall k$ ($4 + 2 K$ parameters).
A function simulating for such model is defined as follows.
<>=
gcam1 <- function(x, y.pars, pk.pars){
options(geoR.messages=F)
require(geoR)
# x is a n x 2 matrix with coordinates for the locations
# y.pars is a vector (mu, sigma, phi, tau)
# pk.pars is a K x 3 matrix with each row given by (beta_0k, beta_1k, tau_k)
n <- nrow(x)
K <- nrow(pk.pars)
R <- grf(grid=x, cov.pars=c(1, y.pars[3]))$data
Y <- exp(y.pars[1] + y.pars[2] * R + rnorm(n, sd=y.pars[4]))
S <- t(pk.pars[,1] + kronecker(pk.pars[,2], t(R)) + rnorm(K*n, sd=pk.pars[,3]))
Q <- exp(S)/(1+exp(S))
P <- Q/rowSums(Q)
colnames(P) <- paste("p", 1:K, sep="")
return(cbind(Y,P))
}
@
Consider now a first simulation.
<<>>=
set.seed(123)
locs <- expand.grid((0:20)/20, (0:20)/20)
Ypars = c(0.5, sqrt(0.9), 0.2, sqrt(0.1))
Ppars <- rbind(c(0, sqrt(0.45), sqrt(0.05)), c(0, sqrt(0.25), sqrt(0.25)), c(0, -sqrt(0.25), sqrt(0.25)))
sim1 <- gcam1(locs, Ypars, Ppars)
apply(sim1,2,range)
plot(as.data.frame(sim1))
cor(sim1)
cor(sim1, meth="sp")
@
<>=
plot(as.data.frame(sim1))
@
And a second one with different parameters.
<<>>=
set.seed(123)
Ppars <- rbind(c(10, sqrt(0.45), sqrt(0.05)), c(-3, sqrt(0.25), sqrt(0.25)), c(0, -sqrt(0.25), sqrt(0.25)))
sim2 <- gcam1(locs, Ypars, Ppars)
apply(sim2,2,range)
plot(as.data.frame(sim2))
cor(sim2)
cor(sim2, meth="sp")
@
<>=
plot(as.data.frame(sim2))
@
\section{Parameter effects and interpretation}
Innitially we considered the parameter effects were as follows:
\begin{itemize}
\item $\beta_{0k}$ controls to the overall level of the proportion for each class
\item $\beta_{1k}$ shape and strength of the association with $Y$
\item $\beta_{1k}/\tau_k$ strength of the association with $Y$
\item $\beta_{1k} + \tau_k$ variance of the proportions
\end{itemize}
however, is not that simple, probably because the non-linearity in the model.
The last one is not true, see the difference:
<<>>=
Ppars[,2]^2+Ppars[,3]^2
apply(sim1,2,var)
@
In fact the signal of $\beta_1$ is also relevant. Besides that, because of the
operation forcing the proportions add to one makes
causes an interaction.
Inspect this noticing the age with signal for $\beta_1$ different from the other two have greater variability:
<<>>=
Ppars <- rbind(c(0, sqrt(0.25), sqrt(0.05)), c(0, sqrt(0.25), sqrt(0.25)), c(0, sqrt(0.25), sqrt(0.25)))
set.seed(1211); apply(gcam1(locs, Ypars, Ppars), 2, range)
Ppars <- rbind(c(0, sqrt(0.45), sqrt(0.05)), c(0, sqrt(0.25), sqrt(0.25)), c(0, sqrt(0.05), sqrt(0.45)))
set.seed(1211); apply(gcam1(locs, Ypars, Ppars), 2, range)
Ppars <- rbind(c(0, sqrt(0.45), sqrt(0.05)), c(0, sqrt(0.25), sqrt(0.25)), c(0, -sqrt(0.25), sqrt(0.25)))
set.seed(1211); apply(gcam1(locs, Ypars, Ppars), 2, range)
Ppars <- rbind(c(0, sqrt(0.45), sqrt(0.05)), c(0, -sqrt(0.25), sqrt(0.25)), c(0, -sqrt(0.25), sqrt(0.25)))
set.seed(1211); apply(gcam1(locs, Ypars, Ppars), 2, range)
Ppars <- rbind(c(0, -sqrt(0.45), sqrt(0.05)), c(0, sqrt(0.25), sqrt(0.25)), c(0, -sqrt(0.45), sqrt(0.05)))
set.seed(1211); apply(gcam1(locs, Ypars, Ppars), 2, range)
@
Also see that $\bar p_k$ change with variance.
<<>>=
set.seed(111)
apply(gcam1(locs, Ypars, rbind(c(1.2, 0.5, 0.05), c(0.1, 0.1, 0.2), c(0.7, -0.2, 0.2))),2,mean)
set.seed(111)
apply(gcam1(locs, Ypars, rbind(c(1.2, 0.05, 0.5), c(0.1, 0.1, 0.2), c(0.7, -0.2, 0.2))),2,mean)
@
but there is no change here. Probably this was before correcting the code.
Also see that $var p_k$ changes even if $\beta_{1i} + \tau_k$ is constant.
<<>>=
set.seed(111)
apply(gcam1(locs, Ypars, rbind(c(1.2, 0.5, 0.05), c(0.1, 0.1, 0.2), c(0.7, -0.2, 0.2))),2,var)
set.seed(111)
apply(gcam1(locs, Ypars, rbind(c(1.2, 0.05, 0.5), c(0.1, 0.1, 0.2), c(0.7, -0.2, 0.2))),2,var)
@
Maybe these are just small issues but I can not fully understand how these parameters interact. Note that trying to induce the same correlation in opposite directions does not work. The example below I'd expect to show the same correlation between $k=1$ and $k=3$ but with opposite signals.
<<>>=
Ppars <- rbind(c(1.2, 0.2, 0.2), c(0.1, 0.1, 0.2), c(1.2, -0.2, 0.2))
set.seed(111)
sim <- gcam1(locs, Ypars, Ppars)
cor(sim)
@
Note that I can not increase correlation in $p_1$ and although I'm simulating independence in $p_2$ and $p_3$ it does not occurr.
<<>>=
Ppars <- rbind(c(1.2, 10, 0.05), c(0.1, 0, 0.2), c(0.7, 0, 0.2))
set.seed(111)
sim <- gcam1(locs, Ypars, Ppars)
cor(sim)
@
\section{Scenarios for the simulation study}
For the simuation study to be included in the paper we consider the following
grid of points and parameters for the total abundance.
<<>>=
gs <- expand.grid((0:40)/40, (0:40)/40)
Ypars = c(0.5, sqrt(0.9), 0.2, sqrt(0.1))
@
We consider three scenarios for simulation reflecting different levels of association
between the proportions at age and totasl abundance
\begin{emumerate}
\item strong association
<<>>=
Ppars <- rbind(c(10, sqrt(0.45), sqrt(0.05)), c(-3, sqrt(0.25),
sqrt(0.25)), c(0, -sqrt(0.25), sqrt(0.25)))
set.seed(123)
sim.strong <- gcam1(locs, Ypars, Ppars)
cor(sim.strong, met="sp")
plot(as.data.frame(sim.strong))
@
\item week association
<<>>=
Ppars <- rbind(c(10, sqrt(0.01), sqrt(0.4)), c(-3, sqrt(0.01),
sqrt(0.4)), c(0, -sqrt(0.01), sqrt(0.4)))
set.seed(123)
sim.week <- gcam1(locs, Ypars, Ppars)
cor(sim.week, met="sp")
plot(as.data.frame(sim.week))
@
\item no association (complete independence)
<<>>=
Ppars <- rbind(c(0, 0, sqrt(0.5)), c(0, 0, sqrt(0.5)), c(0, 0, sqrt(0.5)))
set.seed(123)
sim.indep <- gcam1(locs, Ypars, Ppars)
cor(sim.indep, met="sp")
plot(as.data.frame(sim.indep))
@
\section{The proposed model II (mechanism)}
This model has the same ingredients as the previous one.
However, instead of computing probabilities as functions of the Gaussian random field.
So the steps are to simulate $R(x)$ and $Y(x)$ as for model~I and them compute
the abundance at age by
$S_k(x) = exp(\beta_{0k} + \beta_{1k } R(x) + e_k)$ for $k = 1, \dots K-1$
$S_K(x) = Y(x) - \sum_k^{K-1} S_k(x)$.
Then the age proportions as simply given by
$p_k(x) = S_k(x)/Y(x)$.
<>=
gcam2 <- function(x, y.pars, pk.pars){
options(geoR.messages=F)
require(geoR)
# x is a n x 2 matrix with coordinates for the locations
# y.pars is a vector (mu, sigma, phi, tau)
# pk.pars is a (K-1) x 3 matrix with each row given by (beta_0k, beta_1k, tau_k)
n <- nrow(x)
K <- nrow(pk.pars)
R <- grf(grid=x, cov.pars=c(1, y.pars[3]^2))$data
Y <- exp(y.pars[1] + y.pars[2] * R + rnorm(n, sd=y.pars[4]))
S <- t(pk.pars[,1] + kronecker(pk.pars[,2], t(R)) + rnorm(K*n, sd=pk.pars[,3]))
S <- cbind(S, Y-apply(S,1,sum))
P <- S/Y
colnames(P) <- paste("S", 1:K, sep="")
return(cbind(Y,P))
}
@
Alternatively we can start by
simulating $R(x)$ followed by computing all $S_k(x)$, k=1, \dots, K$.
From these, the total abundance is $Y(x) = \sum_k S_k(x)$ .
<>=
gcam3 <- function(x, phi, pk.pars){
options(geoR.messages=F)
require(geoR)
# x is a n x 2 matrix with coordinates for the locations
# phi is the parameter for the exponential correlation function
# pk.pars is a K x 3 matrix with each row given by (beta_0k, beta_1k, tau_k)
n <- nrow(x)
K <- nrow(pk.pars)
R <- grf(grid=x, cov.pars=c(1, phi))$data
S <- t(pk.pars[,1] + kronecker(pk.pars[,2], t(R)) + rnorm(K*n, sd=pk.pars[,3]))
Y <- rowSums(S)
P <- S/Y
colnames(P) <- paste("p", 1:K, sep="")
return(cbind(Y,P))
}
@
For both case the model has $3K + 1 $ parameters and therefore more parcimonius them model I.
\end{document}