Hello, Is it possible to create normal.bayes models in Zelig that have random effects? Here's my application: I would like to create specifications for an electrical motor on how much electrical current they draw. Eighty samples from five batches of motor were measured, with different numbers of motors from each batch. So the batches are a random effect. I've been advised that a Bayesian analysis will handle the unbalanced data better. Batches 1 & 2 had 20 motors measured apiece, Batch 3 had 15 motors, Batch 4 had 4 motors, leaving the last motor from the 5th batch. My objective is to get an interval that will capture a tolerance range of the current measurements. Let's say an interval that predicts 99% of future Current Draw measurements. Does Zelig have this capability, I didn't find it in the few vignettes I read through. Regards, Paul Paul Prew | Statistician 651-795-5942   |   fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

Dear Gregor, thank you for your very helpful reply. I made some progress, but wonder if you have some advice for a sticking point I've reached. I apologize for not providing a self-contained example, that would take me a bit to figure out how to create it. Here's the steps I've followed: Unbalanced data (previous post was incorrect: 60 motors, not 80). Two batches/DateCodes = 20 motors each, 1 DateCode = 15 motors, 1 DateCode = 4 motors, leaving the last motor from the 5th batch/DateCode. 'data.frame': 60 obs. of 2 variables: $ Date.Code : Factor w/ 5 levels "22209","23109",..: 3 3 3 3 3 3 3 3 3 3 ... $ Current.NoLoad: num 0.17 0.18 0.23 0.21 0.22 0.22 0.23 0.25 0.2 0.19 ... I used the syntax below, that seems to parallel that used in the Dyestuff example in the lme4 documentation. model="ls.mixed", + data=AmpNoLoad.data) A quick look at batch-to-batch (DateCode-to-DateCode) mean differences list(Date.Code=AmpNoLoad.data$Date.Code), mean, na.rm=TRUE) Date.Code 22209 23109 26509 29408 30009 0.2070 0.2300 0.2173 0.2275 0.2145 However, in the mixed model, the Date.Code factor has no calculated Variance. I have only a dim view of what I'm talking about here, but am aware that REML employs shrinkage on random effects that appear to be weak. Has the Date.Code variance been shrunk to zero for the amount of precision? Below, I compare this REML analysis to the Bayes analysis performed using BugsXLA (an Excel-based interface to WinBUGS that I came across). --- So yes, I have my Bayesian analysis, with a distribution on the variance components --- but am still curious as to whether I'm using Zelig correctly. Also, the BugsXLA doesn't allow me to specify the width of the HPD for individual predictions, it just uses a default of 95%. I can't set specifications expecting 5% failures. Linear mixed model fit by REML Formula: Current.NoLoad ~ tag(1 | Date.Code) Data: AmpNoLoad.data AIC BIC logLik deviance REMLdev -232 -226 119 -247 -238 Random effects: Groups Name Variance Std.Dev. Date.Code (Intercept) 0.000000 0.0000 Residual 0.000966 0.0311 Number of obs: 60, groups: Date.Code, 5 Fixed effects: Estimate Std. Error t value (Intercept) 0.21383 0.00401 53.3 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ BugsXLA output Label Mean St.Dev. CONSTANT 0.2152 6.719E-3 SD(DateCode) 8.680E-3 7.732E-3 SD(residual) 0.0322 3.072E-3 Exchangeable (Random) Effects DateCode 26509 4.485E-4 6.944E-3 DateCode 29408 2.439E-3 8.460E-3 DateCode 23109 1.815E-3 0.0101 DateCode 22209 -4.105E-3 7.423E-3 DateCode 30009 -7.154E-4 6.845E-3 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ BACK TO R R version 2.9.1 (2009-06-26) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] tcltk stats graphics grDevices utils datasets methods base other attached packages: [1] lme4_0.999375-32 Matrix_0.999375-31 lattice_0.17-26 Zelig_3.4-5 boot_1.2-41 MASS_7.2-49 Rcmdr_1.5-3 [8] car_1.2-16 loaded via a namespace (and not attached): [1] grid_2.9.1 Paul Prew | Statistician 651-795-5942 | fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 ________________________________ From: gregor.gorjanc(a)gmail.com [mailto:gregor.gorjanc@gmail.com] On Behalf Of Gregor GORJANC Sent: Monday, December 21, 2009 3:34 PM To: Prew, Paul Cc: zelig(a)lists.gking.harvard.edu Subject: Re: [zelig] Bayesian normal models with random effects Paul, you can fit so called mixed effect models with Zelig via lme4 package. However, that is not full Bayesian analysis, as realised values of random effects are evaluated using the mode of posterior for variance components instead of full distribution, i.e., REML analysis. But, this is a first and easy step in the direction you want to go. gg 2009/12/21 Prew, Paul &lt;Paul.Prew(a)ecolab.com&gt; Hello, Is it possible to create normal.bayes models in Zelig that have random effects? Here's my application: I would like to create specifications for an electrical motor on how much electrical current they draw. Eighty samples from five batches of motor were measured, with different numbers of motors from each batch. So the batches are a random effect. I've been advised that a Bayesian analysis will handle the unbalanced data better. Batches 1 & 2 had 20 motors measured apiece, Batch 3 had 15 motors, Batch 4 had 4 motors, leaving the last motor from the 5th batch. My objective is to get an interval that will capture a tolerance range of the current measurements. Let's say an interval that predicts 99% of future Current Draw measurements. Does Zelig have this capability, I didn't find it in the few vignettes I read through. Regards, Paul Paul Prew | Statistician 651-795-5942 | fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/ -- -- Lep pozdrav / With regards, Gregor Gorjanc ---------------------------------------------------------------------- University of Ljubljana PhD student Biotechnical Faculty www: http://gregor.gorjanc.googlepages.com Zootechnical Department blog: http://ggorjan.blogspot.com Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si SI-1230 Domzale fax: +386 (0)1 72 17 888 Slovenia, Europe tel: +386 (0)1 72 17 861 ---------------------------------------------------------------------- CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

Dear Gregor, Thank you for pointing me to the half-Cauchy as a possibility for the prior on the between-group variance. I got one Gelman paper on that topic, and it helped me think through my choices. Further, your point about the small number of groups is well-taken, especially in light of how unbalanced they are. The engineer I'm advising assures me that there's more verification work to be done, right now we're floating these as trial specs to compare vendors. I'm sure you will sleep much more easily tonight knowing that. I may jump this topic over to the R-Mixed Models thread, to see if anyone wants to weigh in on how the mode of the between group variance could be zero, given the group averages below for the 5 DateCodes/lot numbers: Date.Code 22209 23109 26509 29408 30009 0.2070 0.2300 0.2173 0.2275 0.2145 By the way, BugsXLA is an Excel GUI that acts as a front-end for WinBUGs, written up in the Journal of Statistical Software. Regards, Paul Paul Prew | Statistician 651-795-5942 | fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 ________________________________ From: gregor.gorjanc(a)gmail.com [mailto:gregor.gorjanc@gmail.com] On Behalf Of Gregor GORJANC Sent: Wednesday, December 23, 2009 2:14 AM To: Prew, Paul Cc: zelig(a)lists.gking.harvard.edu Subject: Re: [zelig] Bayesian normal models with random effects Hi Paul, You have 5 groups and 60 records, right? You can not expect to get meaningful estimates of between group variance from only 5 groups - that is why Zelig (via lme4) reports zero variance. I guess that the mode of posterior for between group variance is zero, while there is some mass above - a rather long tail as REML estimates give the mode of joint posterior of variances where prior for variances is assumed uniform on (0,+inf). I am not familiar with what BugsXLA does, but the later fact (long tail) is likely the reason that you get larger estimates from BugsXLS, as mean is reported and not the mode. An important fact here is also a prior on between group variance as there is little information from the data (5 groups). If I would be in your place I would put group as fixed effect, i.e., whithout hierarchical prior. Another thing is to try Gelman's half-Cauchy prior for between group variance - look on his website or free internet Bayesian Analysis journal. Good luck! gg 2009/12/22 Prew, Paul &lt;Paul.Prew(a)ecolab.com&gt; Dear Gregor, thank you for your very helpful reply. I made some progress, but wonder if you have some advice for a sticking point I've reached. I apologize for not providing a self-contained example, that would take me a bit to figure out how to create it. Here's the steps I've followed: Unbalanced data (previous post was incorrect: 60 motors, not 80). Two batches/DateCodes = 20 motors each, 1 DateCode = 15 motors, 1 DateCode = 4 motors, leaving the last motor from the 5th batch/DateCode. 'data.frame': 60 obs. of 2 variables: $ Date.Code : Factor w/ 5 levels "22209","23109",..: 3 3 3 3 3 3 3 3 3 3 ... $ Current.NoLoad: num 0.17 0.18 0.23 0.21 0.22 0.22 0.23 0.25 0.2 0.19 ... I used the syntax below, that seems to parallel that used in the Dyestuff example in the lme4 documentation. model="ls.mixed", + data=AmpNoLoad.data) A quick look at batch-to-batch (DateCode-to-DateCode) mean differences list(Date.Code=AmpNoLoad.data$Date.Code), mean, na.rm=TRUE) Date.Code 22209 23109 26509 29408 30009 0.2070 0.2300 0.2173 0.2275 0.2145 However, in the mixed model, the Date.Code factor has no calculated Variance. I have only a dim view of what I'm talking about here, but am aware that REML employs shrinkage on random effects that appear to be weak. Has the Date.Code variance been shrunk to zero for the amount of precision? Below, I compare this REML analysis to the Bayes analysis performed using BugsXLA (an Excel-based interface to WinBUGS that I came across). --- So yes, I have my Bayesian analysis, with a distribution on the variance components --- but am still curious as to whether I'm using Zelig correctly. Also, the BugsXLA doesn't allow me to specify the width of the HPD for individual predictions, it just uses a default of 95%. I can't set specifications expecting 5% failures. Linear mixed model fit by REML Formula: Current.NoLoad ~ tag(1 | Date.Code) Data: AmpNoLoad.data AIC BIC logLik deviance REMLdev -232 -226 119 -247 -238 Random effects: Groups Name Variance Std.Dev. Date.Code (Intercept) 0.000000 0.0000 Residual 0.000966 0.0311 Number of obs: 60, groups: Date.Code, 5 Fixed effects: Estimate Std. Error t value (Intercept) 0.21383 0.00401 53.3 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ BugsXLA output Label Mean St.Dev. CONSTANT 0.2152 6.719E-3 SD(DateCode) 8.680E-3 7.732E-3 SD(residual) 0.0322 3.072E-3 Exchangeable (Random) Effects DateCode 26509 4.485E-4 6.944E-3 DateCode 29408 2.439E-3 8.460E-3 DateCode 23109 1.815E-3 0.0101 DateCode 22209 -4.105E-3 7.423E-3 DateCode 30009 -7.154E-4 6.845E-3 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ BACK TO R R version 2.9.1 (2009-06-26) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] tcltk stats graphics grDevices utils datasets methods base other attached packages: [1] lme4_0.999375-32 Matrix_0.999375-31 lattice_0.17-26 Zelig_3.4-5 boot_1.2-41 MASS_7.2-49 Rcmdr_1.5-3 [8] car_1.2-16 loaded via a namespace (and not attached): [1] grid_2.9.1 Paul Prew | Statistician 651-795-5942 | fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 ________________________________ From: gregor.gorjanc(a)gmail.com [mailto:gregor.gorjanc@gmail.com] On Behalf Of Gregor GORJANC Sent: Monday, December 21, 2009 3:34 PM To: Prew, Paul Cc: zelig(a)lists.gking.harvard.edu Subject: Re: [zelig] Bayesian normal models with random effects Paul, you can fit so called mixed effect models with Zelig via lme4 package. However, that is not full Bayesian analysis, as realised values of random effects are evaluated using the mode of posterior for variance components instead of full distribution, i.e., REML analysis. But, this is a first and easy step in the direction you want to go. gg 2009/12/21 Prew, Paul &lt;Paul.Prew(a)ecolab.com&gt; Hello, Is it possible to create normal.bayes models in Zelig that have random effects? Here's my application: I would like to create specifications for an electrical motor on how much electrical current they draw. Eighty samples from five batches of motor were measured, with different numbers of motors from each batch. So the batches are a random effect. I've been advised that a Bayesian analysis will handle the unbalanced data better. Batches 1 & 2 had 20 motors measured apiece, Batch 3 had 15 motors, Batch 4 had 4 motors, leaving the last motor from the 5th batch. My objective is to get an interval that will capture a tolerance range of the current measurements. Let's say an interval that predicts 99% of future Current Draw measurements. Does Zelig have this capability, I didn't find it in the few vignettes I read through. Regards, Paul Paul Prew | Statistician 651-795-5942 | fax 651-204-7504 Ecolab Research Center | Mail Stop ESC-F4412-A 655 Lone Oak Drive | Eagan, MN 55121-1560 CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/ -- -- Lep pozdrav / With regards, Gregor Gorjanc ---------------------------------------------------------------------- University of Ljubljana PhD student Biotechnical Faculty www: http://gregor.gorjanc.googlepages.com Zootechnical Department blog: http://ggorjan.blogspot.com Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si SI-1230 Domzale fax: +386 (0)1 72 17 888 Slovenia, Europe tel: +386 (0)1 72 17 861 ---------------------------------------------------------------------- CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/ -- -- Lep pozdrav / With regards, Gregor Gorjanc ---------------------------------------------------------------------- University of Ljubljana PhD student Biotechnical Faculty www: http://gregor.gorjanc.googlepages.com Zootechnical Department blog: http://ggorjan.blogspot.com Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si SI-1230 Domzale fax: +386 (0)1 72 17 888 Slovenia, Europe tel: +386 (0)1 72 17 861 ---------------------------------------------------------------------- CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message.