RE: [R] regression on a matrix

From: Huntsinger, Reid <reid_huntsinger_at_merck.com>
Date: Fri 04 Mar 2005 - 09:24:22 EST


You might use lsfit instead and just do the whole Y matrix at once. That saves all the recalculation of things involving only X.

Reid Huntsinger

-----Original Message-----
From: r-help-bounces@stat.math.ethz.ch
[mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Eduardo Leoni Sent: Thursday, March 03, 2005 5:16 PM
To: r-help@stat.math.ethz.ch
Subject: [R] regression on a matrix

Hi -

I am doing a monte carlo experiment that requires to do a linear regression of a matrix of vectors of dependent variables on a fixed set of covariates (one regression per vector). I am wondering if anyone has any idea of how to speed up the computations in R. The code follows:

#regression function
#Linear regression code
qreg <- function(y,x) {
  X=cbind(1,x)
  m<-lm.fit(y=y,x=X)
  p<-m$rank

  r <- m$residuals
  n <- length(r)
  rss <- sum(r^2)
  resvar <- rss/(n - p)   

  Qr <- m$qr
  p1 <- 1:p
  R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])   se <- sqrt(diag(R) * resvar)
  b <- m$coefficients
  return(c(b[2],se[2]))
}

#simulate

a <- c(1,.63,.63,1)
a <- matrix(a,2,2)
c <- chol(a)
C <- 0.7

theta <- 0.8
sims <- 1000
n<-20
u <- rnorm(n,0,sqrt(1-C))
w <- rgamma(n,C/theta,1/theta)
e <- rnorm(n,0,sqrt(w))
  

x1 <- rnorm(n)

x <- x1*c[2,2]+c[1,2]*w
v <- e+u
y <- 1+x+v
w <- rgamma(n,C/theta,1/theta)

#create matrix of dep variable
newdep <- matrix(rnorm(length(y)*sims,y,sqrt(w)),c(length(y),sims))

monte <- apply(newdep,2,qreg,x=x)



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Mar 04 09:49:09 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:30:40 EST