RE: [R] regression on a matrix

From: Liaw, Andy <andy_liaw_at_merck.com>
Date: Fri 04 Mar 2005 - 22:33:49 EST


> From: Martin Maechler
>
> >>>>> "ReidH" == Huntsinger, Reid <reid_huntsinger@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@stat.math.ethz.ch
> ReidH> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf
> ReidH> Of Eduardo Leoni Sent: Thursday, March 03, 2005 5:16
> ReidH> PM To: r-help@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@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@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@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 23:52:43 2005

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