Naoki Takebayashi
2001-Jun-21 17:24 UTC
[R] problem with generalized least square, and lm.gls() bug
I''m very confused about generalized least square method (GLS), so please help me! I''m trying to use GLS with R to estimate an intercept and slope; vector, b = (intercept, slope). My understanding of GLS is like this: y = X %*% b + e, where Var(e) is symmetric dispersion matrix (covariance-variance matrix); V (which is theoretically known in my case). We can decompose V into t(R) %*% R, where R is an upper triangular, square-root matrix of V, and t(R) is tranpose of R. Then if we call T = Inverse(t(R)), we can multply both side of equation to make the problem as a ordinary least square. T %*% y = T %*% X %*% b + T * e To find T, "solve(t(chol(V))" should do it. However, when I looked at lm.gls() from MASS library, it uses other method of obtaining T. From my understanding, when I give V (my dispersion matrix) to the argument "W=" of lm.gls, I should also give "inverse = TRUE" (can anyone confirm that I should be using this option?). Then it uses the following method to obtain T. eW <- eigen(V, TRUE) T <- diag(eW$values ^ (-0.5)) %*% t(eW$vector) I was checking if these two methods give same answer with a small matrix; so I tried this:> tempSqRt <- matrix(c(7,0,0, 19, 76, 0, 31, 94, 51), c(3,3)) > temp <- t(tempSqRt) %*% tempSqRt # this is a symmetric dispersionmatrix> temp.e <- eigen(temp, TRUE) > temp.T <- diag(temp.e$values ^ (-0.5)) %*% t(temp.e$vector) > temp.T[,1] [,2] [,3] [1,] 0.0001090669 0.0042107823 0.006247486 [2,] 0.0004108533 -0.0272660081 0.018370016 [3,] 0.1487442294 0.0003382315 -0.002824703> solve(t(chol(temp)))[,1] [,2] [,3] [1,] 0.14285714 -5.569319e-18 3.101497e-17 [2,] -0.03571429 1.315789e-02 -6.217046e-18 [3,] -0.02100840 -2.425181e-02 1.960784e-02 So these two methods do NOT give the same answer. I guess I''m missunderstanding something. Could someone point out what''s wrong? The reason I started to compare the two methods is that when I estimate the slope and intercept with lm.gls(), the estimate was way off. A linear regression with OLS gives that slope of 26 and intercept of 5, and the regression line on scatter plot fits well, but with lm.gls(), the estimate is slope of -5 and intercept of 0.89. So the regression line doesn''t even touch the cloud of data points in the scatter plot. I was giving dispersion matrix (var-cov matrix) for "W", so I used "inverse=T". But when I use "inverse=F" (which should be wrong) using the same dispersion matrix without inversing, I get more decent estimates of slope and intercept. This is why I started to compare the two methods of obtaining T. By the way, I got MASS library from VR_6.2-5.tar.gz, and lm.gls() has a bug. In the 2nd half of the code, there is fit <- lm.fit(A %*% X, A %*% Y, method, ...) It needs to be changed to (the 3rd arg of lm.fit isn''t method): fit <- lm.fit(A %*% X, A %*% Y, method=method, ...) Also I had to change class(fit) <- c("lm.gls", class(fit)) to class(fit) <- c("lm", class(fit)) Without this change I couldn''t use summary() (I still can''t use anova() on the lm.gls() output). But I''m not sure if this is a correct way. Thank you, Naoki Naoki Takebayashi <ntakebay at bio.indiana.edu> --- Dept. of Biology, Box 90338, Duke University, Durham, NC 27708-0338 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard BSA
2001-Jun-21 18:06 UTC
[R] problem with generalized least square, and lm.gls() bug
Naoki Takebayashi <ntakebay at bio.indiana.edu> writes:> I''m very confused about generalized least square method (GLS), so > please help me! I''m trying to use GLS with R to estimate an intercept and > slope; vector, b = (intercept, slope). > > My understanding of GLS is like this: > y = X %*% b + e, where Var(e) is symmetric dispersion matrix > (covariance-variance matrix); V (which is theoretically known in my > case). We can decompose V into t(R) %*% R, where R is an upper > triangular, square-root matrix of V, and t(R) is tranpose of R. Then > if we call T = Inverse(t(R)), we can multply both side of equation to > make the problem as a ordinary least square. > > T %*% y = T %*% X %*% b + T * e > > To find T, "solve(t(chol(V))" should do it. However, when I looked at > lm.gls() from MASS library, it uses other method of obtaining T. From > my understanding, when I give V (my dispersion matrix) to the argument > "W=" of lm.gls, I should also give "inverse = TRUE" (can anyone > confirm that I should be using this option?). Then it uses the > following method to obtain T. > > eW <- eigen(V, TRUE) > T <- diag(eW$values ^ (-0.5)) %*% t(eW$vector) > > I was checking if these two methods give same answer with a small > matrix; so I tried this: > > > tempSqRt <- matrix(c(7,0,0, 19, 76, 0, 31, 94, 51), c(3,3)) > > temp <- t(tempSqRt) %*% tempSqRt # this is a symmetric dispersion > matrix > > temp.e <- eigen(temp, TRUE) > > temp.T <- diag(temp.e$values ^ (-0.5)) %*% t(temp.e$vector) > > temp.T > [,1] [,2] [,3] > [1,] 0.0001090669 0.0042107823 0.006247486 > [2,] 0.0004108533 -0.0272660081 0.018370016 > [3,] 0.1487442294 0.0003382315 -0.002824703 > > > solve(t(chol(temp))) > [,1] [,2] [,3] > [1,] 0.14285714 -5.569319e-18 3.101497e-17 > [2,] -0.03571429 1.315789e-02 -6.217046e-18 > [3,] -0.02100840 -2.425181e-02 1.960784e-02 > > So these two methods do NOT give the same answer. I guess I''m > missunderstanding something. Could someone point out what''s wrong?The T matrix is not uniquely defined. Any matrix for which crossprod(T) = solve(temp) will do. It doesn''t have to be the Choleski triangular factor. (And since R/S aren''t good at using the upper-lower diagonal matrix properties, any speed advantage of a Choleski-based method is lost anyway.) Check:> solve(temp)[,1] [,2] [,3] [1,] 2.212503e-02 3.956691e-05 -0.0004119295 [2,] 3.956691e-05 7.612803e-04 -0.0004755256 [3,] -4.119295e-04 -4.755256e-04 0.0003844675> crossprod(temp.T)[,1] [,2] [,3] [1,] 2.212503e-02 3.956691e-05 -0.0004119295 [2,] 3.956691e-05 7.612803e-04 -0.0004755256 [3,] -4.119295e-04 -4.755256e-04 0.0003844675> temp.C <- solve(t(chol(temp))) > crossprod(temp.C)[,1] [,2] [,3] [1,] 2.212503e-02 3.956691e-05 -0.0004119295 [2,] 3.956691e-05 7.612803e-04 -0.0004755256 [3,] -4.119295e-04 -4.755256e-04 0.0003844675 -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /''_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._