diff -ur a/Zelig/R/bootfn.default.R b/Zelig/R/bootfn.default.R
--- a/Zelig/R/bootfn.default.R	2011-10-04 12:58:06.000000000 +0200
+++ b/Zelig/R/bootfn.default.R	2011-11-04 15:41:42.000000000 +0100
@@ -1,12 +1,12 @@
 bootfn.default <- function(data, i, object) {
   d <- data[i,]
-  object$call$data <- d
+  object@call$data <- d
   l <- length(param(object, bootstrap = TRUE))
-  fit <- eval(object$call, sys.parent())
+  fit <- eval(object@call, parent.frame())
   l1 <- length(param(fit, bootstrap = TRUE))
   while (l > l1) {
-    object$call$data <- data[sample(nrow(data), replace=TRUE),]
-    fit <- eval(object$call, sys.parent())
+    object@call$data <- data[sample(nrow(data), replace=TRUE),]
+    fit <- eval(object@call, parent.frame())
     l1 <- length(param(fit, bootstrap = TRUE))
   }
   param(fit, bootstrap = TRUE)
diff -ur a/Zelig/R/param.mixed.R b/Zelig/R/param.mixed.R
--- a/Zelig/R/param.mixed.R	2011-10-04 12:58:06.000000000 +0200
+++ b/Zelig/R/param.mixed.R	2011-11-04 21:07:59.000000000 +0100
@@ -22,21 +22,41 @@
 			J <- dim(vars.m)[1]
 			gammas[[m]] <- mvrnorm(num, rep(0, J), V.beta)
 		}
-	}
-	else {
-		object <- summary(object)
+                # betaF is a matrix, so there is a real need for a (nested) list
+		names(gammas) <- names(vars)
+		betas <- betasF
+		scale <- object@sigma
+		return(list(betas=betas, gammas=gammas, scale=scale))
+	} else {
+		object <- selectMethod("summary", "mer")(object)
 		# set fixed effects
-		betasF <- fixed
+		betasF <- fixef(object)
 		# sample random effects
 		for (m in 1:n.G){
 			V.beta <- VarCorr(object)[[m]]
 			gammas[[m]]<- mvrnorm(1, 0, V.beta)
 		}
+
+		# boot::boot() collects results in a matrix (named "t.star").
+		# Returning a list will therefore create havoc with boot(), 
+		# also because boot() does only reserve space for length(list),
+		# which would result in three for a list like 
+		# list(betas=betas, gammas=gammas, scale=scale)
+		# 
+		# One solution to this situation is to return a vector, and 
+		# pass the name attributes via names(). However, the beta and 
+		# gamma namespaces may intersect, and certainly user models may
+		# also use 'scale' within fixed or random effects. Hence prepend 
+		# a local 'namespace' to all variable names in order to emulate 
+		# a list return type here.
+		scale          <- object@sigma
+		betas          <- betasF
+
+		names(scale)   <- sprintf('a::%s','')
+		names(betas)   <- sprintf('b::%s',names(betasF))
+		names(gammas)  <- sprintf('g::%s',names(vars))
+
+		return(c(betas,gammas,scale))
 	}
-	
-	names(gammas) <- names(vars)
-	betas <- betasF
-	scale <- object@sigma
-	
-	list(betas=betas, gammas=gammas, scale=scale)
-}
\ Kein Zeilenumbruch am Dateiende.
+	stopifnot(F)
+}
diff -ur a/Zelig/R/qi.mixed.R b/Zelig/R/qi.mixed.R
--- a/Zelig/R/qi.mixed.R	2011-10-04 12:58:06.000000000 +0200
+++ b/Zelig/R/qi.mixed.R	2011-11-04 20:28:48.000000000 +0100
@@ -15,11 +15,45 @@
   D <- eval(object@call$data, envir = parent.frame())
   fml <- eval(object@call$formula)
   parsefml <- .getRandAndFixedTerms(fml)
+
+  #
+  # boot::boot() seems to have changed from a list 
+  # result type into an array. Therefore the original
+  # code below will not work  with bootstrap=T:
+  #
+  # Hence, while bootstraping, the related function 
+  # param.mer() (in param.mixed.R) will return an array
+  # prepending "a::","b::","g::" to the names of scale 
+  # (=mer@sigma), betas(=fixef(mer)) and gammas, respectively.
+
+  #print("qi.mixed.R.simpar=");print(simpar)
+
+  if ( is.list(simpar) )  {
   
-  betas <- simpar[["betas"]]
-  gammas <- simpar[["gammas"]]
-  alpha <- simpar[["scale"]]
-  
+	  betas  <- simpar[["betas"]]
+	  gammas <- simpar[["gammas"]]
+	  alpha  <- simpar[["scale"]]
+
+  } else {
+
+	  alpha.idx  <- grep('a::',colnames(simpar))
+	  betas.idx  <- grep('b::',colnames(simpar))
+	  gammas.idx <- grep('g::',colnames(simpar))
+
+	  alpha  <- as.matrix(simpar[,alpha.idx])
+	  betas  <- as.matrix(simpar[,betas.idx])
+	  gammas <- as.matrix(simpar[,gammas.idx])
+
+	  colnames(alpha)  <- gsub("a::","",colnames(simpar)[alpha.idx])
+	  colnames(betas)  <- gsub("b::","",colnames(simpar)[betas.idx])
+	  colnames(gammas) <- gsub("g::","",colnames(simpar)[gammas.idx])
+
+	  alpha <- as.vector(unique(alpha))
+  }
+
+  #print("qi.mixed.R.alpha=");print(alpha)
+  #print("qi.mixed.R.betas=");print(betas)
+  #print("qi.mixed.R.gammas=");print(gammas)
   
   fTermsNames <- colnames(model.matrix(parsefml$fixed, data = D))
                                         #cat(fTermsNames, "\n")
@@ -65,10 +99,18 @@
     
     eta <- betas %*% t(as.matrix(fTerms))
     mu <- eta
-    ## For predicted values, add in random effects draws
-    for (i in 1:length(rTerms)){
-      mu <- mu + gammas[[names(rTerms[i])]] %*% t(as.matrix(rTerms[[i]]))
+
+    ## For predicted values, add in random effects draws, and
+    ## account for the now Janus-faced return type of param.mer()
+    ## This could be done without code duplication, but for now it works at least.
+    if (is.list(simpar)) { 
+         for (i in 1:length(rTerms))
+              mu <- mu + gammas[[names(rTerms[i])]] %*% t(as.matrix(rTerms[[i]]))
+    } else { 
+         for (i in 1:length(rTerms))
+              mu <- mu + gammas[,names(rTerms[i])] %*% t(as.matrix(rTerms[[i]]))
     }
+         
     if (is.null(family)){
       ev <- eta
       n <- length(mu[,1])