# [R] fastest R platform

Jason Liao jg_liao at yahoo.com
Mon Apr 9 04:46:37 CEST 2001

Hello, everyone! I picked up R several months ago and have adopted it
as my choice for statistical programming. Coming from a Java
background, I can honestly say that R is not only free, it is better
tha S-plus: the lexical scope in R makes it very simple to simulate
Java's object model. For this, I encourage everyone to read the artcle:

Robert Gentleman and Ross Ihaka (2000) "Lexical Scope and Statistical
Computing", Journal of Computational and Graphical Statistics, 9,
491-508.

I have coded a simulation program which runs at 1/5 the speed of
minimum usability on a Pentium 700 PC (windows me). I have spent quite
some time trying to optimize the program. So my question is: on which
platform does R run the fastest? I would be willing to invest on a

My program follows in case someone can help me improve it.

Thanks.

Jason

# Liao and Lipsitz (2001)
# A REML type estimate of variance components in generalized linear
mixed models
# This program runs on R. It does not run on S-plus.

rm(list=ls(all=TRUE))

###small and common functions##################################
glmm.sample.y <- function(x, z, b, D.u)
{
pp <- matrix(0, nrow = m, ncol = n);
D.u.half <- t( chol(D.u) );
zero <- numeric(n.random);

for(i in 1:m)
{
obj <- rmulti.norm(zero, D.u.half);
pp[i, ] <- as.vector( x[i, , ] %*% b + z[i, , ] %*% obj\$ran );
}
pp <- exp(pp)/( 1+exp(pp) );

y <- runif(m*n);
dim(y) <- c(m, n);
ifelse(y < pp, 1, 0);    #returning vector
}
################
rmulti.norm <- function(mean, var.half)
{
ran <- rnorm( length(mean) );
log.prob <-  -sum(ran*ran)/2;

ran <- mean + as.vector( var.half %*% ran );
return(ran, log.prob);
}
#################
my.outer <- function(x) #this is much faster than R's outer()
function
{
dim(x) <- c(length(x), 1);
x %*% t(x);
}

quadratic.form <- function(A, x) {  as.vector(x %*% A %*% x)  }

#######################################################################

REML.b.D.u <- function(b, D.u, y, x, z)
{
TT <- matrix(0, n.random, n.random);
for(i in 1:30)
{
obj <- sample.u(b, D.u, x, y, z);
b <- b - solve(obj\$Hessian, obj\$score);
D.u <- obj\$uu;

yy <- glmm.sample.y(x, z, b, D.u);  #this is the sampling stage
obj <- mle.b(b, D.u, yy, x, z);     #given D.u, solves for the
mle of b

TT <-  TT*(1-1/i) + obj\$uu.diff/i;
D.u <- D.u + TT;
print(i); print(date());
print("inside REML");
print(D.u);print("TT"); print(TT);
}

return(b, D.u);
}
mle.b <- function(b, D.u, y, x, z)  #b here is the initial value
{
for(i in 1:30)
{
obj <- sample.u(b, D.u, x, y, z);
b <- b - solve(obj\$Hessian, obj\$score);
if(i==1) uu.first <- obj\$uu;
}

uu.diff <- uu.first - obj\$uu;
return(b, uu=obj\$uu, uu.diff);
}

mle.b.D.u <- function(b, D.u, y, x, z)
{
for(i in 1:30)
{
obj <- sample.u(b, D.u, x, y, z);
b <- b - solve(obj\$Hessian, obj\$score);
D.u <- obj\$uu;
}
return(b, D.u);
}

sample.u <- function(b, D.u, x, y, z)
{
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)
{
obj <-  sample.u.one(D.u.inv, ada.part=as.vector(x[i, ,] %*% b),
x[i, ,], y[i,], z[i, ,]);
uu <- uu + obj\$uu;
score <- score + obj\$score;
Hessian <- Hessian + obj\$Hessian;
}

uu <- uu/m;
return(uu, score, Hessian);
}

z for one subject only
{
func <- one.log.like(D.u.inv, ada.part, y, z);  #function to be
maximized
obj <- optim( numeric(n.random), func, control=list(fnscale=-1) )

obj <- optim(obj\$par, func, method = c("BFGS"),
control=list(fnscale=-1), hessian = T )

mean <- obj\$par;
var.half <-  solve(-obj\$hessian)*1.2;
var.half <- t( chol(var.half) );

n.simu <- 1000  #sample size of importance sampling

sum.weight <- 0;
uu <- matrix(0, n.random, n.random);
pi <- numeric(n);
product <- numeric(n);

func <- one.score.Hessian(D.u.inv, ada.part, x, y, z);
for(i in 1:n.simu)
{
obj1 <- rmulti.norm(mean, var.half);
obj2 <- func(obj1\$ran);
weight <- exp(obj2\$log.likelihood - obj1\$log.prob);

uu <- uu + my.outer(obj1\$ran) * weight;
pi <- pi + obj2\$pi * weight;
product <- product + obj2\$pi*(1-obj2\$pi) * weight;

sum.weight <- sum.weight + weight;
}

uu <- uu/sum.weight;
pi <- pi/sum.weight;
product <- product/sum.weight;

score <- as.vector( (y - pi) %*% x );
Hessian <- -t(x) %*% diag(product) %*% x;

return(uu, score, Hessian);
}

one.log.like <- function(D.u.inv, ada.part, y, z)
{
function(u)
{
pp <- exp(  ada.part + as.vector( z %*% u )  );
pp <- pp/(1+pp);
-quadratic.form(D.u.inv, u)/2 + sum( y*log(pp) + (1 -
y)*log(1-pp) );
}
}

one.score.Hessian <- function(D.u.inv, ada.part, x, y, z)
{
function(u)
{
pi <- exp(  ada.part + as.vector( z %*% u )  );
pi <- pi/(1 + pi);
log.likelihood <- -quadratic.form(D.u.inv, u)/2 + sum(
y*log(pi) + (1-y)*log(1-pi) );

return(pi, log.likelihood);
}
}

main <- function()
{
b <- c(.5 , 1, -1, -.5);
if(n.fixed > 4) b <- c( b, numeric(n.fixed-4) );
D.u <- c(.5, 0, 0, .25);
D.u <- matrix(D.u, nrow=2);

x <- array(0, c(m, n, n.fixed) );
x[, , 1] <- 1;

for(j in 1:n) x[, j, 2] <- j-5;
x[1:trunc(m/2), , 3] <- 0;
x[trunc(m/2+1):m, , 3] <- 1;

x[, , 4] <- x[, , 2]*x[, , 3];
if( n.fixed > 4 ) x[, , 5:n.fixed] <- rnorm( m*n*(n.fixed-4) );

z <- x[, , 1:2];

for(i in 1:50)
{
y <- glmm.sample.y(x, z, b, D.u);

obj <- mle.b.D.u(b, D.u, y, x, z); print("mle"); print(obj\$D.u);

write(obj\$D.u, file = "simu1.dat", append=T);

obj <- REML.b.D.u(b, D.u, y, x, z); print("REML");
print(obj\$D.u);
write(obj\$D.u, file = "simu2.dat", append=T);
}
}
###################################################
#global variable m, n, n.fixed, n.random

n <- 10;
m <- 25;

n.fixed <- 4;
n.random <- 2;

main();
print(date());

=====
Jason G. Liao
Department of Biometry and Epidemiology
Medical University of South Carolina
135 Rutledge Ave., STE 1148, Charleston, SC 29425
phone (843) 876-1114, fax (843) 876-1126

http://www.geocities.com/jg_liao/index.html

__________________________________________________
Do You Yahoo!?
Get email at your own domain with Yahoo! Mail.
http://personal.mail.yahoo.com/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._