// GLMM with negative binomial response // Example 2 in Booth, Casella, Friedl, Hobert (2003) Statistical Modelling, vol. 3, no. 3, pp. 179-191 // Numerical results meant to match rightmost column in Table 2 of Booth et al. 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_vector I(1,n) // Grouping variable init_int cor_flag init_int easy_flag init_int zi_flag init_int like_type_flag // 0 neg bin 1 poisson init_int no_rand_flag // 0 have random effects 1 no random effects init_int intermediate_maxfn 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); cout << I << endl; 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 log_alpha 1 pz .001 PARAMETER_SECTION LOC_CALCS int rand_phase=3; int rand_sd_phase; if (no_rand_flag==1) { rand_phase=-1; rand_sd_phase=-1; } else { rand_phase=3; rand_sd_phase=4; } int kludge_phase=-1; if (easy_flag) kludge_phase=-1; else kludge_phase=15; int zi_phase=5; if (zi_flag) zi_phase=5; else zi_phase=-5; int cor_phase=4; if (cor_flag && no_rand_flag==0) cor_phase=4; else cor_phase=-4; int alpha_phase=5; if(like_type_flag==0) alpha_phase=6; else alpha_phase=-6; global_rand_phase=rand_phase; END_CALCS init_bounded_number pz(.000001,0.999,zi_phase) 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 init_bounded_number log_alpha(-5.,6.,alpha_phase) // phase 2 by dave //sdreport_number sigma // sqrt(variance component) sdreport_number alpha // 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); BETWEEN_PHASES_SECTION maxfn=1000; if (current_phase()>global_rand_phase && global_rand_phase>0 ) { if (!last_phase()) { maxfn=intermediate_maxfn; } } PROCEDURE_SECTION g=0.0; int i; for (i=1;i<=q;i++) { fu(u(i)); } switch(like_type_flag) { case 0: // neg binomial if (easy_flag) { for (i=1;i<=n;i++) { if (cor_flag) { easy_f1_nb(i,tmpL,tmpL1,u((int) I(i)),b,log_alpha,pz); } else { easy_f1_nb(i,tmpL,u((int) I(i)),b,log_alpha,pz); } } } else { for (i=1;i<=n;i++) { if (cor_flag) { hard_f1_nb(i,tmpL,tmpL1,u((int) I(i)),b,log_alpha,pz); } else { hard_f1_nb(i,tmpL,u((int) I(i)),b,log_alpha,pz); } } } break; case 1: // poisson if (easy_flag) { for (i=1;i<=n;i++) { if (cor_flag) { easy_f1_poisson(i,tmpL,tmpL1,u((int) I(i)),b,pz); } else { easy_f1_poisson(i,tmpL,u((int) I(i)),b,pz); } } } else { for (i=1;i<=n;i++) { if (cor_flag) { hard_f1_poisson(i,tmpL,tmpL1,u((int) I(i)),b,pz); } else { hard_f1_poisson(i,tmpL,u((int) I(i)),b,pz); } } } break; default: cerr << "Illegal value for like_type_flag" << endl; ad_exit(1); } if (sd_phase()) { alpha = exp(log_alpha); real_b=b*phi; dvar_matrix L(1,m,1,m); L.initialize(); int i,j; if (tmpL.indexmax() != m) { cerr << "size error in get_choleski" << endl; } int ii=1; L(1,1)=1.0; for (i=2;i<=m;i++) { L(i,i)=1; for (int j=1;j #include #include int controlword=0; /* BOOL DllMain(HANDLE hinstDLL, DWORD reason, LPVOID res) { // for floating point so that R doesn't complain //controlword=control87(0,0); return TRUE; } */ #if defined(PRINT_DEBUG) #include ofstream ofs("model.out"); #endif int global_rand_phase; void normal_negbin_mixture(int i,double e,double pn,const prevariable& normal,const prevariable& pz,const prevariable& lambda, const prevariable& tau,double yi,prevariable& g) { if (yi==0) { g-=log(e+pn*exp(-normal)+pz+(1.0-pz) *exp(log_negbinomial_density(yi,lambda,tau))); } else { g-=log(e+pn*exp(-normal)+(1.0-pz)*exp(log_negbinomial_density(yi,lambda,tau))); } } void normal_negbin_mixture(int i,double e,double pn,const df1b2variable& normal,const df1b2variable& pz,const df1b2variable& lambda, const df1b2variable& tau,double yi,df1b2variable& g) { if (yi==0) { g-=log(e+pn*exp(-normal)+pz+(1.0-pz) *exp(log_negbinomial_density(yi,lambda,tau))); } else { g-=log(e+pn*exp(-normal)+(1.0-pz)*exp(log_negbinomial_density(yi,lambda,tau))); } } FINAL_SECTION //control87(controlword,0xfffff);