Revision 1016 trunk/examples/admb-re/glmmadmb/glmmadmb.tpl

glmmadmb.tpl (revision 1016)
1
// 2011-03-24 version from Dave Fournier
2
// modified 2011-04-27 Hans Skaug, BMB
3
DATA_SECTION
4

  
5
  init_int n					// Number of observations
6
  init_int p_y					// Dimension of y(i) (multivariate response)
7
  init_matrix y(1,n,1,p_y)			// Observation matrix
8
  init_int p					// Number of fixed effects
9
  init_matrix X(1,n,1,p)			// Design matrix for fixed effects
10
  init_int M					// Number of RE blocks (crossed terms)
11
  init_ivector q(1,M)				// Number of levels of the grouping variable per RE block; Can be skipped
12
  init_ivector m(1,M)				// Number of random effects parameters within each block
13
  int sum_mq					// sum(m*q), calculated below: should be read from R
14
  init_int ncolZ				
15
  init_matrix Z(1,n,1,ncolZ)			// Design matrix for random effects
16
  init_imatrix I(1,n,1,ncolZ)			// Index vectors into joint RE vector "u" for each
17
  init_ivector cor_flag(1,M)			// Indicator for whether each RE block should be correlated
18
  init_ivector cor_block_start(1,M) 		// Not used: remove
19
  init_ivector cor_block_stop(1,M) 		// Not used: remove
20
  init_int numb_cor_params			// Total number of correlation parameters to be estimated
21
  init_int like_type_flag   			// 0 poisson 1 binomial 2 negative binomial 3 Gamma 4 beta 5 gaussian 6 truncated poisson 7 trunc NB 8 logistic 9 betabinomial
22
  init_int link_type_flag   			// 0 log 1 logit 2 probit 3 inverse 4 cloglog 5 identity
23
  init_int rlinkflag                            // robust link function?
24
  init_int no_rand_flag   			// 0 have random effects 1 no random effects
25
  init_int zi_flag				// Zero inflation (zi) flag: 1=zi, 0=no zi
26
  // init_int zi_model_flag                        // ZI varies among groups/covariates?
27
  // init_matrix G(1,n,1,ncolG)                    // Design matrix for zero-inflation (fixed effects)
28
  // TESTING: remove eventually?
29
  init_int zi_kluge				// apply zi=0.001?
30
  init_int nbinom1_flag				// 1=NBinom1, 0=NBinom2
31
  init_int intermediate_maxfn			// Not used
32
  init_int has_offset				// Offset in linear predictor: 0=no offset, 1=with offset
33
  init_vector offset(1,n)			// Offset vector
34

  
35

  
36
  // Makes design matrix X orthogonal to improve numeric stability
37
  matrix rr(1,n,1,6)
38
  matrix phi(1,p,1,p)
39
 LOC_CALCS
40
  int i,j;
41
  phi.initialize();
42
  for (i=1;i<=p;i++)
43
  {
44
    phi(i,i)=1.0;
45
  }
46
  dmatrix trr=trans(rr);
47
  trr(6).fill_seqadd(1,1);
48
  rr=trans(trr);
49

  
50
  dmatrix TX(1,p,1,n);
51
  TX=trans(X);
52
  for (i=1;i<=p;i++)
53
  {
54
    double tmp=norm(TX(i));
55
    TX(i)/=tmp;
56
    phi(i)/=tmp;
57
    for (j=i+1;j<=p;j++)
58
    {
59
      double a=TX(j)*TX(i);
60
      TX(j)-=a*TX(i);
61
      phi(j)-=a*phi(i);
62
    }
63
  }
64
  X=trans(TX);
65

  
66
  sum_mq = 0;
67
  for (i=1;i<=M;i++)
68
    sum_mq += m(i)*q(i);
69

  
70
  ofstream ofs("phi.rep");
71
  for (i=1; i<=p; i++)
72
  {
73
    for (j=1; j<=p; j++)
74
    {
75
      ofs << phi(i,j) << " ";
76
    }
77
    ofs << endl;
78
  }
79
  ofs << endl;
80

  
81
INITIALIZATION_SECTION
82
 tmpL 1.0
83
 tmpL1 0.0
84
 log_alpha 1
85
 pz .001
86

  
87
PARAMETER_SECTION
88
 LOC_CALCS
89

  
90
  // BMB: FIXME: do we need this?  formerly disallowed for binomial (was like_type_flag 2);
91
  //     would be problematic for binary data but otherwise OK.  Should test in R code
92
  //  if(zi_flag && (like_type_flag>=2))
93
  // {
94
  //  cerr << "Zero inflation not allowed for this response type" << endl;
95
  //  ad_exit(1);
96
  // }
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff