------------------------------------------------------------------------------------------------------------ name: log: c:\vhm812-data\l11b_glmm_discrete.txt log type: text opened on: 18 Mar 2016, 10:39:42 . . *Random effects binary data . *Pig respiratory disease data . use pig_adg.dta, clear (Growth Performance and Abattoir Surveillance Data on Pigs, Real Data) . * generate dichotomous atrophic rhinits variable . egen ar_g1=cut(ar), at(0, 1.5, 99) icodes . . * 2x2 table analysis . cc pn ar_g1 Proportion | Exposed Unexposed | Total Exposed -----------------+------------------------+------------------------ Cases | 109 77 | 186 0.5860 Controls | 66 89 | 155 0.4258 -----------------+------------------------+------------------------ Total | 175 166 | 341 0.5132 | | | Point estimate | [95% Conf. Interval] |------------------------+------------------------ Odds ratio | 1.908894 | 1.21155 3.009556 (exact) Attr. frac. ex. | .4761365 | .1746111 .6677251 (exact) Attr. frac. pop | .2790262 | +------------------------------------------------- chi2(1) = 8.69 Pr>chi2 = 0.0032 . . * ordinary logistic regression . logit pn ar_g1 Iteration 0: log likelihood = -234.95215 Iteration 1: log likelihood = -230.59287 Iteration 2: log likelihood = -230.59173 Iteration 3: log likelihood = -230.59173 Logistic regression Number of obs = 341 LR chi2(1) = 8.72 Prob > chi2 = 0.0031 Log likelihood = -230.59173 Pseudo R2 = 0.0186 ------------------------------------------------------------------------------ pn | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | .6465241 .2203379 2.93 0.003 .2146697 1.078378 _cons | -.1448309 .1556373 -0.93 0.352 -.4498744 .1602125 ------------------------------------------------------------------------------ . predict pa_noclus (option pr assumed; Pr(pn)) . . *data from 15 farms - outcome varies from 17% to 96% . preserve . collapse (sum) pn (count) total=ar_g1, by(farm) . gen prev=pn/total . summ prev Variable | Obs Mean Std. Dev. Min Max -------------+--------------------------------------------------------- prev | 15 .5493906 .2270394 .1785714 .9583333 . graph bar prev, over(farm) . restore . . * random effect logistic regression . melogit pn ar_g1 || farm: Fitting fixed-effects model: Iteration 0: log likelihood = -230.79889 Iteration 1: log likelihood = -230.59177 Iteration 2: log likelihood = -230.59173 Iteration 3: log likelihood = -230.59173 Refining starting values: Grid node 0: log likelihood = -213.91338 Fitting full model: Iteration 0: log likelihood = -213.91338 Iteration 1: log likelihood = -213.52 Iteration 2: log likelihood = -213.51177 Iteration 3: log likelihood = -213.51176 Mixed-effects logistic regression Number of obs = 341 Group variable: farm Number of groups = 15 Obs per group: min = 14 avg = 22.7 max = 28 Integration method: mvaghermite Integration pts. = 7 Wald chi2(1) = 2.86 Log likelihood = -213.51176 Prob > chi2 = 0.0905 ------------------------------------------------------------------------------ pn | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | .4369322 .2581451 1.69 0.091 -.0690229 .9428872 _cons | .0196548 .3009262 0.07 0.948 -.5701497 .6094593 -------------+---------------------------------------------------------------- farm | var(_cons)| .8774246 .4325507 .333877 2.305861 ------------------------------------------------------------------------------ LR test vs. logistic model: chibar2(01) = 34.16 Prob >= chibar2 = 0.0000 . melogit pn ar_g1 || farm:, or Fitting fixed-effects model: Iteration 0: log likelihood = -230.79889 Iteration 1: log likelihood = -230.59177 Iteration 2: log likelihood = -230.59173 Iteration 3: log likelihood = -230.59173 Refining starting values: Grid node 0: log likelihood = -213.91338 Fitting full model: Iteration 0: log likelihood = -213.91338 Iteration 1: log likelihood = -213.52 Iteration 2: log likelihood = -213.51177 Iteration 3: log likelihood = -213.51176 Mixed-effects logistic regression Number of obs = 341 Group variable: farm Number of groups = 15 Obs per group: min = 14 avg = 22.7 max = 28 Integration method: mvaghermite Integration pts. = 7 Wald chi2(1) = 2.86 Log likelihood = -213.51176 Prob > chi2 = 0.0905 ------------------------------------------------------------------------------ pn | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | 1.547951 .3995959 1.69 0.091 .9333053 2.567383 _cons | 1.019849 .3068993 0.07 0.948 .5654408 1.839437 -------------+---------------------------------------------------------------- farm | var(_cons)| .8774246 .4325507 .333877 2.305861 ------------------------------------------------------------------------------ LR test vs. logistic model: chibar2(01) = 34.16 Prob >= chibar2 = 0.0000 . . *ICC // approximation to real ICC by pi^2/3 = 3.29 . estat icc Residual intraclass correlation ------------------------------------------------------------------------------ Level | ICC Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ farm | .2105503 .0819422 .0921359 .4120752 ------------------------------------------------------------------------------ . . **Predictions - not in the notes . * prob. for a pig from a specific farm . predict re_cs, mu (predictions based on fixed effects and posterior means of random effects) (using 7 quadrature points) . . * prediction random effects . predict u_farm, reffects (calculating posterior means of random effects) (using 7 quadrature points) . scatter u_farm farm, xlabel(1(1)15) . . * PA logistic regression . xtgee pn ar_g1, fam(binomial) link(logit) i(farm) robust Iteration 1: tolerance = .17657314 Iteration 2: tolerance = .0013339 Iteration 3: tolerance = .00001851 Iteration 4: tolerance = 2.251e-07 GEE population-averaged model Number of obs = 341 Group variable: farm Number of groups = 15 Link: logit Obs per group: Family: binomial min = 14 Correlation: exchangeable avg = 22.7 max = 28 Wald chi2(1) = 2.71 Scale parameter: 1 Prob > chi2 = 0.1000 (Std. Err. adjusted for clustering on farm) ------------------------------------------------------------------------------ | Robust pn | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | .3539583 .2151855 1.64 0.100 -.0677976 .7757142 _cons | .0183926 .2714184 0.07 0.946 -.5135777 .5503629 ------------------------------------------------------------------------------ . predict pa, mu /*mean prob from any farm */ . . * Random effects models for count data . * open the tb_real dataset . use tb_real, clear . * Poisson model with no random effects . glm reactors i.type i.sex i.age, exp(par) link(log) fam(poisson) Iteration 0: log likelihood = -338.47585 Iteration 1: log likelihood = -246.77274 Iteration 2: log likelihood = -239.32758 Iteration 3: log likelihood = -238.67153 Iteration 4: log likelihood = -238.662 Iteration 5: log likelihood = -238.662 Generalized linear models No. of obs = 134 Optimization : ML Residual df = 127 Scale parameter = 1 Deviance = 348.3526568 (1/df) Deviance = 2.742934 Pearson = 1105.702579 (1/df) Pearson = 8.70632 Variance function: V(u) = u [Poisson] Link function : g(u) = ln(u) [Log] AIC = 3.666597 Log likelihood = -238.6620031 BIC = -273.673 ------------------------------------------------------------------------------ | OIM reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | .4422157 .2364116 1.87 0.061 -.0211426 .905574 cervid | 1.066241 .2333805 4.57 0.000 .6088238 1.523658 other | .4383852 .6149186 0.71 0.476 -.7668331 1.643603 | sex | male | -.3618833 .1953974 -1.85 0.064 -.7448552 .0210886 | age | 12-24 mo | 2.673416 .721739 3.70 0.000 1.258834 4.087999 >24 mont | 2.601192 .7136354 3.64 0.000 1.202493 3.999892 | _cons | -11.68987 .7397843 -15.80 0.000 -13.13982 -10.23992 ln(par) | 1 (exposure) ------------------------------------------------------------------------------ . estimates store pois . . * negative binomial model no random effects . glm reactors i.type i.sex i.age, exp(par) link(log) fam(nbin ml) Iteration 0: log likelihood = -168.6314 Iteration 1: log likelihood = -157.89374 Iteration 2: log likelihood = -157.73611 Iteration 3: log likelihood = -157.73596 Iteration 4: log likelihood = -157.73596 Generalized linear models No. of obs = 134 Optimization : ML Residual df = 127 Scale parameter = 1 Deviance = 99.35977885 (1/df) Deviance = .7823605 Pearson = 374.8585955 (1/df) Pearson = 2.951642 Variance function: V(u) = u+(1.7404)u^2 [Neg. Binomial] Link function : g(u) = ln(u) [Log] AIC = 2.458746 Log likelihood = -157.735955 BIC = -522.6659 ------------------------------------------------------------------------------ | OIM reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | .6046115 .674193 0.90 0.370 -.7167824 1.926005 cervid | .6657202 .6821121 0.98 0.329 -.671195 2.002635 other | .8000348 1.11801 0.72 0.474 -1.391225 2.991295 | sex | male | -.0574791 .4044774 -0.14 0.887 -.8502403 .735282 | age | 12-24 mo | 2.253036 .8976472 2.51 0.012 .4936798 4.012392 >24 mont | 2.480955 .8779162 2.83 0.005 .7602705 4.201639 | _cons | -11.18145 1.052762 -10.62 0.000 -13.24482 -9.118069 ln(par) | 1 (exposure) ------------------------------------------------------------------------------ Note: Negative binomial parameter estimated via ML and treated as fixed once estimated. . estimates store nb . * Pearson dispersion parameter still large (2.95) . . * Poisson model with normal distributed random effects . mepoisson reactors i.type i.sex i.age, exp(par) || farm_id: Fitting fixed-effects model: Iteration 0: log likelihood = -602.90391 Iteration 1: log likelihood = -247.78547 Iteration 2: log likelihood = -239.70423 Iteration 3: log likelihood = -238.67796 Iteration 4: log likelihood = -238.66201 Iteration 5: log likelihood = -238.662 Refining starting values: Grid node 0: log likelihood = -150.38337 Fitting full model: Iteration 0: log likelihood = -150.38337 Iteration 1: log likelihood = -144.07026 Iteration 2: log likelihood = -143.57132 Iteration 3: log likelihood = -143.55978 Iteration 4: log likelihood = -143.55988 Iteration 5: log likelihood = -143.55988 Mixed-effects Poisson regression Number of obs = 134 Group variable: farm_id Number of groups = 30 Obs per group: min = 1 avg = 4.5 max = 13 Integration method: mvaghermite Integration pts. = 7 Wald chi2(6) = 18.42 Log likelihood = -143.55988 Prob > chi2 = 0.0053 ------------------------------------------------------------------------------ reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- | type | beef | -.3941614 .3320059 -1.19 0.235 -1.044881 .2565583 cervid | -.2394465 .4847969 -0.49 0.621 -1.189631 .7107379 other | -.1058558 .7985319 -0.13 0.895 -1.67095 1.459238 | sex | male | -.3392105 .2083118 -1.63 0.103 -.7474941 .0690731 | age | 12-24 mo | 2.716639 .7470885 3.64 0.000 1.252372 4.180906 >24 mont | 2.466631 .7256478 3.40 0.001 1.044388 3.888875 | _cons | -11.05318 .8283607 -13.34 0.000 -12.67674 -9.429625 ln(par) | 1 (exposure) -------------+---------------------------------------------------------------- farm_id | var(_cons)| 1.697739 .6022648 .8470544 3.402754 ------------------------------------------------------------------------------ LR test vs. Poisson model: chibar2(01) = 190.20 Prob >= chibar2 = 0.0000 . mepoisson reactors i.type i.sex i.age, exp(par) irr || farm_id: Fitting fixed-effects model: Iteration 0: log likelihood = -602.90391 Iteration 1: log likelihood = -247.78547 Iteration 2: log likelihood = -239.70423 Iteration 3: log likelihood = -238.67796 Iteration 4: log likelihood = -238.66201 Iteration 5: log likelihood = -238.662 Refining starting values: Grid node 0: log likelihood = -150.38337 Fitting full model: Iteration 0: log likelihood = -150.38337 Iteration 1: log likelihood = -144.07026 Iteration 2: log likelihood = -143.57132 Iteration 3: log likelihood = -143.55978 Iteration 4: log likelihood = -143.55988 Iteration 5: log likelihood = -143.55988 Mixed-effects Poisson regression Number of obs = 134 Group variable: farm_id Number of groups = 30 Obs per group: min = 1 avg = 4.5 max = 13 Integration method: mvaghermite Integration pts. = 7 Wald chi2(6) = 18.42 Log likelihood = -143.55988 Prob > chi2 = 0.0053 ------------------------------------------------------------------------------ reactors | IRR Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- | type | beef | .6742453 .2238534 -1.19 0.235 .3517337 1.292474 cervid | .7870634 .3815659 -0.49 0.621 .3043336 2.035493 other | .8995544 .7183229 -0.13 0.895 .1880684 4.30268 | sex | male | .7123325 .1483872 -1.63 0.103 .4735517 1.071514 | age | 12-24 mo | 15.12939 11.30299 3.64 0.000 3.498633 65.42507 >24 mont | 11.78269 8.550082 3.40 0.001 2.841659 48.85589 | _cons | .0000158 .0000131 -13.34 0.000 3.12e-06 .0000803 ln(par) | 1 (exposure) -------------+---------------------------------------------------------------- farm_id | var(_cons)| 1.697739 .6022648 .8470544 3.402754 ------------------------------------------------------------------------------ LR test vs. Poisson model: chibar2(01) = 190.20 Prob >= chibar2 = 0.0000 . . * Negative binomial with Normal random effects . meglm reactors i.type i.sex i.age, exp(par) || farm_id:, fam(nbin) Fitting fixed-effects model: Iteration 0: log likelihood = -175.65303 Iteration 1: log likelihood = -158.06607 Iteration 2: log likelihood = -157.73916 Iteration 3: log likelihood = -157.73596 Iteration 4: log likelihood = -157.73596 Refining starting values: Grid node 0: log likelihood = -153.88154 Fitting full model: Iteration 0: log likelihood = -153.88154 Iteration 1: log likelihood = -147.66002 Iteration 2: log likelihood = -144.59242 Iteration 3: log likelihood = -143.77217 Iteration 4: log likelihood = -143.59593 Iteration 5: log likelihood = -143.5688 Iteration 6: log likelihood = -143.56194 Iteration 7: log likelihood = -143.56037 Iteration 8: log likelihood = -143.55999 Iteration 9: log likelihood = -143.55991 Iteration 10: log likelihood = -143.55989 Iteration 11: log likelihood = -143.55988 Mixed-effects GLM Number of obs = 134 Family: negative binomial Link: log Overdispersion: mean Group variable: farm_id Number of groups = 30 Obs per group: min = 1 avg = 4.5 max = 13 Integration method: mvaghermite Integration pts. = 7 Wald chi2(6) = 18.42 Log likelihood = -143.55988 Prob > chi2 = 0.0053 ------------------------------------------------------------------------------ reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- | type | beef | -.3941584 .3320062 -1.19 0.235 -1.044879 .2565617 cervid | -.2394449 .484797 -0.49 0.621 -1.18963 .7107398 other | -.1058542 .798532 -0.13 0.895 -1.670948 1.45924 | sex | male | -.3392105 .2083121 -1.63 0.103 -.7474948 .0690738 | age | 12-24 mo | 2.716639 .7470886 3.64 0.000 1.252372 4.180905 >24 mont | 2.466631 .7256478 3.40 0.001 1.044388 3.888875 | _cons | -11.05318 .8283609 -13.34 0.000 -12.67674 -9.429626 ln(par) | 1 (exposure) -------------+---------------------------------------------------------------- /lnalpha | -15.62553 1798.877 -0.01 0.993 -3541.36 3510.109 -------------+---------------------------------------------------------------- farm_id | var(_cons)| 1.697739 .6022654 .8470541 3.402757 ------------------------------------------------------------------------------ LR test vs. nbinomial model: chibar2(01) = 28.35 Prob >= chibar2 = 0.0000 . * the overdispersion lnalpha is non signfincant which suggests that the . * overdispersion is due to clustering since this model is exactly . * the same as the Poisson with random effects . . estimates store pois_norm . estimates table pois nb pois_norm, se(%4.3f) b(%4.3f) -------------------------------------------- Variable | pois nb pois_~m -------------+------------------------------ reactors | type | beef | 0.442 0.605 -0.394 | 0.236 0.674 0.332 cervid | 1.066 0.666 -0.239 | 0.233 0.682 0.485 other | 0.438 0.800 -0.106 | 0.615 1.118 0.799 | sex | male | -0.362 -0.057 -0.339 | 0.195 0.404 0.208 | age | 12-24 mo | 2.673 2.253 2.717 | 0.722 0.898 0.747 >24 mont | 2.601 2.481 2.467 | 0.714 0.878 0.726 | _cons | -11.690 -11.181 -11.053 | 0.740 1.053 0.828 -------------+------------------------------ lnalpha | _cons | -15.626 | 1798.877 -------------+------------------------------ var(_cons[~) | _cons | 1.698 | 0.602 -------------------------------------------- legend: b/se . estimates stats pois nb pois_norm Akaike's information criterion and Bayesian information criterion ----------------------------------------------------------------------------- Model | Obs ll(null) ll(model) df AIC BIC -------------+--------------------------------------------------------------- pois | 134 . -238.662 7 491.324 511.6089 nb | 134 . -157.736 7 329.4719 349.7568 pois_norm | 134 . -143.5599 9 305.1198 331.2003 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note. . . **PA Poisson model . xtgee reactors i.type i.sex i.age, exp(par) i(farm_id) family(poisson) link(log) robust Iteration 1: tolerance = .05719918 Iteration 2: tolerance = .00962239 Iteration 3: tolerance = .00017782 Iteration 4: tolerance = .0000172 Iteration 5: tolerance = 7.799e-07 GEE population-averaged model Number of obs = 134 Group variable: farm_id Number of groups = 30 Link: log Obs per group: Family: Poisson min = 1 Correlation: exchangeable avg = 4.5 max = 13 Wald chi2(6) = 14.21 Scale parameter: 1 Prob > chi2 = 0.0274 (Std. Err. adjusted for clustering on farm_id) ------------------------------------------------------------------------------ | Robust reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | .5192248 .6667257 0.78 0.436 -.7875335 1.825983 cervid | 1.139405 .6818041 1.67 0.095 -.1969064 2.475717 other | .4687865 .624348 0.75 0.453 -.7549131 1.692486 | sex | male | -.4109311 .2359891 -1.74 0.082 -.8734613 .0515991 | age | 12-24 mo | 2.510354 .8547058 2.94 0.003 .8351612 4.185546 >24 mont | 2.442831 .9865273 2.48 0.013 .5092732 4.376389 | _cons | -11.62611 .9211882 -12.62 0.000 -13.4316 -9.820612 ln(par) | 1 (exposure) ------------------------------------------------------------------------------ . end of do-file . exit, clear