# Simulate count data with underlying dynamics

I'm trying to simulating time-series of count data given an underlying population model in order to get a better feel for the processes shaping the data. As I am modelling the actual time-series of counts following eq (1) and (2) (substitute m for t), the aim of the below simData() is to generate count data given the underlying model.

Here's what I need help with: As I'm modelling the actual count data as a poisson distribution, I use a log link function to "link" the data and process model. However, it's not clear to me how to make this "link" when simulating the data. Below I'm using an inverse link transformation i.e. lambda[t+1,i] <- exp( mu[t+1,i] ), but this makes the simulated data experience exponential growth. If I remove the exponential, the simulated data 'looks' good. But at the same time, haven't I removed the necessary link function?

Model description: If we denote n[i,t] as the natural logarithm of N[i,t], then on a natural logarithmic scale we have the number of species i in time t described by where r[i] and k[i] represent the intrinsic growth rate and the carrying capacity of species i, respectively. α[i,j] represent the interaction coefficient between species i and j and expresses the per capita effect of species j on the growth rate of species j from time t−1 to time t. Finally, ɛ[i,t] represents the effect of unexplained, latent, stochastic noise on the population dynamics of species i. The time-series is modelled using a Poisson distribution, where y[i,t] denotes the count of species i in time t In simData, I have deliberately excluded the offset (i.e. log(N_t) + \Lambda_(m)).

simData <- function(r, k, alpha, TIME, NSpecies) {

#process error
sd <- runif(n = NSpecies, min = 0, max = 1000)^-2
# interaction matrix
alpha <- matrix(ncol=NSpecies,nrow=NSpecies, data=alpha, byrow = T)

#Create empty matrices
N <- matrix(numeric(length = (TIME+1)*NSpecies),ncol=NSpecies,nrow=TIME+1)
lambda <- matrix(numeric(length = (TIME+1)*NSpecies),ncol=NSpecies,nrow=TIME+1)
mu <- matrix(numeric(length = (TIME+1)*NSpecies),ncol=NSpecies,nrow=TIME+1)
mu.star <- matrix(numeric(length = (TIME+1)*NSpecies),ncol=NSpecies,nrow=TIME+1)

#Set initial values
N[1,] <- k
lambda[1,] <- k
mu[1,] <- c( 0.1, 0.1, 0.1 )
mu.star[1,] <- c( 0.1, 0.1, 0.1 )

for( t in 1:TIME ) {
for( i in 1:NSpecies ) {
mu.star[t+1,i] <- mu[t,i] + r[i] * ( 1 - (alpha[i,] %*% mu[t,]) / k[i] )
mu[t+1,i] <- rnorm( n = 1, mean = mu.star[t+1,i], sd = sd[i] )
lambda[t+1,i] <- exp( mu[t+1,i] )
N[t+1,i] <- rpois( n = 1, lambda = lambda[t+1,i] )
}
}

out <- list( N, lambda, mu, mu.star )
return(out)
}

# same r and k, and intra-specific across species. No inter-specific interactions
simData_out <- simData(r = c(0.1,0.1,0.1), k = c(100,100,100), alpha = c(1,0,0,0,1,0,0,0,1), TIME = 100, NSpecies = 3)

• Counts can be 0. How can you take the log of 0? – Glen_b Dec 29 '16 at 9:18
• @Glen_b Thanks for your comment. Indeed, so where would you change the code? Remove the exp()? Though, I reckoned it needed an inverse link transformation. – jO. Dec 29 '16 at 16:51
• Note that in general code review is off topic here (see the help center). If you can rephrase to ask about a primarily statistical issue it may be on topic – Glen_b Dec 31 '16 at 1:50
• @Glen_b Ok, sorry about that. I have now tried to update the question. – jO. Jan 2 '17 at 16:56