I care a lot about R's speed. So I decided to give REvolution's R (http://revolution-computing.com/) a try, which bills itself as an optimized R. Note that I used the free version.
My machine is a Intel core 2 duo under Windows XP professional. The code I run is in the end of this post. First, the regular R 1.9. It takes 2 minutes and 6 seconds, CPU usage 50% Next, REvolution's R. It takes 2 minutes and 10 seconds, CPU usage 100%. In other words, REvolution's R consumes double the CPU with slightly less speed. The above has been replicated a few times (as a statistician of course). Anyone has any insight on this? Anyway, my high hope was dashed. rm(list=ls(all=TRUE)) library(MASS); ###small and common functions################################## glmm.sample.y <- function(b, D.u, x, z) { pp <- matrix(0, nrow = n, ncol = m); zero <- numeric(n.random); ran <- t( mvrnorm(m, zero, D.u) ); for(j in 1:m) pp[, j] <- as.vector( x[, , j] %*% b + z[, , j] %*% ran[ ,j] ); pp <- exp(pp)/( 1+exp(pp) ); y <- runif(m*n); ifelse(y<pp, 1, 0); } ################# quadratic.form <- function(A, x) { return( as.vector(x %*% A %*% x) ) } quadratic.form2 <- function(A, x) { x <- chol(A) %*% x; colSums(x*x) } ####################################################################### REML.b.D.u <- function(b, D.u, x, y, z, n.iter) { u.mean.initial <- array( 0, c(n.random, m) ); TT <- matrix(0, n.random, n.random); for(i in 1:n.iter) { TT <- TT*(1-1/i) + bias(b, D.u, x, z)/i; obj <- sample.u(b, D.u, x, y, z, u.mean.initial); b <- b - solve(obj$Hessian, obj$score); D.u <- obj$uu + TT; u.mean.initial <- obj$u.mean.initial; print(i); print(date()); print("inside REML"); print(TT); if(i==n.iter) write(TT, file="c:/liao/reml/simu50_8.dat", append=T); print(D.u); } list(b=b, D.u=D.u); } bias <- function(b, D.u, x, z) { yy <- glmm.sample.y(b, D.u, x, z); #this is the sampling stage obj1 <- mle(b, D.u, x, yy, z, F, F, n.iter=50); obj2 <- mle(b, D.u, x, yy, z, T, F, n.iter=50); obj1$uu - obj2$uu } ################################################## mle <- function(b, D.u, x, y, z, indx1, indx2, n.iter) { u.mean.initial <- array( 0, c(n.random, m) ); for(i in 1:n.iter) { obj <- sample.u(b, D.u, x, y, z, u.mean.initial); if(indx1) b <- b - solve(obj$Hessian, obj$score); if(indx2) D.u <- obj$uu; u.mean.initial <- obj$u.mean.initial; } list(b=b, D.u=D.u, uu=obj$uu, u.mean.initial=u.mean.initial); } ############################################## sample.u <- function(b, D.u, x, y, z, u.mean.initial) { D.u.inv <- solve(D.u); uu <- matrix(0, n.random, n.random); score <- numeric(n.fixed); Hessian <- matrix(0, n.fixed, n.fixed); for(i in 1:m) { ada.part <- as.vector(x[, , i] %*% b); obj <- sample.u.one( D.u.inv, ada.part, x[, , i], y[, i], z[, ,i], u.mean.initial[, i] ); uu <- uu + obj$uu; score <- score + obj$score; Hessian <- Hessian + obj$Hessian; u.mean.initial[, i] <- obj$u.mean.initial; } uu <- uu/m; list(uu=uu, score=score, Hessian=Hessian, u.mean.initial=u.mean.initial); } ########################################## sample.u.one <- function(D.u.inv, ada.part, x, y, z, u.mean.initial) #ada.part, x, y, z for one subject only { fn <- log.like(y, z, D.u.inv, ada.part); obj <- optim(u.mean.initial, fn, control=list(fnscale=-1), hessian = T) u.mean.initial <- obj$par; var.root <- solve(-obj$hessian); var.root <- t( chol(var.root) ); u <- matrix( rnorm(n.random*n.simu), nrow=n.random, ncol=n.simu ); log.prob1 <- -colSums(u*u)/2; u <- u.mean.initial + var.root %*% u; ada <- exp( ada.part + z %*% u ); ada <- ada/(1+ada); #it is probability now log.prob2 <- colSums( y*log(ada) + (1-y)*log(1-ada) ); log.prob2 <- -quadratic.form2(D.u.inv, u)/2 + log.prob2; weight <- exp(log.prob2 - log.prob1); weight <- weight/sum(weight); ada <- t(ada); pi <- colSums(ada*weight); product <- colSums( ada*(1-ada)*weight ); score <- as.vector( (y - pi) %*% x ); Hessian <- -t(x) %*% ( rep(product, n.fixed)*x ); obj <- cov.wt( t(u), weight, center=T ); uu <- obj$cov + obj$center %*% t(obj$center); list(uu=uu, score=score, Hessian=Hessian, u.mean.initial=obj$center); } ######################################### log.like <- function(y, z, D.u.inv, ada.part) { function(u) { ada <- exp( ada.part + z %*% u ); ada <- ada/(1+ada); sum( y*log(ada) + (1-y)*log(1-ada) ) -quadratic.form(D.u.inv, u)/2; } } main <- function() { b <- c(.5, .5 , -1, -.5); b <- c(.5, .5 , -1, -.5, 0, 0, 0, 0); D.u <- diag( c(.5, .5) ); x <- array(0, c(n, n.fixed, m) ); x[ ,1, ] <- 1; for(i in 1:n) x[i, 2, ] <- (i-5)/4; x[, 3, 1:trunc(m/2)] <- 0; x[, 3, trunc(m/2+1):m] <- 1; x[, 4, ] <- x[, 2, ]*x[, 3, ]; x[1, 5:n.fixed, ] <- rnorm(4*m); for(i in 2:n) x[i, 5:n.fixed,] <- x[1, 5:n.fixed, ]; z <- x[, 1:2, ]; for(i in 1:200) { print(i); print(date()); y <- glmm.sample.y(b, D.u, x, z); obj <- mle(b, D.u, x, y, z, F, T, n.iter=50); print("estimate with b known") print(obj$D.u); obj <- mle(b, D.u, x, y, z, T, T, n.iter=50); print("profile MLE"); print(obj$D.u); obj <- REML.b.D.u(b, D.u, x, y, z, n.iter=50); print("REML type estimate") print(obj$D.u); } } ################################################### n <- 10; #number of observation for each cluster m <- 50; #number of clusters n.fixed <- 8; n.random <- 2; n.simu <- 2000; main(); [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.