[R] Problem with lm.resid() when weights are provided

Fox, John j|ox @end|ng |rom mcm@@ter@c@
Mon Sep 17 19:16:19 CEST 2018


Dear Hamed,

> -----Original Message-----
> From: Hamed Ha [mailto:hamedhaseli using gmail.com]
> Sent: Monday, September 17, 2018 3:56 AM
> To: Fox, John <jfox using mcmaster.ca>
> Cc: r-help using r-project.org
> Subject: Re: [R] Problem with lm.resid() when weights are provided
> 
> H i John,
> 
> 
> Thank you for your reply.
> 
> 
> I see your point, thanks. I checked lm.wfit() and realised that there is a tol
> parameter that is already set to 10^-7. This is not even the half decimal to the
> machine precision. Furthermore, plying with tol parameter does not solve the
> problem, as far as I checked.

tol plays a different role in lm.wfit(). It's for the QR decomposition (done in C code), I suppose to determine the rank of the weighted model matrix. Generally in this kind of context, you'd use something like the square root of the machine double epsilon to define a number that's effectively 0, and the tolerance used here isn't too far off that -- about an order of magnitude larger.
 
I'm not an expert in computer arithmetic or numerical linear algebra, so I don't have anything more to say about this.

> 
> 
> I still see this issue as critical and we should report it to the R core team to be
> investigated more. What do you think?

I don't think that it's a critical issue because it isn't sensible to specify nonzero weights so close to 0. A simple solution is to change these weights to 0 in your code calling lm().

That said, I suppose that it might be better to make lm.wfit() more robust to near-zero weights. If you feel strongly about this, you can file a bug report, but I'm not interested in pursuing it.

Best,
 John

> 
> 
> Regards,
> Hamed.
> 
> 
> On Fri, 14 Sep 2018 at 22:46, Fox, John <jfox using mcmaster.ca
> <mailto:jfox using mcmaster.ca> > wrote:
> 
> 
> 	Dear Hamed,
> 
> 	When you post a question to r-help, generally you should cc
> subsequent messages there as well, as I've done to this response.
> 
> 	The algorithm that lm() uses is much more numerically stable than
> inverting the weighted sum-of-squares-and-product matrix. If you want to see
> how the computations are done, look at lm.wfit(), in which the residuals and
> fits are computed as
> 
> 	    z$residuals <- z$residuals/wts
> 	    z$fitted.values <- y - z$residuals
> 
> 	Zero weights are handled specially, and your tiny weights are thus the
> source of the problem. When you divide by a number less than the machine
> double-epsilon, you can't expect numerically stable results. I suppose that
> lm.wfit() could check for 0 weights to a tolerance rather than exactly.
> 
> 	John
> 
> 	> -----Original Message-----
> 	> From: Hamed Ha [mailto:hamedhaseli using gmail.com
> <mailto:hamedhaseli using gmail.com> ]
> 	> Sent: Friday, September 14, 2018 5:34 PM
> 	> To: Fox, John <jfox using mcmaster.ca <mailto:jfox using mcmaster.ca> >
> 	> Subject: Re: [R] Problem with lm.resid() when weights are provided
> 	>
> 	> Hi John,
> 	>
> 	> Thank you for your reply.
> 	>
> 	> I agree that the small weights are the potential source of the
> instability in the
> 	> result. I also suspected that there are some failure/bugs in the actual
> 	> algorithm that R uses for fitting the model. I remember that at some
> points I
> 	> checked the theoretical estimation of the parameters, solve(t(x)
> %*% w %*%
> 	> x) %*% t(x) %*% w %*% y, (besides the point that I had to set tol
> parameter in
> 	> solve() to a super small value) and realised  that lm() and the
> theoretical
> 	> results match together. That is the parameter estimation is right in
> R.
> 	> Moreover, I checked the predictions, predict(lm.fit), and it was right.
> Then the
> 	> only source of error remained was resid() function. I further checked
> this
> 	> function and it is nothing more than calling a sub-element from and
> lm() fit.
> 	> Putting all together, I think that there is something wrong/bug/miss-
> 	> configuration in the lm() algorithm and I highly recommend the R
> core team to
> 	> fix that.
> 	>
> 	> Please feel free to contact me for more details if required.
> 	>
> 	> Warm regards,
> 	> Hamed.
> 	>
> 	>
> 	>
> 	>
> 	>
> 	>
> 	>
> 	>
> 	>
> 	> On Fri, 14 Sep 2018 at 13:35, Fox, John <jfox using mcmaster.ca
> <mailto:jfox using mcmaster.ca>
> 	> <mailto:jfox using mcmaster.ca <mailto:jfox using mcmaster.ca> > > wrote:
> 	>
> 	>
> 	>       Dear Hamed,
> 	>
> 	>       I don't think that anyone has picked up on this problem.
> 	>
> 	>       What's peculiar about your weights is that several are 0 within
> 	> rounding error but not exactly 0:
> 	>
> 	>       > head(df)
> 	>                  y          x       weight
> 	>       1  1.5115614  0.5520924 2.117337e-34
> 	>       2 -0.6365313 -0.1259932 2.117337e-34
> 	>       3  0.3778278  0.4209538 4.934135e-31
> 	>       4  3.0379232  1.4031545 2.679495e-24
> 	>       5  1.5364652  0.4607686 2.679495e-24
> 	>       6 -2.3772787 -0.7396358 6.244160e-21
> 	>
> 	>       I can reproduce the results that you report:
> 	>
> 	>       > (mod.1 <- lm(y ~ x, data=df))
> 	>
> 	>       Call:
> 	>       lm(formula = y ~ x, data = df)
> 	>
> 	>       Coefficients:
> 	>       (Intercept)            x
> 	>          -0.04173      2.03790
> 	>
> 	>       > max(resid(mod.1))
> 	>       [1] 1.14046
> 	>       > (mod.2 <- lm(y ~ x, data=df, weights=weight))
> 	>
> 	>       Call:
> 	>       lm(formula = y ~ x, data = df, weights = weight)
> 	>
> 	>       Coefficients:
> 	>       (Intercept)            x
> 	>          -0.05786      1.96087
> 	>
> 	>       > max(resid(mod.2))
> 	>       [1] 36.84939
> 	>
> 	>       But the problem disappears when the tiny nonzero weight are set
> to 0:
> 	>
> 	>       > df2 <- df
> 	>       > df2$weight <- zapsmall(df2$weight)
> 	>       > head(df2)
> 	>                  y          x weight
> 	>       1  1.5115614  0.5520924      0
> 	>       2 -0.6365313 -0.1259932      0
> 	>       3  0.3778278  0.4209538      0
> 	>       4  3.0379232  1.4031545      0
> 	>       5  1.5364652  0.4607686      0
> 	>       6 -2.3772787 -0.7396358      0
> 	>       > (mod.3 <- update(mod.2, data=df2))
> 	>
> 	>       Call:
> 	>       lm(formula = y ~ x, data = df2, weights = weight)
> 	>
> 	>       Coefficients:
> 	>       (Intercept)            x
> 	>          -0.05786      1.96087
> 	>
> 	>       > max(resid(mod.3))
> 	>       [1] 1.146663
> 	>
> 	>       I don't know exactly why this happens, but suspect numerical
> 	> instability produced by the near-zero weights, which are smaller
> than the
> 	> machine double-epsilon
> 	>
> 	>       > .Machine$double.neg.eps
> 	>       [1] 1.110223e-16
> 	>
> 	>       The problem also disappears, e.g., if the tiny weight are set to 1e-
> 15
> 	> rather than 0.
> 	>
> 	>       I hope this helps,
> 	>        John
> 	>
> 	>       -----------------------------------------------------------------
> 	>       John Fox
> 	>       Professor Emeritus
> 	>       McMaster University
> 	>       Hamilton, Ontario, Canada
> 	>       Web: https://socialsciences.mcmaster.ca/jfox/
> 	>
> 	>
> 	>
> 	>       > -----Original Message-----
> 	>       > From: R-help [mailto:r-help-bounces using r-project.org <mailto:r-
> help-bounces using r-project.org>  <mailto:r-help- <mailto:r-help->
> 	> bounces using r-project.org <mailto:bounces using r-project.org> > ] On
> Behalf Of Hamed Ha
> 	>       > Sent: Tuesday, September 11, 2018 8:39 AM
> 	>       > To: r-help using r-project.org <mailto:r-help using r-project.org>
> <mailto:r-help using r-project.org <mailto:r-help using r-project.org> >
> 	>       > Subject: [R] Problem with lm.resid() when weights are provided
> 	>       >
> 	>       > Dear R Help Team.
> 	>       >
> 	>       > I get some weird results when I use the lm function with weight.
> The
> 	> issue can
> 	>       > be reproduced by the example below:
> 	>       >
> 	>       >
> 	>       > The input data is (weights are intentionally designed to reflect
> some
> 	>       > structures in the data)
> 	>       >
> 	>       >
> 	>       > > df
> 	>       > y x weight
> 	>       >  1.51156139  0.55209240 2.117337e-34
> 	>       > -0.63653132 -0.12599316 2.117337e-34
> 	>       >  0.37782776  0.42095384 4.934135e-31
> 	>       >  3.03792318  1.40315446 2.679495e-24
> 	>       >  1.53646523  0.46076858 2.679495e-24
> 	>       > -2.37727874 -0.73963576 6.244160e-21
> 	>       >  0.37183065  0.20407468 1.455107e-17
> 	>       > -1.53917553 -0.95519361 1.455107e-17
> 	>       >  1.10926675  0.03897129 3.390908e-14
> 	>       > -0.37786333 -0.17523593 3.390908e-14
> 	>       >  2.43973603  0.97970095 7.902000e-11
> 	>       > -0.35432394 -0.03742559 7.902000e-11
> 	>       >  2.19296613  1.00355263 4.289362e-04
> 	>       >  0.49845532  0.34816207 4.289362e-04
> 	>       >  1.25005260  0.76306225 5.000000e-01
> 	>       >  0.84360691  0.45152356 5.000000e-01
> 	>       >  0.29565993  0.53880068 5.000000e-01
> 	>       > -0.54081334 -0.28104525 5.000000e-01
> 	>       >  0.83612836 -0.12885659 9.995711e-01
> 	>       > -1.42526769 -0.87107631 9.999998e-01
> 	>       >  0.10204789 -0.11649899 1.000000e+00
> 	>       >  1.14292898  0.37249631 1.000000e+00
> 	>       > -3.02942081 -1.28966997 1.000000e+00
> 	>       > -1.37549764 -0.74676145 1.000000e+00
> 	>       > -2.00118016 -0.55182759 1.000000e+00
> 	>       > -4.24441674 -1.94603608 1.000000e+00
> 	>       >  1.17168144  1.00868008 1.000000e+00
> 	>       >  2.64007761  1.26333069 1.000000e+00
> 	>       >  1.98550114  1.18509599 1.000000e+00
> 	>       > -0.58941683 -0.61972416 9.999998e-01
> 	>       > -4.57559611 -2.30914920 9.995711e-01
> 	>       > -0.82610544 -0.39347576 9.995711e-01
> 	>       > -0.02768220  0.20076910 9.995711e-01
> 	>       >  0.78186399  0.25690215 9.995711e-01
> 	>       > -0.88314153 -0.20200148 5.000000e-01
> 	>       > -4.17076452 -2.03547588 5.000000e-01
> 	>       >  0.93373070  0.54190626 4.289362e-04
> 	>       > -0.08517734  0.17692491 4.289362e-04
> 	>       > -4.47546619 -2.14876688 4.289362e-04
> 	>       > -1.65509103 -0.76898087 4.289362e-04
> 	>       > -0.39403030 -0.12689705 4.289362e-04
> 	>       >  0.01203300 -0.18689898 1.841442e-07
> 	>       > -4.82762639 -2.31391121 1.841442e-07
> 	>       > -0.72658380 -0.39751171 3.397282e-14
> 	>       > -2.35886866 -1.01082109 0.000000e+00
> 	>       > -2.03762707 -0.96439902 0.000000e+00
> 	>       >  0.90115123  0.60172286 0.000000e+00
> 	>       >  1.55999194  0.83433953 0.000000e+00
> 	>       >  3.07994058  1.30942776 0.000000e+00
> 	>       >  1.78871462  1.10605530 0.000000e+00
> 	>       >
> 	>       >
> 	>       >
> 	>       > Running simple linear model returns:
> 	>       >
> 	>       > > lm(y~x,data=df)
> 	>       >
> 	>       > Call:
> 	>       > lm(formula = y ~ x, data = df)
> 	>       >
> 	>       > Coefficients:
> 	>       > (Intercept)            x
> 	>       >    -0.04173      2.03790
> 	>       >
> 	>       > and
> 	>       > > max(resid(lm(y~x,data=df)))
> 	>       > [1] 1.14046
> 	>       >
> 	>       >
> 	>       > *HOWEVER if I use the weighted model then:*
> 	>       >
> 	>       > lm(formula = y ~ x, data = df, weights = df$weights)
> 	>       >
> 	>       > Coefficients:
> 	>       > (Intercept)            x
> 	>       >    -0.05786      1.96087
> 	>       >
> 	>       > and
> 	>       > > max(resid(lm(y~x,data=df,weights=df$weights)))
> 	>       > [1] 60.91888
> 	>       >
> 	>       >
> 	>       > as you see, the estimation of the coefficients are nearly the
> same
> 	> but the
> 	>       > resid() function returns a giant residual (I have some cases
> where
> 	> the value is
> 	>       > much much higher). Further, if I calculate the residuals by
> simply
> 	>       > predict(lm(y~x,data=df,weights=df$weights))-df$y then I get the
> true
> 	> value for
> 	>       > the residuals.
> 	>       >
> 	>       >
> 	>       > Thanks.
> 	>       >
> 	>       > Please do not hesitate to contact me for more details.
> 	>       > Regards,
> 	>       > Hamed.
> 	>       >
> 	>       >       [[alternative HTML version deleted]]
> 	>       >
> 	>       > ______________________________________________
> 	>       > R-help using r-project.org <mailto:R-help using r-project.org>
> <mailto:R-help using r-project.org <mailto:R-help using r-project.org> >  mailing list --
> 	> To UNSUBSCRIBE and more, see
> 	>       > https://stat.ethz.ch/mailman/listinfo/r-help
> 	>       > PLEASE do read the posting guide http://www.R-
> project.org/posting-
> 	>       > guide.html
> 	>       > and provide commented, minimal, self-contained, reproducible
> 	> code.
> 	>
> 
> 



More information about the R-help mailing list