This file shows how to convert diet.tpl into diet_sk.tpl, i.e. how go from normal random effects (re's). For normally distributed re's the log-likelihood is // Exposure model g -= -0.9189385 - 0.5*square(u_j); // prior on the u_j's dvariable eta_j = X1(j)*gamma + mfexp(0.5*log_psi)*u_j; which is converted to the skewed re model by transforming the u_j to the skewed zeta_j and using zeta_j as the re in the code. // Exposure model g -= -0.9189385 - 0.5*square(u_j); // normal prior on the u_j's dvariable a2 = square(a)+square(1-a); dvariable zeta_j = 1.0/sqrt(a2)*(a*u_j + (1-a)*(mfexp(u_j)-exp(.5))/2.161197); dvariable eta_j = X1(j)*gamma + mfexp(0.5*log_psi)*zeta_j; Ths involves the extra parameter a which must be added into the function call so that SEPARABLE_FUNCTION void loglik_j(const dvar_vector& gamma, const dvariable& log_psi, const dvar_vector& alpha, const dvariable& log_theta, const dvariable& lambda, const dvar_vector& beta, const dvariable& u_j, const int& j) is changed to SEPARABLE_FUNCTION void loglik_j(const dvar_vector& gamma, const dvariable& log_psi, const dvar_vector& alpha, const dvariable& log_theta, const dvariable& lambda, const dvar_vector& beta, const dvariable& u_j, const dvariable& a, const int& j) Then the parameter a is added to the PARAMETER_SECTION so that PARAMETER_SECTION init_bounded_vector gamma(1,3,-10.,10.,2) init_bounded_number log_psi(-1.,4.,2) init_bounded_vector alpha(0,1,0.,30.,1) init_bounded_number log_theta(-1.,4.,1) init_bounded_number lambda(-6.,5.,2) init_bounded_vector beta(0,3,-10.,10.,1) random_effects_vector u(1,n,2) // N(0,1) distributed r objective_function_value g becomes PARAMETER_SECTION init_bounded_vector gamma(1,3,-3.0,3.0,2) init_bounded_number log_psi(0.0,4.0,2) init_bounded_vector alpha(0,1,0.0,30.0,1) init_bounded_number log_theta(0.0,4.0,1) init_bounded_number lambda(-0.3,0.3,2) init_bounded_vector beta(0,3,-3.0,3.0,1) init_bounded_number a(0.3,0.999,3) //add the parameter a random_effects_vector u(1,n,2) // N(0,1) distributed random effects objective_function_value g // (negative) Log-likelihood Finally the parameter a is aded to the function call so that for(j=1;j<=n;j++) loglik_j(gamma,log_psi,alpha,log_theta,lambda,beta,u(j),j); becomes for(j=1;j<=n;j++) loglik_j(gamma,log_psi,alpha,log_theta,lambda,beta,u(j),a,j); and that's it!