crayfish data analysis using TMB

gistfile1.r

# this is what i have so far

require(TMB)

rak <- read.table("./data/master_all_r3.txt", dec = ".", header = TRUE, sep = "\t")

# prepare the data
rakdat <- na.omit(rak[, c("vrsta", "tw", "ETS_t_1000")])
xvar <- rakdat[, "tw"]
yvar <- rakdat[, "ETS_t_1000"]
inter.new <- seq(from = min(xvar), to = max(xvar), by = 0.1)
extra.new <- seq(from = 0, to = 40, by = 0.1)

# construct model matrix
rakdat.mm <- model.matrix( ~ 0 + vrsta, data = rakdat)

# compile code
compile("./cpp_source//ricker_compare_lvls.cpp")
dyn.load("./cpp_source/ricker_compare_lvls")

# prepare parameters for the model
model.data <- list(x = xvar, y = yvar, inter_new = inter.new,
                   extra_new = extra.new, mm = rakdat.mm)
nspp <- ncol(rakdat.mm)
model.parameters <- list(a = rep(0.001, nspp), b = rep(-0.5, nspp), c = rep(100, nspp), logsigma = 3)

# full model
full.model <- MakeADFun(data = model.data, parameters = model.parameters,
                           DLL = "ricker_compare_lvls")

This is the C++ function that crashes R.

#include <TMB.hpp>                                // Links in the TMB libraries

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(x);                                // Data vector transmitted from R
  DATA_VECTOR(y);
  DATA_VECTOR(inter_new);
  DATA_MATRIX(mm);

  PARAMETER_VECTOR(a);                                  // Parameter value transmitted from R
  PARAMETER_VECTOR(b);
  PARAMETER_VECTOR(c);
  PARAMETER(logsigma);

  vector<Type> xmax = -(Type(1)/(mm*b).array());
  vector<Type> ymax = -((mm*a).array()/(mm*b).array()) * Type(exp(-1)) + (mm*c).array();

  vector<Type> intra_fit = (mm*a).array() * inter_new * exp((mm*b).array() * inter_new).array() + (mm*c).array();
  //extra_fit = (mm*a).array() * extra_new * exp((mm*b).array() * extra_new) + mm*c;

  ADREPORT(intra_fit);
  //ADREPORT(extra_fit);
  ADREPORT(xmax);
  ADREPORT(ymax);

  Type nll;                                                     
  Type sigma = exp(logsigma);

  nll = -sum(dnorm(y, (mm*a).array() * x * exp((mm*b).array() * x) + mm*c, sigma, true));    // Use R-style call to normal density

  return nll;
}