Set-up (1) Reading data

data <- read.table("linreg.dat", header=TRUE)
data
##     x    y
## 1  -1  1.4
## 2   0  4.7
## 3   1  5.1
## 4   2  8.3
## 5   3  9.0
## 6   4 14.5
## 7   5 14.0
## 8   6 13.4
## 9   7 19.2
## 10  8 18.0

Set-up (2) Writing cpp file as “linreg.cpp”"

#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(x);
  DATA_VECTOR(y);
  int n = y.size();

  PARAMETER(b0);
  PARAMETER(b1);
  PARAMETER(logSigma);
  vector<Type> yfit(n);

  Type neglogL = 0.0;

  yfit = b0 + b1*x;
  neglogL = -sum(dnorm(y, yfit, exp(logSigma), true));

  return neglogL;
}

Minimizing an objective function defined

model <- MakeADFun(data, parameters, DLL="linreg")
fit <- nlminb(model$par, model$fn, model$gr)
## outer mgc:  1468.4 
## outer mgc:  162.4893 
## outer mgc:  91.03071 
## outer mgc:  98.22766 
## outer mgc:  2.172514 
## outer mgc:  5.116226 
## outer mgc:  6.162891 
## outer mgc:  5.261146 
## outer mgc:  8.460502 
## outer mgc:  3.897932 
## outer mgc:  3.234142 
## outer mgc:  0.9433855 
## outer mgc:  0.3326979 
## outer mgc:  0.01294656 
## outer mgc:  0.009230682 
## outer mgc:  0.001157467 
## outer mgc:  0.0001261349 
## outer mgc:  1.271994e-05

Extracting estimates and their SEs and showing results

best <- model$env$last.par.best
rep <- sdreport(model)
## outer mgc:  1.271994e-05 
## outer mgc:  0.017538 
## outer mgc:  0.01756344 
## outer mgc:  0.1027843 
## outer mgc:  0.1028098 
## outer mgc:  0.01999046 
## outer mgc:  0.02000952
print(best)
##        b0        b1  logSigma 
## 4.0781824 1.9090907 0.3451266
print(rep)
## sdreport(.) result
##           Estimate Std. Error
## b0       4.0781824  0.7039414
## b1       1.9090907  0.1554747
## logSigma 0.3451266  0.2236068
## Maximum gradient component: 1.271994e-05