[EMAIL PROTECTED] wrote:
This looks more serious. 100 times machine precision is quite a large margin in these matters. Could you perhaps stick in a printout of the two terms and their difference?
I have an ATLAS build on AMD64 and it passes all the checks, but it is using ATLAS 3.7.8, so you might want to try an upgrade.
Attached ... you actually weren't very far off:
> print (max(abs(f1 - f2[common]))) [1] 2.842171e-14 > print (100*.Machine$double.eps) [1] 2.220446e-14 > stopifnot(max(abs(f1 - f2[common])) < 100*.Machine$double.eps) Error: max(abs(f1 - f2[common])) < 100 * .Machine$double.eps is not TRUE Execution halted
I'll go ahead and try ATLAS 3.7.8, but that takes a couple of hours to build. I've also got a Pentium III I can test this on.
R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.0 beta (2005-04-08), ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ## Tests of functions handling NAs in fits > ## These functions were introduced in 1.3.0. > ## They are used by lm and glm in base R, and by > ## packages MASS, rpart and survival. > > dim(airquality) [1] 153 6 > nd <- airquality[c(6,25:27), ] > > sm <- function(x) cat("length", length(x), "with", sum(is.na(x)), "NAs\n") > > # default is to omit some rows > fit <- lm(Ozone ~ ., data=airquality, na.action=na.omit) > summary(fit) Call: lm(formula = Ozone ~ ., data = airquality, na.action = na.omit) Residuals: Min 1Q Median 3Q Max -37.014 -12.284 -3.302 8.454 95.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -64.11632 23.48249 -2.730 0.00742 ** Solar.R 0.05027 0.02342 2.147 0.03411 * Wind -3.31844 0.64451 -5.149 1.23e-06 *** Temp 1.89579 0.27389 6.922 3.66e-10 *** Month -3.03996 1.51346 -2.009 0.04714 * Day 0.27388 0.22967 1.192 0.23576 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 20.86 on 105 degrees of freedom Multiple R-Squared: 0.6249, Adjusted R-squared: 0.6071 F-statistic: 34.99 on 5 and 105 DF, p-value: < 2.2e-16 > sm(fitted(fit)) length 111 with 0 NAs > sm(resid(fit)) length 111 with 0 NAs > sm(predict(fit)) length 111 with 0 NAs > (pp <- predict(fit, nd)) 6 25 26 27 NA -16.177404 1.688479 NA > > fit2 <- lm(Ozone ~ ., data=airquality, na.action=na.exclude) > summary(fit2) # same as before Call: lm(formula = Ozone ~ ., data = airquality, na.action = na.exclude) Residuals: Min 1Q Median 3Q Max -37.014 -12.284 -3.302 8.454 95.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -64.11632 23.48249 -2.730 0.00742 ** Solar.R 0.05027 0.02342 2.147 0.03411 * Wind -3.31844 0.64451 -5.149 1.23e-06 *** Temp 1.89579 0.27389 6.922 3.66e-10 *** Month -3.03996 1.51346 -2.009 0.04714 * Day 0.27388 0.22967 1.192 0.23576 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 20.86 on 105 degrees of freedom Multiple R-Squared: 0.6249, Adjusted R-squared: 0.6071 F-statistic: 34.99 on 5 and 105 DF, p-value: < 2.2e-16 > sm(fitted(fit2)) length 153 with 42 NAs > sm(resid(fit2)) length 153 with 42 NAs > sm(predict(fit2)) length 153 with 42 NAs > (pp2 <- predict(fit2, nd)) 6 25 26 27 NA -16.177404 1.688479 NA > > ## same as before: napredict is only applied to predictions on the > ## original data, following Therneau's original code (and S-PLUS). > ## However, as from R 1.8.0 there is a separate na.action arg to predict.lm() > stopifnot(all.equal(pp, pp2)) > > ## should fail > try(fit3 <- lm(Ozone ~ ., data=airquality, na.action=na.fail)) Error in na.fail.default(list(Ozone = c(41, 36, 12, 18, NA, 28, 23, 19, : missing values in object > > ## more precise tests. > f1 <- fitted(fit) > f2 <- fitted(fit2) > common <- match(names(f1), names(f2)) > stopifnot(max(abs(f1 - f2[common])) < 100*.Machine$double.eps) > stopifnot(all(is.na(f2[-common]))) > > r1 <- resid(fit) > r2 <- resid(fit2) > common <- match(names(r1), names(r2)) > stopifnot(max(abs(r1 - r2[common])) < 100*.Machine$double.eps) > stopifnot(all(is.na(r2[-common]))) > > p1 <- predict(fit) > p2 <- predict(fit2) > common <- match(names(p1), names(p2)) > stopifnot(max(abs(p1 - p2[common])) < 100*.Machine$double.eps) > stopifnot(all(is.na(p2[-common]))) > > > ### now try out glm > gfit <- glm(Ozone ~ ., data=airquality, na.action=na.omit) > summary(gfit) Call: glm(formula = Ozone ~ ., data = airquality, na.action = na.omit) Deviance Residuals: Min 1Q Median 3Q Max -37.014 -12.284 -3.302 8.454 95.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -64.11632 23.48249 -2.730 0.00742 ** Solar.R 0.05027 0.02342 2.147 0.03411 * Wind -3.31844 0.64451 -5.149 1.23e-06 *** Temp 1.89579 0.27389 6.922 3.66e-10 *** Month -3.03996 1.51346 -2.009 0.04714 * Day 0.27388 0.22967 1.192 0.23576 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for gaussian family taken to be 435.0755) Null deviance: 121802 on 110 degrees of freedom Residual deviance: 45683 on 105 degrees of freedom AIC: 997.22 Number of Fisher Scoring iterations: 2 > sm(fitted(gfit)) length 111 with 0 NAs > sm(resid(gfit)) length 111 with 0 NAs > sm(predict(gfit)) length 111 with 0 NAs > predict(gfit, nd) 6 25 26 27 NA -16.177404 1.688479 NA > (pp <- predict(gfit, nd)) 6 25 26 27 NA -16.177404 1.688479 NA > > gfit2 <- glm(Ozone ~ ., data=airquality, na.action=na.exclude) > summary(gfit2) # same as before Call: glm(formula = Ozone ~ ., data = airquality, na.action = na.exclude) Deviance Residuals: Min 1Q Median 3Q Max -37.014 -12.284 -3.302 8.454 95.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -64.11632 23.48249 -2.730 0.00742 ** Solar.R 0.05027 0.02342 2.147 0.03411 * Wind -3.31844 0.64451 -5.149 1.23e-06 *** Temp 1.89579 0.27389 6.922 3.66e-10 *** Month -3.03996 1.51346 -2.009 0.04714 * Day 0.27388 0.22967 1.192 0.23576 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for gaussian family taken to be 435.0755) Null deviance: 121802 on 110 degrees of freedom Residual deviance: 45683 on 105 degrees of freedom AIC: 997.22 Number of Fisher Scoring iterations: 2 > sm(fitted(gfit2)) length 153 with 42 NAs > sm(resid(gfit2)) length 153 with 42 NAs > sm(predict(gfit2)) length 153 with 42 NAs > (pp2 <- predict(gfit2, nd)) 6 25 26 27 NA -16.177404 1.688479 NA > stopifnot(all.equal(pp, pp2)) > > ## more precise tests. > f1 <- fitted(gfit) > f2 <- fitted(gfit2) > common <- match(names(f1), names(f2)) > print (common) [1] 1 2 3 4 7 8 9 12 13 14 15 16 17 18 19 20 21 22 [19] 23 24 28 29 30 31 38 40 41 44 47 48 49 50 51 62 63 64 [37] 66 67 68 69 70 71 73 74 76 77 78 79 80 81 82 85 86 87 [55] 88 89 90 91 92 93 94 95 99 100 101 104 105 106 108 109 110 111 [73] 112 113 114 116 117 118 120 121 122 123 124 125 126 127 128 129 130 131 [91] 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 [109] 151 152 153 > print (f1) 1 2 3 4 7 8 9 32.971099 37.113091 27.472204 16.891921 32.320560 -6.091053 -26.953745 12 13 14 15 16 17 18 35.461009 33.416076 31.035783 -5.787958 25.025785 26.074607 -23.464453 19 20 21 22 23 24 28 32.827271 13.723369 6.500011 26.103224 11.694003 7.703839 16.202396 29 30 31 38 40 41 44 45.409358 70.963391 62.723917 49.211501 59.564916 63.392633 57.551881 47 48 49 50 51 62 63 28.159115 4.382598 17.130318 29.110834 39.908867 74.042091 58.231930 64 66 67 68 69 70 71 50.319370 66.856807 53.212618 80.301978 83.724400 86.240715 70.309270 73 74 76 77 78 79 80 22.101594 31.076282 25.334479 62.718783 54.309578 72.201821 77.218308 81 82 85 86 87 88 89 54.121624 38.098263 70.456707 67.256388 52.712887 49.337475 75.253703 90 91 92 93 94 95 99 74.853015 68.333499 58.892879 46.672109 21.082539 47.154782 81.752335 100 101 104 105 106 108 109 61.708671 68.508934 49.378753 46.141704 42.765387 31.311125 47.644865 110 111 112 113 114 116 117 41.798651 40.734935 40.285066 24.876177 8.442082 46.373463 72.652242 118 120 121 122 123 124 125 65.983901 81.140660 101.389698 92.784665 86.803528 66.813059 76.464152 126 127 128 129 130 131 132 85.562396 80.164720 55.046451 22.602751 38.602215 35.466808 28.565003 133 134 135 136 137 138 139 30.487396 27.515347 17.475536 49.119123 11.994736 14.701687 49.795201 140 141 142 143 144 145 146 5.664599 24.711067 20.426534 53.013693 5.758723 19.324367 41.190110 147 148 149 151 152 153 14.189862 -19.275130 35.155615 20.525269 40.584670 18.702940 > print (f2[common]) 1 2 3 4 7 8 9 32.971099 37.113091 27.472204 16.891921 32.320560 -6.091053 -26.953745 12 13 14 15 16 17 18 35.461009 33.416076 31.035783 -5.787958 25.025785 26.074607 -23.464453 19 20 21 22 23 24 28 32.827271 13.723369 6.500011 26.103224 11.694003 7.703839 16.202396 29 30 31 38 40 41 44 45.409358 70.963391 62.723917 49.211501 59.564916 63.392633 57.551881 47 48 49 50 51 62 63 28.159115 4.382598 17.130318 29.110834 39.908867 74.042091 58.231930 64 66 67 68 69 70 71 50.319370 66.856807 53.212618 80.301978 83.724400 86.240715 70.309270 73 74 76 77 78 79 80 22.101594 31.076282 25.334479 62.718783 54.309578 72.201821 77.218308 81 82 85 86 87 88 89 54.121624 38.098263 70.456707 67.256388 52.712887 49.337475 75.253703 90 91 92 93 94 95 99 74.853015 68.333499 58.892879 46.672109 21.082539 47.154782 81.752335 100 101 104 105 106 108 109 61.708671 68.508934 49.378753 46.141704 42.765387 31.311125 47.644865 110 111 112 113 114 116 117 41.798651 40.734935 40.285066 24.876177 8.442082 46.373463 72.652242 118 120 121 122 123 124 125 65.983901 81.140660 101.389698 92.784665 86.803528 66.813059 76.464152 126 127 128 129 130 131 132 85.562396 80.164720 55.046451 22.602751 38.602215 35.466808 28.565003 133 134 135 136 137 138 139 30.487396 27.515347 17.475536 49.119123 11.994736 14.701687 49.795201 140 141 142 143 144 145 146 5.664599 24.711067 20.426534 53.013693 5.758723 19.324367 41.190110 147 148 149 151 152 153 14.189862 -19.275130 35.155615 20.525269 40.584670 18.702940 > print (f1 - f2[common]) 1 2 3 4 7 -2.131628e-14 -1.421085e-14 -1.065814e-14 -2.486900e-14 -2.842171e-14 8 9 12 13 14 -1.776357e-14 -7.105427e-15 -2.131628e-14 -2.131628e-14 -1.776357e-14 15 16 17 18 19 -1.865175e-14 -2.486900e-14 -2.131628e-14 -1.065814e-14 -2.131628e-14 20 21 22 23 24 -1.953993e-14 -2.131628e-14 -1.065814e-14 -2.131628e-14 -1.865175e-14 28 29 30 31 38 -1.065814e-14 -7.105427e-15 -1.421085e-14 -1.421085e-14 -7.105427e-15 40 41 44 47 48 0.000000e+00 -7.105427e-15 -7.105427e-15 -7.105427e-15 -3.552714e-15 49 50 51 62 63 -1.776357e-14 -1.421085e-14 -1.421085e-14 -1.421085e-14 -7.105427e-15 64 66 67 68 69 -1.421085e-14 -1.421085e-14 -7.105427e-15 -1.421085e-14 -1.421085e-14 70 71 73 74 76 -1.421085e-14 0.000000e+00 -1.065814e-14 -3.552714e-15 -3.552714e-15 77 78 79 80 81 -1.421085e-14 -1.421085e-14 -2.842171e-14 -1.421085e-14 -7.105427e-15 82 85 86 87 88 -1.421085e-14 -1.421085e-14 -1.421085e-14 -7.105427e-15 0.000000e+00 89 90 91 92 93 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.421085e-14 94 95 99 100 101 -3.552714e-15 -7.105427e-15 -1.421085e-14 -7.105427e-15 0.000000e+00 104 105 106 108 109 -7.105427e-15 -7.105427e-15 -1.421085e-14 -1.065814e-14 -1.421085e-14 110 111 112 113 114 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.065814e-14 -8.881784e-15 116 117 118 120 121 -1.421085e-14 -2.842171e-14 -1.421085e-14 0.000000e+00 -1.421085e-14 122 123 124 125 126 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.421085e-14 127 128 129 130 131 0.000000e+00 -7.105427e-15 0.000000e+00 -7.105427e-15 -1.421085e-14 132 133 134 135 136 -1.776357e-14 -1.776357e-14 -3.552714e-15 -1.065814e-14 -2.131628e-14 137 138 139 140 141 -1.421085e-14 -1.421085e-14 -1.421085e-14 -1.687539e-14 -1.065814e-14 142 143 144 145 146 -1.776357e-14 -1.421085e-14 -2.042810e-14 -1.421085e-14 -7.105427e-15 147 148 149 151 152 -1.598721e-14 -1.065814e-14 -2.842171e-14 -1.065814e-14 -2.131628e-14 153 -1.776357e-14 > print (abs(f1 - f2[common])) 1 2 3 4 7 8 2.131628e-14 1.421085e-14 1.065814e-14 2.486900e-14 2.842171e-14 1.776357e-14 9 12 13 14 15 16 7.105427e-15 2.131628e-14 2.131628e-14 1.776357e-14 1.865175e-14 2.486900e-14 17 18 19 20 21 22 2.131628e-14 1.065814e-14 2.131628e-14 1.953993e-14 2.131628e-14 1.065814e-14 23 24 28 29 30 31 2.131628e-14 1.865175e-14 1.065814e-14 7.105427e-15 1.421085e-14 1.421085e-14 38 40 41 44 47 48 7.105427e-15 0.000000e+00 7.105427e-15 7.105427e-15 7.105427e-15 3.552714e-15 49 50 51 62 63 64 1.776357e-14 1.421085e-14 1.421085e-14 1.421085e-14 7.105427e-15 1.421085e-14 66 67 68 69 70 71 1.421085e-14 7.105427e-15 1.421085e-14 1.421085e-14 1.421085e-14 0.000000e+00 73 74 76 77 78 79 1.065814e-14 3.552714e-15 3.552714e-15 1.421085e-14 1.421085e-14 2.842171e-14 80 81 82 85 86 87 1.421085e-14 7.105427e-15 1.421085e-14 1.421085e-14 1.421085e-14 7.105427e-15 88 89 90 91 92 93 0.000000e+00 1.421085e-14 1.421085e-14 1.421085e-14 1.421085e-14 1.421085e-14 94 95 99 100 101 104 3.552714e-15 7.105427e-15 1.421085e-14 7.105427e-15 0.000000e+00 7.105427e-15 105 106 108 109 110 111 7.105427e-15 1.421085e-14 1.065814e-14 1.421085e-14 1.421085e-14 1.421085e-14 112 113 114 116 117 118 1.421085e-14 1.065814e-14 8.881784e-15 1.421085e-14 2.842171e-14 1.421085e-14 120 121 122 123 124 125 0.000000e+00 1.421085e-14 1.421085e-14 1.421085e-14 1.421085e-14 1.421085e-14 126 127 128 129 130 131 1.421085e-14 0.000000e+00 7.105427e-15 0.000000e+00 7.105427e-15 1.421085e-14 132 133 134 135 136 137 1.776357e-14 1.776357e-14 3.552714e-15 1.065814e-14 2.131628e-14 1.421085e-14 138 139 140 141 142 143 1.421085e-14 1.421085e-14 1.687539e-14 1.065814e-14 1.776357e-14 1.421085e-14 144 145 146 147 148 149 2.042810e-14 1.421085e-14 7.105427e-15 1.598721e-14 1.065814e-14 2.842171e-14 151 152 153 1.065814e-14 2.131628e-14 1.776357e-14 > print (max(abs(f1 - f2[common]))) [1] 2.842171e-14 > print (100*.Machine$double.eps) [1] 2.220446e-14 > stopifnot(max(abs(f1 - f2[common])) < 100*.Machine$double.eps) Error: max(abs(f1 - f2[common])) < 100 * .Machine$double.eps is not TRUE Execution halted
______________________________________________ R-devel@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-devel