Hugo, thank you for your patch! Additionally, I will be reviewing the submitted bugs, and hopefully will have a fix in the very near future. Cheers - Matt Owen On Nov 6, 2011, at 3:53 PM, Hugo Mildenberger wrote: - 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/

This is a result a change in VGAM that actually was just patched today. This morning, Zelig incremented to version 3.5-2 (for Linux/unix and Windows), in order to be consistent with the VGAM changes. You can find the new version at CRAN: http://cran.r-project.org/web/packages/Zelig/index.html On Nov 8, 2011, at 2:05 PM, Hugo Mildenberger wrote: > Matt, > > please skip also version 2. It is still not the end of the story. I've > encountered another severe problem when testing the mixed model examples > in Zelig's inst/doc directory with bootstrap=T. When running, as I > recall, the example from poisson.mixed.Rnw, I got the following: > > > Error in mvrnorm(1, 0, V.beta) : incompatible arguments > Calls: sim ... sim.default -> boot -> statistic -> param -> > param.mer -> mvrnorm > > The mvrnorm call is again in R/param.mixed.R (using original code here) > > 30 # sample random effects > 31 for (m in 1:n.G){ > 32 V.beta <- VarCorr(object)[[m]] > 33 gammas[[m]]<- mvrnorm(1, 0, V.beta) > 34 } > > > This had to do with random effect term tag(player - 1| month). In the > homerun dataset, there are two players. I concluded that line 33 should > probably written as: > > gammas[[name]] <- mvrnorm(1,rep(0,dim(vars.m)[1]),V.beta) > > Even so this gives apparently acceptable results, I'm a rather unsure > about this. > > Another problem is with gamma.mixed.Rnw. lmer/glmer had deprecated the > method argument and also dropped support for PQL, but introduced nAGQ as > a new argument for use with glmer. Anyway, whatever you try, the model > does not converge any more. I've reported this to R-sig-ME, but there > was no response until now. > > Finally, I tried to run R CMD check with the original version. To my > knowledge, the check command must complete successfully for a package to > be included in the cran repository. But the check command failed in just > another way: > > > > I've also attached the complete output of "R CMD check > Zelig_3.5.1.tar.gz", because there are some formal problems too. > > > > Kind regards > > Hugo Mildenberger > > > > Am 08.11.2011 16:58, schrieb Matt Owen: >> Hugo, thank you for your patch! Additionally, I will be reviewing the submitted bugs, and hopefully will have a fix in the very near future. >> >> Cheers - >> Matt Owen >> >> >> On Nov 6, 2011, at 3:53 PM, Hugo Mildenberger wrote: >> >>> Please use the new patch attached here instead of the previous one. The >>> previous version would likely break non-S4 packages used with >>> bootfn.default(). Besides this patch also fixes the plot problem. >>> Using set.seed(), regression testing was done for the previously working >>> non-bootstrapping code and resulted in no differences with >>> "logit.mixed", at least. >>> >>> Best >>> >>> >>> Am 05.11.2011 14:36, schrieb Hugo Mildenberger: >>>> Dear list members, >>>> >>>> in a recent thread of the R-sig-ME mailing list[1], Christoph Scherber >>>> communicated the following script which demonstrated a problem within >>>> Zelig's bootstrap code: >>>> >>>> require(Zelig) >>>> data(voteincome) >>>> >>>> z.out1 <- zelig(vote ~ education + age + female + tag(1 | state), >>>> data = voteincome, model = "logit.mixed") >>>> >>>> x <- setx(z.out1, education = 4) >>>> s.out1 <- sim(z.out1, x = x, bootstrap=TRUE) >>>> >>>> #Error in `*tmp*`$call : $ operator not defined for this S4 class >>>> >>>> >>>> This error stemmed from the implementation of bootfn.default() in >>>> Zelig/R/bootfn.default.R, which generally uses "object$call" instead of >>>> "object@call" (and therefore never worked with lme4) >>>> >>>> After having fixed this, the next problem became manifest when >>>> zelig::sim() was in bootstrapping mode. R could not find a summary >>>> method for the object passed to param.mer(). Why this happens, is still >>>> unclear to me. >>>> >>>> Anyway, when having fixed that bug along the lines of the >>>> non-bootstrapping branch of param.mer() by using selectMethod(), the >>>> next problem surfaced: >>>> >>>> Error in t.star[r, ] <- res[[r]] : >>>> incorrect number of subscripts on matrix >>>> >>>> After a while it became clear to me that the list returned by >>>> param.mer() was incompatible with the type boot::boot() expects (in >>>> /boot/R/bootfuns.q as of boot version 1.3-3): >>>> >>>> 193 } else lapply(seq_len(RR), fn) >>>> 194 t.star <- matrix(, RR, length(t0)) >>>> 195 for(r in seq_len(RR)) t.star[r, ] <- res[[r]] >>>> >>>> "t.star" is the name of the matrix boot::boot() uses to collect results >>>> from the statistic function. >>>> >>>> "fn" is the "statistic" parameter passed to boot::boot(), which in turn >>>> translates to bootfn.default() as defined in Zelig/R/bootfn.default.R. >>>> >>>> The function bootfn.default() in turn passes the result of param.mer(), >>>> which generally was "list(betas=betas, gammas=gammas, scale=scale)". >>>> >>>> length(t0) was three, because the list returned by param.mer() had three >>>> elements. The implicit type conversion of "t.star" (from matrix to list) >>>> R does after the first assignment to "t.star" then triggered the index >>>> error with the next assignment. >>>> >>>> A change of boot::boot() code in order to become compatible with Zelig >>>> appears to be out of question. >>>> >>>> However, in order to keep the default configuration working, param.mer() >>>> needs to return a list, because part of the result is a complete matrix >>>> (e.g."betas") then, not only a single row of a matrix. >>>> >>>> I therefore made the return type of param.mer(object,num,bootstrap,) >>>> dependent on "bootstrap", returning a list in the default case, and a >>>> one-dimensional array while bootstrapping. The related function qi.mer() >>>> in Zelig/R/qi.mixed.R is of course now able to handle both types. >>>> >>>> Since the column names appearing in the array returned by >>>> param.mer(..,bootstrap=T) are determined by the user controlled model >>>> specification, the name space separation a list provides implicitly had >>>> to be coded explicitly. >>>> >>>> The rudimentary regression testing I did for sim(,bootstrap=F) using >>>> patched code against the original version resulted in a quite close >>>> numerical agreement. I do attribute the differences (generally in the >>>> order of 1E-3) to a different RNG state, but I haven't examined this >>>> very closely. Generally, the effective code path of the default case >>>> should have remained unchanged. >>>> >>>> >>>> The attached patch is against the sources found in the Cran package >>>> Zelig_3.5.1.tar.gz >>>> >>>> >>>> Unrelated to this patch, plotting of simulated results generally does >>>> not work as intended: >>>> >>>>> plot(s.out1) >>>> Error in plot(density(qi), main = x$qi.name[[i]], xlab = xlab, ...): >>>> error in evaluating the argument 'x' in selecting a method for >>>> function 'plot': Error in density.default(qi) : argument 'x' >>>> must be numeric >>>> >>>> No plot generated for model logit.mixed >>>> >>>> >>>> Also, excluding the intercept using the ~ -1 notation is currently >>>> disregarded. >>>> >>>> >>>> Best regards >>>> >>>> Hugo Mildenberger >>>> >>>> >>>> [1]h ttps://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q4/006933.html >>> <Zelig_3.5.1.fix-boot-v2.patch> >> >> > <Zelig_3.5.1.fix-boot-v3.patch><R-CMD-check-Zelig_3.5.1.tar.gz.txt>