DATA_SECTION init_int n // Number of observations init_vector y(1,n) // Observation vector init_int p // Number of fixed effects init_matrix X(1,n,1,p) // Covariate matrix for fixed effects init_int q // Number of Clusters init_int m // Number of random effects within clusters init_matrix Z(1,n,1,m) // Covariate matrix for random effects init_ivector group(1,q+1) // where the qth group ends init_int cor_flag init_int like_type_flag // 0 bvprobit 1 bvlogit init_int no_rand_flag // 0 have random effects 1 no random effects matrix rr(1,n,1,6) matrix phi(1,p,1,p) LOC_CALCS int i,j; phi.initialize(); for (i=1;i<=p;i++) { phi(i,i)=1.0; } dmatrix trr=trans(rr); trr(6).fill_seqadd(1,1); rr=trans(trr); dmatrix TX(1,p,1,n); TX=trans(X); for (i=1;i<=p;i++) { double tmp=norm(TX(i)); TX(i)/=tmp; phi(i)/=tmp; for (j=i+1;j<=p;j++) { double a=TX(j)*TX(i); TX(j)-=a*TX(i); phi(j)-=a*phi(i); } } X=trans(TX); INITIALIZATION_SECTION tmpL 1.0 tmpL1 0.0 PARAMETER_SECTION LOC_CALCS int rand_phase=3; int rand_sd_phase; int kludge_phase=-1; if (no_rand_flag==1) { rand_phase=-1; rand_sd_phase=-1; } else { rand_phase=2; rand_sd_phase=2; } int cor_phase=4; if (cor_flag && no_rand_flag==0) cor_phase=3; else cor_phase=-3; global_rand_phase=rand_phase; END_CALCS init_vector b(1,p,1) // Fixed effects sdreport_vector real_b(1,p) !! int mm = (m*(m+1))/2; !! int m1 = (m*(m-1))/2; init_bounded_vector tmpL(1,m,-10,10.5,rand_sd_phase) // Elements of cholesky-factor L of correlation matrix init_bounded_vector tmpL1(1,m1,-10,10.5,cor_phase) // Elements of cholesky-factor L of correlation matrix //sdreport_number sigma // sqrt(variance component) sdreport_matrix S(1,m,1,m) vector lambda(1,n) init_number kkludge(kludge_phase) random_effects_matrix u(1,q,1,m,rand_phase) // Random effects objective_function_value g // Log-likelihood PRELIMINARY_CALCS_SECTION cout << setprecision(4); PROCEDURE_SECTION g=0.0; int i,j; switch(like_type_flag) { case 0: // binomial probit for (j=1;j<=q;j++) { if (cor_flag) { bin_probit(j,tmpL,tmpL1,u(j),b); } else { bin_probit(j,tmpL,u(j),b); } } break; case 1: // logit for (j=1;j<=q;j++) { if (cor_flag) { bin_logit(j,tmpL,tmpL1,u(j),b); } else { bin_logit(j,tmpL,u(j),b); } } break; default: cerr << "Illegal value for like_type_flag" << endl; ad_exit(1); } if (sd_phase()) { real_b=b*phi; dvar_matrix L(1,m,1,m); L.initialize(); if (tmpL.indexmax() != m) { cerr << "size error in get_choleski" << endl; } int ii=1; if (cor_flag) { int ii=1; L(1,1)=1; for (i=1;i<=m;i++) { L(i,i)=1; for (int j=1;j #include int global_rand_phase; #ifndef PI const double PI = 3.14159265358979323846; #endif