R help - Mar 2005 - 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)

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 at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Eduardo Leoni
Sent: Thursday, March 03, 2005 5:16 PM
To: r-help at 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 at 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

> From: Martin Maechler
> 
> >>>>> "ReidH" == Huntsinger, Reid
<reid_huntsinger at merck.com>
> >>>>>     on Thu, 3 Mar 2005 17:24:22 -0500 writes:
> 
>     ReidH> You might use lsfit instead and just do the whole Y
>     ReidH> matrix at once. That saves all the recalculation of
>     ReidH> things involving only X.
> 
> yes,  but in these cases, we have been recommending
> lm.fit() instead -- just so you use the identical internal
> numeric code as lm() and still have the `benefit' of not having
> to re-build the design matrix X .
Yes, but ?lm.fit says y has to be a vector, while ?lsfit says y can be a
matrix.  If lm.fit can handle y as a matrix, its help page should be
updated.

Andy

 
> Martin Maechler, ETH Zurich
> 
> 
>     ReidH> Reid Huntsinger
> 
>     ReidH> -----Original Message----- From:
>     ReidH> r-help-bounces at stat.math.ethz.ch
>     ReidH> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf
>     ReidH> Of Eduardo Leoni Sent: Thursday, March 03, 2005 5:16
>     ReidH> PM To: r-help at stat.math.ethz.ch Subject: [R]
>     ReidH> regression on a matrix
> 
> 
>     ReidH> Hi -
> 
>     ReidH> I am doing a monte carlo experiment that requires to
>     ReidH> do a linear regression of a matrix of vectors of
>     ReidH> dependent variables on a fixed set of covariates (one
>     ReidH> regression per vector). I am wondering if anyone has
>     ReidH> any idea of how to speed up the computations in
>     ReidH> R. The code follows:
> 
>     ReidH> #regression function #Linear regression code qreg <-
>     ReidH> function(y,x) { X=cbind(1,x) m<-lm.fit(y=y,x=X)
>     ReidH> p<-m$rank
> 
>     ReidH>   r <- m$residuals n <- length(r) rss <- sum(r^2)
>     ReidH> resvar <- rss/(n - p)
>   
>     ReidH>   Qr <- m$qr p1 <- 1:p R <- chol2inv(Qr$qr[p1, p1,
>     ReidH> drop = FALSE]) se <- sqrt(diag(R) * resvar) b <-
>     ReidH> m$coefficients return(c(b[2],se[2])) }
> 
> 
>     ReidH> #simulate a <- c(1,.63,.63,1) a <- matrix(a,2,2) c
<-
>     ReidH> chol(a) C <- 0.7 theta <- 0.8 sims <- 1000 n<-20
> 
>     ReidH> u <- rnorm(n,0,sqrt(1-C)) w <-
>     ReidH> rgamma(n,C/theta,1/theta) e <- rnorm(n,0,sqrt(w))
>   
>     ReidH> x1 <- rnorm(n) x <- x1*c[2,2]+c[1,2]*w v <- e+u y
<-
>     ReidH> 1+x+v w <- rgamma(n,C/theta,1/theta)
> 
>     ReidH> #create matrix of dep variable newdep <-
>     ReidH> matrix(rnorm(length(y)*sims,y,sqrt(w)),c(length(y),sims))
> 
> 
>     ReidH> monte <- apply(newdep,2,qreg,x=x)
> 
>     ReidH> ______________________________________________
>     ReidH> R-help at stat.math.ethz.ch mailing list
>     ReidH> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
>     ReidH> do read the posting guide!
>     ReidH> http://www.R-project.org/posting-guide.html
> 
>     ReidH> ______________________________________________
>     ReidH> R-help at stat.math.ethz.ch mailing list
>     ReidH> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
>     ReidH> do read the posting guide!
>     ReidH> http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at 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
> 
> 
>