----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- name: log: c:\vhm812-data\l11b-cluster_discrete.txt log type: text opened on: 23 Mar 2015, 10:31:55 . . *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)) . . * random effect logistic regression . meqrlogit pn ar_g1 || farm: Refining starting values: Iteration 0: log likelihood = -213.91344 Iteration 1: log likelihood = -213.51368 Iteration 2: log likelihood = -213.5118 Performing gradient-based optimization: Iteration 0: log likelihood = -213.5118 Iteration 1: log likelihood = -213.5118 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 points = 7 Wald chi2(1) = 2.86 Log likelihood = -213.5118 Prob > chi2 = 0.0905 ------------------------------------------------------------------------------ pn | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | .4369353 .2581461 1.69 0.091 -.0690217 .9428923 _cons | .0196485 .3009417 0.07 0.948 -.5701864 .6094834 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ farm: Identity | var(_cons) | .8773258 .4324896 .3338486 2.305537 ------------------------------------------------------------------------------ LR test vs. logistic regression: chibar2(01) = 34.16 Prob>=chibar2 = 0.0000 . meqrlogit pn ar_g1 || farm:, or Refining starting values: Iteration 0: log likelihood = -213.91344 Iteration 1: log likelihood = -213.51368 Iteration 2: log likelihood = -213.5118 Performing gradient-based optimization: Iteration 0: log likelihood = -213.5118 Iteration 1: log likelihood = -213.5118 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 points = 7 Wald chi2(1) = 2.86 Log likelihood = -213.5118 Prob > chi2 = 0.0905 ------------------------------------------------------------------------------ pn | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- ar_g1 | 1.547956 .3995987 1.69 0.091 .9333065 2.567396 _cons | 1.019843 .3069132 0.07 0.948 .56542 1.839481 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ farm: Identity | var(_cons) | .8773258 .4324896 .3338486 2.305537 ------------------------------------------------------------------------------ LR test vs. logistic regression: chibar2(01) = 34.16 Prob>=chibar2 = 0.0000 . . **Predictions - not for the notes . * prob. for a pig from a specific farm . predict re_cs, mu . . * prediction random effects . predict u_farm, reffects relevel(farm) . 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: min = 14 Family: binomial avg = 22.7 Correlation: exchangeable 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 . gen logpar=log(par) . * Poisson model with no random effects . glm reactors i.type i.sex i.age, off(logpar) link(log) fam(poisson) Iteration 0: log likelihood = -338.47586 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.3526588 (1/df) Deviance = 2.742934 Pearson = 1105.70265 (1/df) Pearson = 8.70632 Variance function: V(u) = u [Poisson] Link function : g(u) = ln(u) [Log] AIC = 3.666597 Log likelihood = -238.6620041 BIC = -273.673 ------------------------------------------------------------------------------ | OIM reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | .4422158 .2364116 1.87 0.061 -.0211424 .9055741 cervid | 1.066241 .2333805 4.57 0.000 .608824 1.523659 other | .4383856 .6149186 0.71 0.476 -.7668327 1.643604 | sex | male | -.3618834 .1953974 -1.85 0.064 -.7448553 .0210885 | 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.13983 -10.23992 logpar | 1 (offset) ------------------------------------------------------------------------------ . estimates store pois . . * negative binomial model no random effects . glm reactors i.type i.sex i.age, off(logpar) 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.73595 Iteration 4: log likelihood = -157.73595 Generalized linear models No. of obs = 134 Optimization : ML Residual df = 127 Scale parameter = 1 Deviance = 99.35977922 (1/df) Deviance = .7823605 Pearson = 374.8586498 (1/df) Pearson = 2.951643 Variance function: V(u) = u+(1.7404)u^2 [Neg. Binomial] Link function : g(u) = ln(u) [Log] AIC = 2.458746 Log likelihood = -157.7359545 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 | .6657203 .6821121 0.98 0.329 -.6711948 2.002635 other | .800035 1.11801 0.72 0.474 -1.391225 2.991295 | sex | male | -.057479 .4044774 -0.14 0.887 -.8502402 .7352821 | age | 12-24 mo | 2.253036 .8976472 2.51 0.012 .4936796 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 logpar | 1 (offset) ------------------------------------------------------------------------------ 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 . meqrpoisson reactors i.type i.sex i.age, off(logpar) || farm_id: Refining starting values: Iteration 0: log likelihood = -150.383 Iteration 1: log likelihood = -143.9602 Iteration 2: log likelihood = -143.57299 Performing gradient-based optimization: Iteration 0: log likelihood = -143.57299 Iteration 1: log likelihood = -143.55957 Iteration 2: log likelihood = -143.55956 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 points = 7 Wald chi2(6) = 18.41 Log likelihood = -143.55956 Prob > chi2 = 0.0053 ------------------------------------------------------------------------------ reactors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | -.3941065 .3326697 -1.18 0.236 -1.046127 .2579141 cervid | -.2383307 .4867283 -0.49 0.624 -1.192301 .7156392 other | -.104708 .800331 -0.13 0.896 -1.673328 1.463912 | sex | male | -.3392047 .2083168 -1.63 0.103 -.747498 .0690887 | age | 12-24 mo | 2.716581 .7471013 3.64 0.000 1.252289 4.180872 >24 mont | 2.46664 .7256875 3.40 0.001 1.044319 3.888962 | _cons | -11.05516 .8300711 -13.32 0.000 -12.68207 -9.428252 logpar | 1 (offset) ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ farm_id: Identity | var(_cons) | 1.689432 .5943486 .8477829 3.36664 ------------------------------------------------------------------------------ LR test vs. Poisson regression: chibar2(01) = 190.20 Prob>=chibar2 = 0.0000 . meqrpoisson reactors i.type i.sex i.age, off(logpar) irr || farm_id: Refining starting values: Iteration 0: log likelihood = -150.383 Iteration 1: log likelihood = -143.9602 Iteration 2: log likelihood = -143.57299 Performing gradient-based optimization: Iteration 0: log likelihood = -143.57299 Iteration 1: log likelihood = -143.55957 Iteration 2: log likelihood = -143.55956 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 points = 7 Wald chi2(6) = 18.41 Log likelihood = -143.55956 Prob > chi2 = 0.0053 ------------------------------------------------------------------------------ reactors | IRR Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- type | beef | .6742823 .2243133 -1.18 0.236 .3512957 1.294228 cervid | .7879421 .3835137 -0.49 0.624 .3035222 2.045494 other | .9005874 .720768 -0.13 0.896 .1876216 4.322837 | sex | male | .7123366 .1483917 -1.63 0.103 .4735499 1.071531 | age | 12-24 mo | 15.12851 11.30253 3.64 0.000 3.498342 65.4229 >24 mont | 11.78279 8.550626 3.40 0.001 2.841462 48.86013 | _cons | .0000158 .0000131 -13.32 0.000 3.11e-06 .0000804 logpar | 1 (offset) ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ farm_id: Identity | var(_cons) | 1.689432 .5943486 .8477829 3.36664 ------------------------------------------------------------------------------ LR test vs. Poisson regression: chibar2(01) = 190.20 Prob>=chibar2 = 0.0000 . . 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.236 0.674 cervid | 1.066 0.666 | 0.233 0.682 other | 0.438 0.800 | 0.615 1.118 | sex | male | -0.362 -0.057 | 0.195 0.404 | age | 12-24 mo | 2.673 2.253 | 0.722 0.898 >24 mont | 2.601 2.481 | 0.714 0.878 | _cons | -11.690 -11.181 | 0.740 1.053 -------------+------------------------------ eq1 | | type | beef | -0.394 | 0.333 cervid | -0.238 | 0.487 other | -0.105 | 0.800 | sex | male | -0.339 | 0.208 | age | 12-24 mo | 2.717 | 0.747 >24 mont | 2.467 | 0.726 | _cons | -11.055 | 0.830 -------------+------------------------------ lns1_1_1 | _cons | 0.262 | 0.176 -------------------------------------------- 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.5596 8 303.1191 326.3018 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note . . **PA Poisson model . xtgee reactors i.type i.sex i.age, off(logpar) 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: min = 1 Family: Poisson avg = 4.5 Correlation: exchangeable 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 | .5192249 .6667257 0.78 0.436 -.7875335 1.825983 cervid | 1.139405 .6818042 1.67 0.095 -.1969062 2.475717 other | .468787 .6243481 0.75 0.453 -.7549128 1.692487 | sex | male | -.4109312 .2359892 -1.74 0.082 -.8734615 .0515992 | age | 12-24 mo | 2.510354 .8547058 2.94 0.003 .8351611 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 logpar | 1 (offset) ------------------------------------------------------------------------------ . end of do-file . exit, clear