// Weibull regression with random effects in the BUGS manual: Example 3 in Skaug and Fournier (2003)

DATA_SECTION

  init_int n			// Number of inidviduals
  init_matrix age(1,n,1,2)	// Age
  init_ivector sex(1,n)		// Sex
  init_ivector disease(1,n)	// Disease type
  init_matrix t(1,n,1,2)	// Survival time
  matrix logt(1,n,1,2)		// Log-survival time
  init_matrix delta(1,n,1,2)	// Cencoring indicator

PARAMETER_SECTION

  // Parameters in related to covariates
  init_bounded_number b0(-8,2)			// Intercept
  init_bounded_number b_age(-.1,.1)		// Age effect
  init_bounded_vector b_disease(1,3,-3,3) 	// Contrast for codes 2,3,4
  init_bounded_number b_sex(-3,3)		// Sex effect

  init_bounded_number r(.5,2)			// Shape parameter in Weibull distribution
  init_bounded_number log_sigma(-6.0,1.0)	// 0.5*log-variance component

  random_effects_vector u(1,n)			// Random effect

  objective_function_value g			// Joint log-likelihood

PRELIMINARY_CALCS_SECTION
  cout << setprecision(4);
  logt = log(t);

PROCEDURE_SECTION

  int i, j;

  dvariable sigma, eta, h0t, ht, Ht;                      // Local variables 

  sigma = exp(log_sigma);

  g = -n*log(sigma)-0.5*norm2(u/sigma); // log[h(u)]

  for(i=1;i<=n;i++)                     
  {
    for(j=1;j<=2;j++)                     
    {
      eta = b0 + b_age*age(i,j) + b_sex*sex(i) + u(i);
      if(disease(i)>1)
        eta += b_disease(disease(i)-1);

      h0t = r*exp(logt(i,j)*(r-1));     // Baseline hazard
      ht = h0t * mfexp(eta);            // Hazard function    
      Ht = ht*t(i,j)/r;                 // Cumulative hazard
 
      g += delta(i,j)*log(ht) - Ht;     // log[f(y|u)] 
    }
  }

  g *= -1;

REPORT_SECTION

TOP_OF_MAIN_SECTION
  arrmblsize = 400000L;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(300000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(20000);
  gradient_structure::set_MAX_NVAR_OFFSET(2000);