R Under development (unstable) (2023-06-09 r84528 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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 an HTML browser interface to help. Type 'q()' to quit R. > > library(robustbase) > > ## instead of relying on system.file("test-tools-1.R", package="Matrix"): > source(system.file("xtraR/test-tools.R", package = "robustbase")) # assert.EQ() etc > > #### Poisson examples from Eva Cantoni's paper > > ### Using Possum Data > ### ================ > > data(possumDiv) > > ## Try to follow closely Cantoni & Ronchetti(2001), JASA > dim(X <- possum.mat[, -1]) # 151 13 [1] 151 13 > str(y <- possum.mat[, "Diversity"]) int [1:151] 3 2 1 2 3 2 3 2 0 0 ... > ##--- reduce the matrix from singularity ourselves: > X. <- possum.mat[, -c(1, match(c("E.nitens", "NW-NE"), colnames(possum.mat)))] > dim(X.)# 151 11 [1] 151 11 > > ## "classical via robust: c = Inf : > Inf. <- 1e5 ## --- FIXME > > ## The following used to fail because glm.fit() returns NA coefficients > ## now fine .. keep this as test! > glm.cr <- glmrob(y ~ X, family = "poisson", tcc = Inf.) initial start 'theta' has NA's; eliminating columns X[, j]; j = 10, 14 > (scr <- summary(glm.cr)) Call: glmrob(formula = y ~ X, family = "poisson", tcc = Inf.) Coefficients: (2 not defined because of singularities) Estimate Std. Error z value Pr(>|z|) (Intercept) -1.32093 0.39820 -3.317 0.000909 *** XShrubs 0.01192 0.02195 0.543 0.587006 XStumps -0.27241 0.28592 -0.953 0.340727 XStags 0.04023 0.01120 3.590 0.000330 *** XBark 0.03989 0.01439 2.772 0.005571 ** XHabitat 0.07173 0.03814 1.881 0.059999 . XBAcacia 0.01764 0.01060 1.664 0.096044 . XE.regnans -0.11492 0.27242 -0.422 0.673132 XE.delegatensis -0.13027 0.30883 -0.422 0.673171 XNW-NE 0.48891 0.24747 1.976 0.048193 * XNW-SE 0.55566 0.23387 2.376 0.017506 * XSE-SW 0.60585 0.22725 2.666 0.007676 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: All 151 weights are ~= 1. Number of observations: 151 Fitted by method 'Mqle' (in 1 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: acc 1e-04 maxit tcc 5e+01 1e+05 test.acc "coef" > > scl <- summary(glm.cl <- glm (Diversity ~ . , data=possumDiv, family=poisson)) > sc2 <- summary(glm.c2 <- glmrob(Diversity ~ . , data=possumDiv, family=poisson, tcc = Inf.)) > MMg <- model.matrix(glm.cl) > > assert.EQ(coef(scl), coef(sc2), tol = 6e-6, giveRE=TRUE) # 1.37e-6 Mean relative difference: 1.363273e-06 > dnms <- list(colnames(MMg), c("Estimate", "Std. Error", "z value", "Pr(>|z|)")) > cf.sc <- array(c(-0.9469439, 0.01192096, -0.2724059, 0.04022862, 0.03988606, 0.07173483, + 0.01763833, -0.01534376, 0.1149216, 0.06675529, 0.1169463, -0.4889071, + ## SE + 0.2655031, 0.02194661, 0.2859216, 0.01120463, 0.01438884, 0.03814053, + 0.01059779, 0.1916126, 0.2724202, 0.1901612, 0.1902903, 0.2474653, + ## z val + -3.566603, 0.5431798, -0.9527294, 3.590356, 2.772014, 1.880803, + 1.664341, -0.08007701, 0.421854, 0.3510457, 0.6145675, -1.975659, + ## P val + 0.0003616393, 0.587006, 0.3407272, 0.0003302263, 0.00557107, 0.05999869, + 0.09604432, 0.936176, 0.6731316, 0.7255541, 0.5388404, 0.04819339), + dim = c(12L, 4L), dimnames = dnms) > assert.EQ(cf.sc, coef(sc2), tol = 4e-7, giveRE=TRUE) # 8.48e-8 Mean relative difference: 8.640753e-08 > > > ## c = 2.0 > summary(g2 <- glmrob(Diversity ~ . , data=possumDiv, family=poisson, tcc = 2.0, trace=TRUE)) Initial theta: (Intercept) Shrubs Stumps -0.94694387 0.01192096 -0.27240588 Stags Bark Habitat 0.04022862 0.03988606 0.07173483 BAcacia eucalyptusdelegatensis eucalyptusnitens 0.01763833 -0.01534376 0.11492155 aspectNW-SE aspectSE-SW aspectSW-NW 0.06675529 0.11694626 -0.48890705 it | d{theta} | rel.change 1 | 0.039 -0.00046 0.012 -0.00032 -0.00038 -0.001 -0.0001 -0.0039 0.0012 -0.002 -0.0074 -0.0025 | 0.0372002 2 | -0.0024 8.4e-5 -0.0023 -3.1e-6 7.9e-6 1e-5 -1.4e-5 0.00028 -0.00067 0.00027 0.0021 0.0016 | 0.0040033 3 | 0.00027 -2.6e-6 0.00019 -1.1e-6 -3.8e-6 -7e-6 -3.2e-7 -2.3e-5 6.6e-5 -3.8e-5 -0.00015 -0.00013 | 0.00036205 4 | -2.2e-5 1.1e-7 -2e-5 -2.5e-8 3.9e-7 7.3e-7 -2.1e-8 1.1e-6 -8.7e-6 4.4e-6 1.6e-5 1.4e-5 | 3.47231e-05 Call: glmrob(formula = Diversity ~ ., family = poisson, data = possumDiv, trace.lev = TRUE, tcc = 2) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.91031 0.26641 -3.417 0.000633 *** Shrubs 0.01154 0.02202 0.524 0.600326 Stumps -0.26273 0.28635 -0.918 0.358869 Stags 0.03990 0.01125 3.548 0.000388 *** Bark 0.03951 0.01443 2.738 0.006177 ** Habitat 0.07073 0.03826 1.849 0.064529 . BAcacia 0.01752 0.01063 1.648 0.099385 . eucalyptusdelegatensis -0.01899 0.19246 -0.099 0.921395 eucalyptusnitens 0.11548 0.27308 0.423 0.672389 aspectNW-SE 0.06500 0.19060 0.341 0.733069 aspectSE-SW 0.11145 0.19087 0.584 0.559266 aspectSW-NW -0.48990 0.24859 -1.971 0.048759 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: 148 weights are ~= 1. The remaining 3 ones are 59 110 139 0.9934 0.8551 0.8867 Number of observations: 151 Fitted by method 'Mqle' (in 4 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: acc 1e-04 maxit tcc 50 2 test.acc "coef" > > ## c = 1.6 > glm.r <- glmrob(Diversity ~ . , data=possumDiv, family=poisson, tcc = 1.6, trace=TRUE) Initial theta: (Intercept) Shrubs Stumps -0.94694387 0.01192096 -0.27240588 Stags Bark Habitat 0.04022862 0.03988606 0.07173483 BAcacia eucalyptusdelegatensis eucalyptusnitens 0.01763833 -0.01534376 0.11492155 aspectNW-SE aspectSE-SW aspectSW-NW 0.06675529 0.11694626 -0.48890705 it | d{theta} | rel.change 1 | 0.056 -0.0013 0.028 -4.5e-5 -0.00031 -0.0011 0.00028 -0.0019 0.0057 -0.0041 -0.026 -0.023 | 0.0641836 2 | -0.0079 0.00028 -0.0065 -9.7e-5 4.9e-5 0.00021 -0.0001 0.00018 -0.0014 0.00084 0.0057 0.0067 | 0.0126651 3 | 0.0013 -3.9e-5 0.00083 3.9e-6 -9.7e-6 1.3e-5 -2.6e-6 -9.2e-5 -0.00015 -0.00016 -0.00073 -0.00093 | 0.0017913 4 | -0.00018 5.5e-6 -0.00013 -2.4e-6 2.5e-6 1.8e-6 -7.2e-7 3.6e-6 2e-5 2.3e-5 0.00013 0.00015 | 0.00027738 5 | 2.8e-5 -7.7e-7 1.8e-5 2.4e-7 -4e-7 2.3e-7 -1.3e-7 -1.3e-6 -7.6e-6 -3.9e-6 -1.7e-5 -2.2e-5 | 4.13599e-05 > (s.16 <- summary(glm.r)) Call: glmrob(formula = Diversity ~ ., family = poisson, data = possumDiv, trace.lev = TRUE, tcc = 1.6) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.89790 0.26823 -3.348 0.000815 *** Shrubs 0.01091 0.02220 0.491 0.623173 Stumps -0.25066 0.28749 -0.872 0.383267 Stags 0.04009 0.01132 3.541 0.000398 *** Bark 0.03962 0.01451 2.730 0.006332 ** Habitat 0.07087 0.03851 1.840 0.065759 . BAcacia 0.01781 0.01070 1.665 0.095914 . eucalyptusdelegatensis -0.01720 0.19374 -0.089 0.929246 eucalyptusnitens 0.11904 0.27466 0.433 0.664712 aspectNW-SE 0.06334 0.19127 0.331 0.740529 aspectSE-SW 0.09620 0.19206 0.501 0.616464 aspectSW-NW -0.50651 0.25065 -2.021 0.043304 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: 146 weights are ~= 1. The remaining 5 ones are 14 59 110 133 139 0.9757 0.7914 0.6840 0.8896 0.7050 Number of observations: 151 Fitted by method 'Mqle' (in 5 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: acc tcc 0.0001 1.6000 maxit 50 test.acc "coef" > str(glm.r) List of 28 $ coefficients : Named num [1:12] -0.8979 0.0109 -0.2507 0.0401 0.0396 ... ..- attr(*, "names")= chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... $ residuals : Named num [1:151] -0.2282 0.2132 -0.7498 -0.0211 -0.2381 ... ..- attr(*, "names")= chr [1:151] "1" "2" "3" "4" ... $ fitted.values : Named num [1:151] 3.42 1.72 2.08 2.03 3.44 ... ..- attr(*, "names")= chr [1:151] "1" "2" "3" "4" ... $ w.r : num [1:151] 1 1 1 1 1 1 1 1 1 1 ... $ w.x : num [1:151] 1 1 1 1 1 1 1 1 1 1 ... $ ni : num [1:151] 1 1 1 1 1 1 1 1 1 1 ... $ dispersion : num 1 $ cov : num [1:12, 1:12] 0.071945 -0.002564 0.005427 -0.000343 -0.000825 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... $ matM : num [1:12, 1:12] 1.3222 6.9774 0.0894 13.3873 12.2805 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... $ matQ : num [1:12, 1:12] 1.204 6.355 0.081 12.243 11.213 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... .. ..$ : chr [1:12] "(Intercept)" "Shrubs" "Stumps" "Stags" ... $ tcc : num 1.6 $ family :List of 13 ..$ family : chr "poisson" ..$ link : chr "log" ..$ linkfun :function (mu) ..$ linkinv :function (eta) ..$ variance :function (mu) ..$ dev.resids:function (y, mu, wt) ..$ aic :function (y, n, mu, wt, dev) ..$ mu.eta :function (eta) ..$ initialize: expression({ if (any(y < 0)) stop("negative values not allowed for the 'Poisson' family") n <- rep.int(1, nobs| __truncated__ ..$ validmu :function (mu) ..$ valideta :function (eta) ..$ simulate :function (object, nsim) ..$ dispersion: num 1 ..- attr(*, "class")= chr "family" $ linear.predictors: Named num [1:151] 1.23 0.543 0.733 0.708 1.236 ... ..- attr(*, "names")= chr [1:151] "1" "2" "3" "4" ... $ deviance : NULL $ iter : int 5 $ y : Named int [1:151] 3 2 1 2 3 2 3 2 0 0 ... ..- attr(*, "names")= chr [1:151] "1" "2" "3" "4" ... $ converged : logi TRUE $ model :'data.frame': 151 obs. of 9 variables: ..$ Diversity : int [1:151] 3 2 1 2 3 2 3 2 0 0 ... ..$ Shrubs : int [1:151] 6 5 7 6 5 3 6 13 5 8 ... ..$ Stumps : int [1:151] 1 0 0 0 0 0 0 0 0 0 ... ..$ Stags : int [1:151] 12 15 6 14 16 16 9 20 7 4 ... ..$ Bark : int [1:151] 29 12 26 16 11 6 10 4 13 9 ... ..$ Habitat : int [1:151] 9 2 2 8 8 10 8 8 3 1 ... ..$ BAcacia : int [1:151] 31 4 8 16 20 31 16 17 0 8 ... ..$ eucalyptus: Factor w/ 3 levels "regnans","delegatensis",..: 1 1 1 1 2 1 2 1 1 1 ... ..$ aspect : Factor w/ 4 levels "NW-NE","NW-SE",..: 4 3 1 4 3 2 4 3 4 3 ... ..- attr(*, "terms")=Classes 'terms', 'formula' language Diversity ~ Shrubs + Stumps + Stags + Bark + Habitat + BAcacia + eucalyptus + aspect .. .. ..- attr(*, "variables")= language list(Diversity, Shrubs, Stumps, Stags, Bark, Habitat, BAcacia, eucalyptus, aspect) .. .. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ... .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. ..$ : chr [1:9] "Diversity" "Shrubs" "Stumps" "Stags" ... .. .. .. .. ..$ : chr [1:8] "Shrubs" "Stumps" "Stags" "Bark" ... .. .. ..- attr(*, "term.labels")= chr [1:8] "Shrubs" "Stumps" "Stags" "Bark" ... .. .. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1 .. .. ..- attr(*, "intercept")= int 1 .. .. ..- attr(*, "response")= int 1 .. .. ..- attr(*, ".Environment")= .. .. ..- attr(*, "predvars")= language list(Diversity, Shrubs, Stumps, Stags, Bark, Habitat, BAcacia, eucalyptus, aspect) .. .. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ... .. .. .. ..- attr(*, "names")= chr [1:9] "Diversity" "Shrubs" "Stumps" "Stags" ... $ call : language glmrob(formula = Diversity ~ ., family = poisson, data = possumDiv, trace.lev = TRUE, tcc = 1.6) $ formula :Class 'formula' language Diversity ~ . .. ..- attr(*, ".Environment")= $ terms :Classes 'terms', 'formula' language Diversity ~ Shrubs + Stumps + Stags + Bark + Habitat + BAcacia + eucalyptus + aspect .. ..- attr(*, "variables")= language list(Diversity, Shrubs, Stumps, Stags, Bark, Habitat, BAcacia, eucalyptus, aspect) .. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ... .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : chr [1:9] "Diversity" "Shrubs" "Stumps" "Stags" ... .. .. .. ..$ : chr [1:8] "Shrubs" "Stumps" "Stags" "Bark" ... .. ..- attr(*, "term.labels")= chr [1:8] "Shrubs" "Stumps" "Stags" "Bark" ... .. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1 .. ..- attr(*, "intercept")= int 1 .. ..- attr(*, "response")= int 1 .. ..- attr(*, ".Environment")= .. ..- attr(*, "predvars")= language list(Diversity, Shrubs, Stumps, Stags, Bark, Habitat, BAcacia, eucalyptus, aspect) .. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ... .. .. ..- attr(*, "names")= chr [1:9] "Diversity" "Shrubs" "Stumps" "Stags" ... $ data :'data.frame': 151 obs. of 9 variables: ..$ Diversity : int [1:151] 3 2 1 2 3 2 3 2 0 0 ... ..$ Shrubs : int [1:151] 6 5 7 6 5 3 6 13 5 8 ... ..$ Stumps : int [1:151] 1 0 0 0 0 0 0 0 0 0 ... ..$ Stags : int [1:151] 12 15 6 14 16 16 9 20 7 4 ... ..$ Bark : int [1:151] 29 12 26 16 11 6 10 4 13 9 ... ..$ Habitat : int [1:151] 9 2 2 8 8 10 8 8 3 1 ... ..$ BAcacia : int [1:151] 31 4 8 16 20 31 16 17 0 8 ... ..$ eucalyptus: Factor w/ 3 levels "regnans","delegatensis",..: 1 1 1 1 2 1 2 1 1 1 ... ..$ aspect : Factor w/ 4 levels "NW-NE","NW-SE",..: 4 3 1 4 3 2 4 3 4 3 ... $ offset : NULL $ control :List of 4 ..$ acc : num 1e-04 ..$ test.acc: chr "coef" ..$ maxit : num 50 ..$ tcc : num 1.6 $ method : chr "Mqle" $ prior.weights : num [1:151] 1 1 1 1 1 1 1 1 1 1 ... $ contrasts :List of 2 ..$ eucalyptus: chr "contr.treatment" ..$ aspect : chr "contr.treatment" $ xlevels :List of 2 ..$ eucalyptus: chr [1:3] "regnans" "delegatensis" "nitens" ..$ aspect : chr [1:4] "NW-NE" "NW-SE" "SE-SW" "SW-NW" - attr(*, "class")= chr [1:2] "glmrob" "glm" > > ## Now with *smaller* X (two variables less): > glm.c2 <- glmrob(y ~ X., family = "poisson", tcc = Inf.) > summary(glm.c2) Call: glmrob(formula = y ~ X., family = "poisson", tcc = Inf.) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.83202 0.35582 -2.338 0.01937 * X.Shrubs 0.01192 0.02195 0.543 0.58701 X.Stumps -0.27241 0.28592 -0.953 0.34073 X.Stags 0.04023 0.01120 3.590 0.00033 *** X.Bark 0.03989 0.01439 2.772 0.00557 ** X.Habitat 0.07173 0.03814 1.881 0.06000 . X.BAcacia 0.01764 0.01060 1.664 0.09604 . X.E.regnans -0.11492 0.27242 -0.422 0.67313 X.E.delegatensis -0.13027 0.30883 -0.422 0.67317 X.NW-SE 0.06676 0.19016 0.351 0.72555 X.SE-SW 0.11695 0.19029 0.615 0.53884 X.SW-NW -0.48891 0.24747 -1.976 0.04819 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: All 151 weights are ~= 1. Number of observations: 151 Fitted by method 'Mqle' (in 1 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: acc 1e-04 maxit tcc 5e+01 1e+05 test.acc "coef" > > ## c = 1.6, x-weights, as in Cantoni-Ronchetti > glm.r2 <- glmrob(y ~ X., family = "poisson", + tcc = 1.6, weights.on.x = "hat") > > ## Now the same, for the direct possum data (no matrix), > ## This indeed gives the same coefficients as in > ## Table 3 of Cantoni+Ronchetti(2001): .. (tech.rep.): > glm.r2. <- glmrob(Diversity ~ ., family = "poisson", data=possumDiv, + tcc = 1.6, weights.on.x = "hat", acc = 1e-15) > ## here iterate till convergence (acc = 10^(-15)) > > (sglm.r2 <- summary(glm.r2.)) Call: glmrob(formula = Diversity ~ ., family = "poisson", data = possumDiv, weights.on.x = "hat", tcc = 1.6, acc = 1e-15) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.898214 0.269307 -3.335 0.000852 *** Shrubs 0.007172 0.022435 0.320 0.749204 Stumps -0.253355 0.288588 -0.878 0.379991 Stags 0.040397 0.011343 3.561 0.000369 *** Bark 0.041110 0.014600 2.816 0.004865 ** Habitat 0.073025 0.038677 1.888 0.059017 . BAcacia 0.017699 0.010741 1.648 0.099399 . eucalyptusdelegatensis -0.028994 0.194215 -0.149 0.881328 eucalyptusnitens 0.149521 0.271649 0.550 0.582030 aspectNW-SE 0.050326 0.191676 0.263 0.792890 aspectSE-SW 0.090987 0.192193 0.473 0.635916 aspectSW-NW -0.512248 0.250764 -2.043 0.041077 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.5665 0.7900 0.8635 0.8431 0.9034 0.9411 Number of observations: 151 Fitted by method 'Mqle' (in 18 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: acc tcc 1.0e-15 1.6e+00 maxit 50 test.acc "coef" > ## This is too accurate for S.E. (but we have converged to end) > cf2 <- matrix(c(-0.898213938628341, 0.269306882951903, + 0.00717220104127189, 0.0224349606070713, + -0.25335520175528, 0.288588183720387, + 0.0403970350911325, 0.0113429514237665, + 0.0411096703375411, 0.0145996036305452, + 0.0730250489306713, 0.0386771060643486, + 0.0176994176433365, 0.0107414247342375, + -0.0289935051669504,0.194215229266707, + 0.149521144883774, 0.271648514202971, + 0.0503262879663932, 0.191675979065398, + 0.0909870068741749, 0.192192515800464, + -0.512247626309172, 0.250763990619973), 12,2, byrow=TRUE) > assert.EQ(cf2, unname(coef(sglm.r2)[, 1:2]), tol = 1e-9, giveRE=TRUE)#-> show : ~ 1.46e-11 Mean relative difference: 1.430199e-15 > stopifnot(abs(glm.r2.$iter - 18) <= 1) # 18 iterations on 32-bit (2008) > > ## MT estimator -- "quick" examples > > if(!robustbase:::doExtras()) { + cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' + quit() + } Time elapsed: 0.31 0.04 0.34 NA NA > proc.time() user system elapsed 0.31 0.04 0.34