> From: Martin Maechler > > >>>>> "ReidH" == Huntsinger, Reid <[EMAIL PROTECTED]> > >>>>> 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> [EMAIL PROTECTED] > ReidH> [mailto:[EMAIL PROTECTED] On Behalf > ReidH> Of Eduardo Leoni Sent: Thursday, March 03, 2005 5:16 > ReidH> PM To: [email protected] 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> [email protected] 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> [email protected] 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 > > ______________________________________________ > [email protected] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > > > ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
