p.dalgaard@biostat.ku.dk 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.
-------------- next part --------------
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