From znmeb at cesmail.net  Sun Apr 10 00:52:15 2005
From: znmeb at cesmail.net (M. Edward (Ed) Borasky)
Date: Sun Apr 10 00:52:39 2005
Subject: [Rd] make check-all fails (PR#7784)
In-Reply-To: <20050409182411.C7B2CCA33@slim.kubism.ku.dk>
References: <20050409182411.C7B2CCA33@slim.kubism.ku.dk>
Message-ID: <42585C9F.7010903@cesmail.net>

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

[Rd] make check-all fails (PR#7784)
M. Edward (Ed