## code taken from Implementation.Rnw in the lmer package, ## showing how to construct a likelihood profile for ## the random effects. This will also show how we ## can find the deviance for 0 (or small) random effects, ## which we can compare with the glm result [which it ## should equal if the computation is being done similarly] ## and with the glmer result) ## *** N.B. THE FOLLOWING SET OF CALLS ACTUALLY MODIFIES THE INTERNAL ## STRUCTURE OF THE ORIGINAL MODEL! MAKE A COPY IF YOU WANT TO ## PRESERVE THE ORIGINAL *** require(lme4) varprof <- function(mm,lower=0,upper=20,n=101) { sg <- seq(lower, upper, len = n) orig.sd <- attr(VarCorr(mm)[[1]],"stddev") dev <- mm@deviance nc <- length(dev) nms <- names(dev) vals <- matrix(0, nrow = length(sg), ncol = nc, dimnames = list(NULL, nms)) update_dev <- function(sd) { .Call("mer_ST_setPars", mm, sd, PACKAGE = "lme4") .Call("mer_update_L", mm, PACKAGE = "lme4") res <- try(.Call("mer_update_RX", mm, PACKAGE = "lme4"), silent = TRUE) if (inherits(res, "try-error")) { val <- NA } else { .Call("mer_update_ranef", mm, PACKAGE = "lme4") .Call("mer_update_dev", mm, PACKAGE = "lme4") ## added for glmmML val <- mm@deviance } val } for (i in seq(along = sg)) { vals[i,] <- update_dev(sg[i]) } update_dev(orig.sd) ## hack! to restore original sd data.frame(sd=sg,vals) } fixprof <- function(mm,n=31,verbose=FALSE) { dev <- mm@deviance nms <- c("parval",names(dev)) npar <- length(fixef(mm)) nc <- length(nms) parnames <- names(fixef(mm)) vals <- list() ## retrieve est sds, use these to set bounds for profiles sds <- sqrt(diag(vcov(mm))) lowvec <- fixef(mm)-3*sds hivec <- fixef(mm)+3*sds mmat <- mm@X ## model matrix colnames(mmat) <- parnames has.intercept <- as.logical(attr(attr(mm@frame,"terms"),"intercept")) if (!has.intercept) warning("not working with no-intercept models yet") mmat1 <- if (has.intercept) mmat[,-1] else mmat ## trying to find a way to combine the response variable, ## model frame, and random effects ## random factors formula <- mm@call$formula reffs <- as.character(sapply(lme4:::findbars(formula),"[[",3)) respvar <- as.character(attr(mm@frame,"terms")[[2]]) x0 <- data.frame(as.data.frame(c(mm@frame))[c(respvar,reffs)], as.data.frame(mmat1)) for (param in 1:npar) { ## set up results matrix vals[[param]] <- matrix(0, nrow = n, ncol = nc, dimnames = list(NULL, nms)) ## eff <- mmat[,param] ## n.b. there should be a better way to manipulate ## the formula, either by operating on the model matrix ## directly and passing it to glmer, or by dropping terms? ## want to keep random effects and all but one fixed effect ## (was previously trying to drop a column of the model ## matrix and use predation~., but this fails to include ## the random effect) f <- paste(respvar,"~", paste(colnames(mmat)[-c(1,param)],collapse="+")) if (param==1) f <- paste(f,"-1") ## explicitly drop intercept term f <- paste(f,"+ (", paste(as.character(lme4:::findbars(formula)), collapse=")+("), ")",sep="") ## not really necessary: check to see that ## we recover original model results ## if (FALSE) { ##g0 <- glmer(as.formula(f), ## family=binomial,nAGQ=8,data=x0, ## offset=fixef(mm)[param]*eff) ## } sg <- seq(lowvec[param], hivec[param], length = n) for (i in 1:n) { ##tmpfit <- glmer(as.formula(f), ## family=binomial,nAGQ=8,data=x0, ## offset=sg[i]*eff) ## ## could probably do this more efficiently/elegantly ## if I knew the magic mer_update formula and ## manipulated the model frame etc. directly ... ## alternative approach (more general?) tmpcall <- mm@call tmpcall$formula <- as.formula(f) tmpcall$offset <- sg[i]*eff tmpcall$data <- x0 tmpfit <- eval(tmpcall) vals[[param]][i,] <- c(sg[i],tmpfit@deviance) ## browser() if (verbose) cat(parnames[param],param,i,sg[i], tmpfit@deviance["ML"],"\n") } ## loop within parameter vals[[param]] <- as.data.frame(vals[[param]]) } ## loop over parameters names(vals) <- parnames vals } dropparm <- function(fit,param) { parnames <- names(fixef(fit)) mmat <- fit@X ## model matrix colnames(mmat) <- parnames has.intercept <- as.logical(attr(attr(fit@frame,"terms"),"intercept")) if (!has.intercept) warning("not working with no-intercept models yet") mmat1 <- if (has.intercept) mmat[,-1] else mmat formula <- fit@call$formula reffs <- as.character(sapply(lme4:::findbars(formula),"[[",3)) respvar <- as.character(attr(fit@frame,"terms")[[2]]) f <- paste(respvar,"~", paste(colnames(mmat)[-c(1,param)],collapse="+")) if (param==1) f <- paste(f,"-1") ## explicitly drop intercept term f <- paste(f,"+ (", paste(as.character(lme4:::findbars(formula)), collapse=")+("), ")",sep="") x0 <- data.frame(as.data.frame(c(fit@frame))[c(respvar,reffs)], as.data.frame(mmat1)) tmpcall <- fit@call tmpcall$formula <- as.formula(f) tmpcall$data <- x0 eval(tmpcall) } ## attempt at a reasonably general simulation function simvals <- function(N=c(block=10,x=10,rep=10), sd=c(block=0.1), pars=c(Intercept=1,x=1), xvals=1:10, family, seed=NULL) { ## how to specify nesting etc.? (1|block) + error [(1|rep:x:block)] ## should look as close as possible to lme4 formula specification ## need -1 to eliminate individual-level error for GLMMs r_eff <- list(block=) ## mapply(rnorm,N,sd,MoreArgs=list(mean=0)) ## generate random effects ## specify covariates -- linear, runif, rnorm, ... ? ## parameters? if (is.function(xvals)) { } zz <- expand.grid(block=factor(1:N["block"]),x=xvals,rep=1:N["rep"]) mz <- model.matrix(~x,data=zz) zz <- within(zz, { eta0 <- with(pars,a+b*x+b_eff[block]) eta1 <- rnorm(length(eta0),eta0,0.1) y <- rpois(length(eta0),exp(eta0)) }) zz2 <- within(zz, { eta0 <- a+b*x+b_eff[block] eta1 <- rnorm(length(eta0),eta0,0.1) y <- rbinom(length(eta0),prob=plogis(eta0),size=20) }) } rqpois <- function(n,lambda,phi,roundvals=FALSE) { sc <- phi sh <- lambda/phi r <- rgamma(n,shape=sh,scale=sc) if (roundvals) round(r) else r } rqpois2 <- function(n,lambda,phi,min.k=1e-3,roundvals=FALSE) { mu <- lambda k <- mu/(phi*mu-1) if (any(k0) { ## GLMM ## n.b. DON'T scale random-effects ## if sigma!=1, it applies to the "quasi"- part of the model etasim <- etasim.fix+etasim.reff family <- object@call$family if(is.symbol(family)) family <- as.character(family) if(is.character(family)) family <- get(family, mode = "function", envir = parent.frame(2)) if(is.function(family)) family <- family() if(is.null(family$family)) stop("'family' not recognized") musim <- family$linkinv(etasim) n <- length(musim) ## FIXME: or could be dims["n"]? vsim <- switch(family$family, poisson=rpois(n,lambda=musim), quasipoisson=rqpois(n,musim,sigma), ## quasipois hack, binomial={ resp <- model.response(object@frame); if (!is.matrix(resp)) { binom1 <- TRUE rbinom(n,prob=musim,size=1) } else { binom1 <- FALSE sizes <- rowSums(resp) Y <- rbinom(n, size = sizes, prob = musim) YY <- cbind(Y, sizes - Y) colnames(YY) <- colnames(resp) yy <- split(as.data.frame(YY), rep(1:nsim,each=length(sizes))) yy <- as.data.frame(yy) } }, stop("simulation not implemented for family", family$family)) if (!(family$family=="binomial" && !binom1)) { vsim <- matrix(vsim,nc=nsim) } return(drop(vsim)) ## note: drop matrix structure } stop("simulate method for NLMMs not yet implemented") } setMethod("refit", signature(object = "mer", newresp = "matrix"), function(object, newresp, ...) { ## newresp <- as.double(newresp[!is.na(newresp)]) wts <- rowSums(newresp) newresp <- newresp[,1]/wts stopifnot(length(newresp) == object@dims["n"]) object@y <- newresp object@pWt <- wts lme4:::mer_finalize(object) }) fixnullsim <- function(model1,model0,drop, stats="deviance",n=1000) { if (is.matrix(model.response(model1))) stop("can't handle two-column binomial responses yet") ## drop: name of parameter(s)? to drop ## model1: full model ## model0: reduced model (may be constructed automatically) if (!missing(model0) && !missing(drop)) stop("specify only one of reduced model or parameter to drop") if (missing(model0) && missing(drop)) stop("must specify either reduced model or parameter to drop") ## which (character): name of factor to drop in null sim if (missing(model0)) { npar <- length(fixef(mm)) parnames <- names(fixef(mm)) if (!drop %in% parnames) stop("parameter ",drop," not in parameter names") w <- which(parnames==drop) mmat <- model1@X ## model matrix colnames(mmat) <- parnames has.intercept <- as.logical(attr(attr(model1@frame,"terms"),"intercept")) if (!has.intercept) warning("not working with no-intercept models yet") mmat1 <- if (has.intercept) mmat[,-1] else mmat ## trying to find a way to combine the response variable, ## model frame, and random effects ## random factors formula <- model1@call$formula reffs <- as.character(sapply(lme4:::findbars(formula),"[[",3)) respvar <- as.character(attr(model1@frame,"terms")[[2]]) x0 <- data.frame(as.data.frame(c(model1@frame))[c(respvar,reffs)], as.data.frame(mmat1)) f <- paste(respvar,"~", paste(colnames(mmat)[-c(1,w)],collapse="+")) if (w==1) f <- paste(f,"-1") ## explicitly drop intercept term f <- paste(f,"+ (", paste(as.character(lme4:::findbars(formula)), collapse=")+("), ")",sep="") call0 <- model1@call call0$formula <- as.formula(f) call0$data <- x0 model0 <- eval(call0) } statm <- match(stats,c("deviance")) if (any(is.na(statm))) stop("unknown statistic(s):",paste(stats[is.na(statm)],",")) vals <- matrix(nrow=n,ncol=length(stats), dimnames=list(NULL,stats)) for (i in 1:n) { x2 <- my.mer.sim(model0) ## FIXME: need more logic for two-column binomial -- ## need to cbind x2 with size m0fit <- refit(model0,x2) m1fit <- refit(model1,x2) ## calculate & save statistics vals[i,] <- deviance(m0fit)-deviance(m1fit) ## progress bar? } } waldF <- function(object,variable,den.df) { if (missing(den.df)) stop("must specify denominator degrees of freedom") f <- object@frame terms <- attr(f,"terms") if (is.character(variable)) { cvar <- variable w <- match(variable,attr(terms,"term.labels")) variable <- which(attr(terms,"order") %in% w) } else cvar <- paste(as.character(variable),collapse=",") cc <- fixef(object)[variable] nn <- length(variable) vv <- vcov(object)[variable,variable] ## better way to do this?? fstat <- as.numeric(cc %*% solve(vv) %*% cc) r <- list(method="Wald F test", data.name=paste(deparse(substitute(object)), ": vars ",cvar,sep=""), statistic=c("Wald F"=fstat), p.value=pf(fstat,nn,den.df,lower.tail=FALSE)) class(r) <- "htest" r } ## wrap glmer to make warnings into errors lmer2 <- function(...,action=c("stop","code")) { action <- match.args(action) op <- options(warn=2) on.exit(options(op)) x <- try(lmer(...),silent=TRUE) if (inherits(x,"try-error")) { if (action=="stop") { stop(x,unclass(x)) } else if (action=="code") { stop("stub: code response not implemented yet") } } ## otherwise OK -- fall through x } glmer2 <- lmer2 ## wrap refit to make warnings into errors refit2 <- function(...) { op <- options(warn=2) on.exit(options(op)) x <- try(refit(...),silent=TRUE) if (inherits(x,"try-error")) { stop(x,unclass(x)) } x } Z.to.t <- function(ztab,df) { ztab[,4] <- 2*pt(abs(ztab[,3]),df=df,lower.tail=FALSE) colnames(ztab)[3:4] <- c("t value","Pr(>|t|)") ztab } ## The lm method now includes "glm" ones : ## SUPERSEDED in R 2.9.x and later if (FALSE) { simulate.lm <- function(object, nsim = 1, seed = NULL, type=c("response","data"), ...) { type <- match.arg(type) if(!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) runif(1) # initialize the RNG if necessary if(is.null(seed)) RNGstate <- get(".Random.seed", envir = .GlobalEnv) else { R.seed <- get(".Random.seed", envir = .GlobalEnv) set.seed(seed) RNGstate <- structure(seed, kind = as.list(RNGkind())) on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv)) } ftd <- fitted(object) # == napredict(*, object$fitted) ntot <- length(ftd) * nsim fam <- if(inherits(object, "glm")) object$family$family else "gaussian" val <- switch(fam, "gaussian" = { vars <- deviance(object)/ df.residual(object) if (!is.null(object$weights)) vars <- vars/object$weights ftd + rnorm(ntot, sd = sqrt(vars)) } , "poisson" = rpois(ntot, ftd) , "binomial" = { wts <- object$prior.weights if (any(wts %% 1 != 0)) stop("cannot simulate from non-integer prior.weights") r <- rbinom(ntot, size = wts, prob = ftd) if (type=="response") { r/wts } else { resp <- model.response(object$model) if (is.matrix(resp)) { cbind(r,wts-r) } else { if (is.factor(resp)) { levels(resp)[1+r] } else r } } }, "Gamma" = { if(is.null(tryCatch(loadNamespace("MASS"), error = function(e) NULL))) stop("Need 'MASS' package for 'Gamma' family") shape <- MASS::gamma.shape(object)$alpha rgamma(ntot, shape = shape, rate = shape/ftd) }, "inverse.gaussian" = { if(is.null(tryCatch(loadNamespace("SuppDists"), error = function(e) NULL))) stop("Need CRAN package 'SuppDists' for 'inverse.gaussian' family") lambda <- 1/summary(object)$dispersion SuppDists::rinvGauss(ntot, nu = ftd, lambda = lambda) }, ## MASS::glm.nb fits if (substr(fam,1,17)=="Negative Binomial") { size <- object$theta rnbinom(ntot, mu=ftd, size=size) } else { stop("family '", fam, "' not yet implemented") }) if (!(fam=="binomial" && type=="data")) { warning("FIXME: haven't decided how to return 2-column binomial data yet") ans <- as.data.frame(matrix(val, ncol = 2*nsim)) } else { ans <- as.data.frame(matrix(val, ncol = nsim)) } attr(ans, "seed") <- RNGstate ans } } if (FALSE) { ## examples: modified from ?simulate set.seed(1001) x <- 1:5 mod1 <- lm(c(1:3,7,6) ~ x) S1 <- simulate(mod1, nsim = 4) ## repeat the simulation: .Random.seed <- attr(S1, "seed") identical(S1, simulate(mod1, nsim = 4)) mod2 <- lm(c(1:3,7,6) ~ x, weights=c(1,2,2,1,1)) S2 <- simulate(mod2, nsim = 4) x <- 1:10 n <- 10 y <- rbinom(length(x),prob=plogis((x-5)/2),size=n) y2 <- c("a","b")[1+rbinom(length(x),prob=plogis((x-5)/2),size=1)] mod1 <- glm(cbind(y,n-y) ~ x,family=binomial) mod2 <- glm(factor(y2) ~ x,family=binomial) S1 <- simulate(mod1, nsim = 4) S1B <- simulate(mod2, nsim = 4) ## repeat the simulation: .Random.seed <- attr(S1, "seed") identical(S1, simulate(mod1, nsim = 4)) S2 <- simulate(mod1, nsim = 200, seed = 101) rowMeans(S2)/10 # after correcting for binomial sample size, should be about fitted(mod1) plot(rowMeans(S2)/10) lines(fitted(mod1)) ## repeat identically: (sseed <- attr(S2, "seed")) # seed; RNGkind as attribute stopifnot(identical(S2, simulate(mod1, nsim = 200, seed = sseed))) x <- 2:20 y3 <- rpois(length(x),exp((x-5)/5)) mod3 <- glm(y3 ~ x,family=poisson) S3 <- simulate(mod3, nsim = 4) if (require(MASS)) { mod4 <- MASS:::glm.nb(y3~x) S4 <- simulate(mod4, nsim = 4) } }