R Under development (unstable) (2024-01-28 r85838 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> suppressMessages(library(cobs))
> 
> data(USArmyRoofs)
> attach(USArmyRoofs)#-> "age" and "fci"
> 
> if(!dev.interactive(orNone=TRUE)) pdf("roof.pdf", width=10)
> 
> ## Compute the quadratic median smoothing B-spline with SIC
> ## chosen lambda
> a50 <- cobs(age,fci,constraint = "decrease",lambda = -1,nknots = 10,
+             degree = 2,pointwise = rbind(c(0,0,100)),
+             trace = 2)# trace > 1 : more tracing

 Searching for optimal lambda. This may take a while.
   While you are waiting, here is something you can consider
   to speed up the process:
       (a) Use a smaller number of knots;
       (b) Set lambda==0 to exclude the penalty term;
       (c) Use a coarser grid by reducing the argument
 	   'lambda.length' from the default value of 25.
fieq=TRUE -> Tnobs = 184, n0 = 29, |ptConstr| = 2

   The algorithm has converged.  You might
   plot() the returned object (which plots 'sic' against 'lambda')
   to see if you have found the global minimum of the information criterion
   so that you can determine if you need to adjust any or all of
   'lambda.lo', 'lambda.hi' and 'lambda.length' and refit the model.

> summary(a50)
COBS smoothing spline (degree = 2) from call:
	 cobs(x = age, y = fci, constraint = "decrease", nknots = 10,     degree = 2, lambda = -1, pointwise = rbind(c(0, 0, 100)),     trace = 2)
{tau=0.5}-quantile;  dimensionality of fit: 5 from {16,11,9,8,6,5,4}
x$knots[1:10]:  0.04998504,  1.06000000,  2.11000000, ... , 15.01001496
lambda = 19.27405, selected via SIC, out of 25 ones.
with 1 pointwise constraints
coef[1:12]: 99.8569569, 98.4177329, 95.9739856, 94.1164468, 93.0604245, ... ,  0.4726201
R^2 = -5.6% ;  empirical tau (over all): 78/153 = 0.5098039 (target tau= 0.5)
> 
> ## Generate Figure 2a (p.22 in online version)
> ## Plot $pp.lambda againt $pp.sic :
> plot(a50$pp.lambda, a50$pp.sic, type = "b", log = "x",
+      xlab = "Lambda", ylab = "SIC") ##<< no lambda in [20, 180] -- bug !?
> lam0 <- a50$pp.lambda[6] ## should be the 2nd smallest one
> abline(v = c(a50$lambda, lam0),lty = 2)
> 
> ## For "testing" (signif() to stay comparable):
> signif(cbind(a50$pp.lambda, a50$pp.sic), 6)
             [,1]    [,2]
 [1,] 1.59089e-03 2.24395
 [2,] 3.11391e-03 2.24395
 [3,] 6.09499e-03 2.24395
 [4,] 1.19300e-02 2.24395
 [5,] 2.33510e-02 2.24395
 [6,] 4.57060e-02 2.24395
 [7,] 8.94622e-02 2.24395
 [8,] 1.75108e-01 2.24424
 [9,] 3.42746e-01 2.24424
[10,] 6.70872e-01 2.24535
[11,] 1.31313e+00 2.18317
[12,] 2.57024e+00 2.15738
[13,] 5.03083e+00 2.14329
[14,] 9.84705e+00 2.11165
[15,] 1.92740e+01 2.09955
[16,] 3.77259e+01 2.11706
[17,] 7.38425e+01 2.10159
[18,] 1.44535e+02 2.10170
[19,] 2.82904e+02 2.10095
[20,] 5.53740e+02 2.12696
[21,] 1.08386e+03 2.12696
[22,] 2.12148e+03 2.12696
[23,] 4.15247e+03 2.12696
[24,] 8.12780e+03 2.12696
[25,] 1.59089e+04 2.12696
> 
> ## Compute the quadratic median smoothing B-spline with lambda
> ## at the the 2nd smallest SIC
> a50.1 <- cobs(age,fci,constraint = "decrease",lambda = lam0, nknots = 10,
+               degree = 2,pointwise = rbind(c(0,0,100)))
> summary(a50.1)
COBS smoothing spline (degree = 2) from call:
	 cobs(x = age, y = fci, constraint = "decrease", nknots = 10,     degree = 2, lambda = lam0, pointwise = rbind(c(0, 0, 100)))
{tau=0.5}-quantile;  dimensionality of fit: 16 from {16}
x$knots[1:10]:  0.04998504,  1.06000000,  2.11000000, ... , 15.01001496
lambda = 0.04570597
with 1 pointwise constraints
coef[1:12]: 99.53956, 95.00000, 95.00000, 95.00000, 95.00000, ... , 76.69617
R^2 = 1.48% ;  empirical tau (over all): 78/153 = 0.5098039 (target tau= 0.5)
> ## Compute the 25th percentile smoothing B-spline
> a25 <- cobs(age,fci,constraint = "decrease",lambda = -1, nknots = 10,
+             tau = .25, pointwise = rbind(c(0,0,100)))

 Searching for optimal lambda. This may take a while.
   While you are waiting, here is something you can consider
   to speed up the process:
       (a) Use a smaller number of knots;
       (b) Set lambda==0 to exclude the penalty term;
       (c) Use a coarser grid by reducing the argument
 	   'lambda.length' from the default value of 25.

   The algorithm has converged.  You might
   plot() the returned object (which plots 'sic' against 'lambda')
   to see if you have found the global minimum of the information criterion
   so that you can determine if you need to adjust any or all of
   'lambda.lo', 'lambda.hi' and 'lambda.length' and refit the model.

> summary(a25)
COBS smoothing spline (degree = 2) from call:
	 cobs(x = age, y = fci, constraint = "decrease", nknots = 10,     tau = 0.25, lambda = -1, pointwise = rbind(c(0, 0, 100)))
{tau=0.25}-quantile;  dimensionality of fit: 5 from {13,12,11,10,8,7,6,5,3}
x$knots[1:10]:  0.04998504,  1.06000000,  2.11000000, ... , 15.01001496
lambda = 73.84247, selected via SIC, out of 25 ones.
with 1 pointwise constraints
coef[1:12]: 99.6189624, 95.7795144, 88.7927299, 82.9207676, 79.1159073, ... ,  0.8113919
empirical tau (over all): 44/153 = 0.2875817 (target tau= 0.25)
> 
> ## Again plot $pp.sic against $pp.lambda
> plot(a25$pp.lambda,a25$pp.sic,type = "l",log = "x")
> abline(v = a25$lambda,lty = 2)
> 
> ## Compute the 75th percentile smoothing B-spline
> a75 <- cobs(age, fci, constraint = "decrease",lambda = -1, nknots = 10,
+             tau = .75, pointwise = rbind(c(0, 0, 100)))

 Searching for optimal lambda. This may take a while.
   While you are waiting, here is something you can consider
   to speed up the process:
       (a) Use a smaller number of knots;
       (b) Set lambda==0 to exclude the penalty term;
       (c) Use a coarser grid by reducing the argument
 	   'lambda.length' from the default value of 25.

 WARNING!  Since the optimal lambda chosen by SIC 
   rests on a flat portion, you might
   plot() the returned object (which plots 'sic' against 'lambda')
   to see if you want to reduce 'lambda.lo' and/or increase 'lambda.ho'

> a75
COBS smoothing spline (degree = 2) from call:
	 cobs(x = age, y = fci, constraint = "decrease", nknots = 10,     tau = 0.75, lambda = -1, pointwise = rbind(c(0, 0, 100)))
{tau=0.75}-quantile;  dimensionality of fit: 70 from {70}
x$knots[1:10]:  0.04998504,  1.06000000,  2.11000000, ... , 15.01001496
lambda = 5.030829, selected via SIC, out of 25 ones.
> zapsmall(a75$coef)
 [1] 100 100 100 100 100 100 100 100 100 100 100   0
> 
> ## We rerun cobs with ... as suggested by warning of scobs (? right ??)
> ## a75 <- cobs(age, fci, constraint = "decrease", lambda = -1, nknots = 10,
> ##             tau = .75,pointwise = rbind(c(0, 0, 100)))
> ## ## still changing tau , ok ?
> ## a75 <- cobs(age, fci, constraint = "decrease", lambda = -1, nknots = 10,
> ##             tau = .75,pointwise = rbind(c(0, 0, 100)))
> ## summary(a75)
> 
> ## Again we plot $pp.sic against $pp.lambda
> plot(a75$pp.lam,a75$pp.sic,type = "l",log = "x")
> ## It seems like the linear fit is really what the data wants
> 
> (pa50   <- predict(a50,  interval = "both"))
                 z      fit    cb.lo     cb.up    ci.lo     ci.up
  [1,]  0.04998516 99.85696 71.02180 128.69211 82.73109 116.98282
  [2,]  0.20109657 99.43170 73.49827 125.36513 84.02923 114.83416
  [3,]  0.35220798 99.01723 75.03391 123.00056 84.77298 113.26149
  [4,]  0.50331939 98.61356 75.96202 121.26510 85.16029 112.06684
  [5,]  0.65443080 98.22068 76.58136 119.86000 85.36859 111.07277
  [6,]  0.80554221 97.83859 77.07493 118.60226 85.50657 110.17061
  [7,]  0.95665362 97.46729 77.45448 117.48011 85.58122 109.35337
  [8,]  1.10776504 97.10679 77.52028 116.69330 85.47391 108.73967
  [9,]  1.25887645 96.75708 77.13096 116.38320 85.10067 108.41349
 [10,]  1.40998786 96.41816 76.61633 116.21998 84.65740 108.17892
 [11,]  1.56109927 96.09003 76.24170 115.93835 84.30165 107.87841
 [12,]  1.71221068 95.77269 76.12812 115.41726 84.10533 107.44006
 [13,]  1.86332209 95.46615 76.26158 114.67072 84.06011 106.87219
 [14,]  2.01443350 95.17040 76.48429 113.85650 84.07228 106.26851
 [15,]  2.16554491 94.88544 76.47657 113.29430 83.95199 105.81889
 [16,]  2.31665632 94.61127 76.13747 113.08507 83.63925 105.58329
 [17,]  2.46776773 94.34789 75.81251 112.88328 83.33930 105.35649
 [18,]  2.61887914 94.09531 75.73439 112.45623 83.19034 105.00029
 [19,]  2.76999056 93.85352 75.98072 111.72632 83.23845 104.46859
 [20,]  2.92110197 93.62252 76.47502 110.77002 83.43822 103.80682
 [21,]  3.07221338 93.40231 76.95365 109.85097 83.63307 103.17155
 [22,]  3.22332479 93.19290 77.13883 109.24697 83.65802 102.72778
 [23,]  3.37443620 92.99428 77.17837 108.81018 83.60084 102.38771
 [24,]  3.52554761 92.80644 77.06394 108.54895 83.45660 102.15629
 [25,]  3.67665902 92.62499 76.50354 108.74644 83.05009 102.19989
 [26,]  3.82777043 92.43415 75.43346 109.43484 82.33704 102.53125
 [27,]  3.97888184 92.23251 74.16977 110.29524 81.50463 102.96039
 [28,]  4.12999325 92.02008 72.94041 111.09975 80.68822 103.35194
 [29,]  4.28110467 91.79686 71.88118 111.71254 79.96847 103.62524
 [30,]  4.43221608 91.56284 71.06304 112.06265 79.38754 103.73815
 [31,]  4.58332749 91.31804 70.51142 112.12465 78.96051 103.67557
 [32,]  4.73443890 91.06244 70.21552 111.90936 78.68097 103.44391
 [33,]  4.88555031 90.79605 70.12967 111.46243 78.52181 103.07029
 [34,]  5.03666172 90.51887 70.16863 110.86911 78.43239 102.60535
 [35,]  5.18777313 90.23090 70.19903 110.26276 78.33351 102.12828
 [36,]  5.33888454 89.93586 70.01785 109.85388 78.10609 101.76564
 [37,]  5.48999595 89.64979 69.48409 109.81549 77.67292 101.62667
 [38,]  5.64110736 89.37451 68.70505 110.04398 77.09844 101.65059
 [39,]  5.79221877 89.11003 67.79542 110.42464 76.45079 101.76927
 [40,]  5.94333019 88.85633 66.85058 110.86209 75.78661 101.92606
 [41,]  6.09444160 88.61343 65.94497 111.28189 75.15011 102.07676
 [42,]  6.24555301 88.38132 65.13462 111.62803 74.57457 102.18808
 [43,]  6.39666442 88.16001 64.46079 111.85922 74.08449 102.23552
 [44,]  6.54777583 87.94948 63.95348 111.94548 73.69770 102.20126
 [45,]  6.69888724 87.74975 63.63413 111.86536 73.42693 102.07257
 [46,]  6.84999865 87.56081 63.51763 111.60398 73.28101 101.84060
 [47,]  7.00111006 87.38266 63.61358 111.15173 73.26565 101.49966
 [48,]  7.15222147 87.21530 63.92703 110.50356 73.38386 101.04674
 [49,]  7.30333288 87.05874 64.45867 109.65880 73.63603 100.48144
 [50,]  7.45444429 86.91296 65.20438 108.62154 74.01973  99.80619
 [51,]  7.60555571 86.77798 66.15405 107.40191 74.52895  99.02701
 [52,]  7.75666712 86.65379 67.28915 106.01844 75.15268  98.15491
 [53,]  7.90777853 86.54040 68.57837 104.50242 75.87233  97.20846
 [54,]  8.05888994 86.43779 69.96982 102.90576 76.65708  96.21850
 [55,]  8.21000135 86.34598 71.37765 101.31431 77.45594  95.23602
 [56,]  8.36111276 86.26496 72.66142  99.86851 78.18550  94.34442
 [57,]  8.51222417 86.19473 73.60378  98.78569 78.71667  93.67279
 [58,]  8.66333558 86.13530 73.94727  98.32332 78.89655  93.37405
 [59,]  8.81444699 86.08665 73.82563  98.34767 78.80455  93.36876
 [60,]  8.96555840 86.04880 73.67561  98.42200 78.70008  93.39753
 [61,]  9.11666981 86.02174 73.71872  98.32476 78.71469  93.32879
 [62,]  9.26778123 86.00548 73.95881  98.05214 78.85068  93.16027
 [63,]  9.41889264 86.00000 74.17580  97.82420 78.97733  93.02267
 [64,]  9.57000405 86.00000 73.87505  98.12495 78.79871  93.20129
 [65,]  9.72111546 86.00000 72.95882  99.04118 78.25454  93.74546
 [66,]  9.87222687 86.00000 71.64259 100.35741 77.47280  94.52720
 [67,] 10.02333828 86.00000 70.10879 101.89121 76.56184  95.43816
 [68,] 10.17444969 86.00000 68.48528 103.51472 75.59760  96.40240
 [69,] 10.32556110 86.00000 66.85512 105.14488 74.62941  97.37059
 [70,] 10.47667251 86.00000 65.27131 106.72869 73.68875  98.31125
 [71,] 10.62778392 86.00000 63.76791 108.23209 72.79584  99.20416
 [72,] 10.77889533 86.00000 62.36714 109.63286 71.96390 100.03610
 [73,] 10.93000675 86.00000 61.08376 110.91624 71.20167 100.79833
 [74,] 11.08111816 86.00000 59.92764 112.07236 70.51502 101.48498
 [75,] 11.23222957 86.00000 58.90536 113.09464 69.90786 102.09214
 [76,] 11.38334098 86.00000 58.02120 113.97880 69.38274 102.61726
 [77,] 11.53445239 86.00000 57.27775 114.72225 68.94119 103.05881
 [78,] 11.68556380 86.00000 56.67629 115.32371 68.58397 103.41603
 [79,] 11.83667521 86.00000 56.21700 115.78300 68.31119 103.68881
 [80,] 11.98778662 86.00000 55.89910 116.10090 68.12238 103.87762
 [81,] 12.13889803 86.00000 55.72086 116.27914 68.01652 103.98348
 [82,] 12.29000944 86.00000 55.67959 116.32041 67.99201 104.00799
 [83,] 12.44112086 86.00000 55.77155 116.22845 68.04662 103.95338
 [84,] 12.59223227 86.00000 55.99177 116.00823 68.17741 103.82259
 [85,] 12.74334368 86.00000 56.33381 115.66619 68.38056 103.61944
 [86,] 12.89445509 86.00000 56.78946 115.21054 68.65118 103.34882
 [87,] 13.04556650 86.00000 57.34829 114.65171 68.98308 103.01692
 [88,] 13.19667791 86.00000 57.99703 114.00297 69.36839 102.63161
 [89,] 13.34778932 86.00000 58.71888 113.28112 69.79711 102.20289
 [90,] 13.49890073 86.00000 59.49252 112.50748 70.25659 101.74341
 [91,] 13.65001214 86.00000 60.29101 111.70899 70.73083 101.26917
 [92,] 13.80112355 86.00000 61.08052 110.91948 71.19974 100.80026
 [93,] 13.95223496 86.00000 61.81913 110.18087 71.63842 100.36158
 [94,] 14.10334638 86.00000 62.45607 109.54393 72.01671  99.98329
 [95,] 14.25445779 86.00000 62.93209 109.06791 72.29944  99.70056
 [96,] 14.40556920 86.00000 63.18182 108.81818 72.44775  99.55225
 [97,] 14.55668061 86.00000 63.13869 108.86131 72.42214  99.57786
 [98,] 14.70779202 86.00000 62.74238 109.25762 72.18676  99.81324
 [99,] 14.85890343 86.00000 61.94676 110.05324 71.71422 100.28578
[100,] 15.01001484 86.00000 60.72548 111.27452 70.98887 101.01113
> (pa50.1 <- predict(a50.1,interval = "both"))
                 z      fit    cb.lo     cb.up    ci.lo     ci.up
  [1,]  0.04998516 99.53956 61.42903 137.65009 84.84354 114.23558
  [2,]  0.20109657 98.28282 64.00741 132.55823 85.06568 111.49995
  [3,]  0.35220798 97.22931 65.53129 128.92733 85.00605 109.45256
  [4,]  0.50331939 96.37902 66.44118 126.31686 84.83452 107.92352
  [5,]  0.65443080 95.73196 67.13194 124.33198 84.70335 106.76058
  [6,]  0.80554221 95.28813 67.84544 122.73083 84.70580 105.87046
  [7,]  0.95665362 95.04753 68.59721 121.49784 84.84787 105.24718
  [8,]  1.10776504 95.00000 69.11311 120.88689 85.01761 104.98239
  [9,]  1.25887645 95.00000 69.06076 120.93924 84.99742 105.00258
 [10,]  1.40998786 95.00000 68.82854 121.17146 84.90788 105.09212
 [11,]  1.56109927 95.00000 68.76708 121.23292 84.88418 105.11582
 [12,]  1.71221068 95.00000 69.03638 120.96362 84.98802 105.01198
 [13,]  1.86332209 95.00000 69.61791 120.38209 85.21227 104.78773
 [14,]  2.01443350 95.00000 70.30315 119.69685 85.47651 104.52349
 [15,]  2.16554491 95.00000 70.66957 119.33043 85.61781 104.38219
 [16,]  2.31665632 95.00000 70.58375 119.41625 85.58471 104.41529
 [17,]  2.46776773 95.00000 70.50236 119.49764 85.55333 104.44667
 [18,]  2.61887914 95.00000 70.73294 119.26706 85.64224 104.35776
 [19,]  2.76999056 95.00000 71.37807 118.62193 85.89101 104.10899
 [20,]  2.92110197 95.00000 72.33668 117.66332 86.26067 103.73933
 [21,]  3.07221338 95.00000 73.26031 116.73969 86.61684 103.38316
 [22,]  3.22332479 95.00000 73.78183 116.21817 86.81794 103.18206
 [23,]  3.37443620 95.00000 74.09661 115.90339 86.93932 103.06068
 [24,]  3.52554761 95.00000 74.19361 115.80639 86.97673 103.02327
 [25,]  3.67665902 95.00000 73.69278 116.30722 86.78360 103.21640
 [26,]  3.82777043 95.00000 72.53071 117.46929 86.33549 103.66451
 [27,]  3.97888184 95.00000 71.12704 118.87296 85.79421 104.20579
 [28,]  4.12999325 95.00000 69.78299 120.21701 85.27593 104.72407
 [29,]  4.28110467 95.00000 68.67806 121.32194 84.84985 105.15015
 [30,]  4.43221608 95.00000 67.90605 122.09395 84.55215 105.44785
 [31,]  4.58332749 95.00000 67.50054 122.49946 84.39578 105.60422
 [32,]  4.73443890 95.00000 67.44726 122.55274 84.37523 105.62477
 [33,]  4.88555031 95.00000 67.68588 122.31412 84.46725 105.53275
 [34,]  5.03666172 95.00000 68.10372 121.89628 84.62837 105.37163
 [35,]  5.18777313 95.00000 68.52450 121.47550 84.79063 105.20937
 [36,]  5.33888454 94.99079 68.66576 121.31582 84.83945 105.14214
 [37,]  5.48999595 94.93286 68.28048 121.58525 84.65529 105.21044
 [38,]  5.64110736 94.82170 67.50351 122.13990 84.28738 105.35603
 [39,]  5.79221877 94.65731 66.48645 122.82817 83.79419 105.52044
 [40,]  5.94333019 94.43969 65.35536 123.52402 83.22432 105.65506
 [41,]  6.09444160 94.16884 64.20863 124.12904 82.61571 105.72196
 [42,]  6.24555301 93.84475 63.12029 124.56920 81.99692 105.69258
 [43,]  6.39666442 93.46743 62.14491 124.78995 81.38897 105.54589
 [44,]  6.54777583 93.03688 61.32210 124.75165 80.80716 105.26659
 [45,]  6.69888724 92.55310 60.68023 124.42596 80.26242 104.84377
 [46,]  6.84999865 92.01608 60.23896 123.79320 79.76233 104.26984
 [47,]  7.00111006 91.42583 60.01098 122.84069 79.31177 103.53989
 [48,]  7.15222147 90.78236 60.00297 121.56174 78.91334 102.65137
 [49,]  7.30333288 90.08565 60.21584 119.95545 78.56738 101.60391
 [50,]  7.45444429 89.33570 60.64415 118.02726 78.27179 100.39962
 [51,]  7.60555571 88.53253 61.27451 115.79054 78.02141  99.04365
 [52,]  7.75666712 87.67612 62.08247 113.26978 77.80681  97.54544
 [53,]  7.90777853 86.76648 63.02663 110.50634 77.61202  95.92094
 [54,]  8.05888994 85.80361 64.03841 107.56882 77.41061  94.19662
 [55,]  8.21000135 84.78751 65.00434 104.57069 77.15881  92.41621
 [56,]  8.36111276 83.71818 65.73879 101.69756 76.78504  90.65131
 [57,]  8.51222417 82.59561 65.95453  99.23669 76.17855  89.01267
 [58,]  8.66333558 82.00000 65.89146  98.10854 75.78830  88.21170
 [59,]  8.81444699 82.00000 65.79499  98.20501 75.75109  88.24891
 [60,]  8.96555840 82.00000 65.64673  98.35327 75.69392  88.30608
 [61,]  9.11666981 82.00000 65.73948  98.26052 75.72969  88.27031
 [62,]  9.26778123 82.00000 66.07830  97.92170 75.86034  88.13966
 [63,]  9.41889264 82.00000 66.37232  97.62768 75.97372  88.02628
 [64,]  9.57000405 82.00000 65.97483  98.02517 75.82044  88.17956
 [65,]  9.72111546 82.00000 64.76387  99.23613 75.35348  88.64652
 [66,]  9.87222687 82.00000 63.02426 100.97574 74.68266  89.31734
 [67,] 10.02333828 82.00000 60.99708 103.00292 73.90095  90.09905
 [68,] 10.17444969 82.00000 58.85133 105.14867 73.07351  90.92649
 [69,] 10.32556110 82.00000 56.69680 107.30320 72.24269  91.75731
 [70,] 10.47667251 82.00000 54.60353 109.39647 71.43549  92.56451
 [71,] 10.62778392 82.00000 52.61653 111.38347 70.66927  93.33073
 [72,] 10.77889533 82.00000 50.76518 113.23482 69.95536  94.04464
 [73,] 10.93000675 82.00000 49.06898 114.93102 69.30128  94.69872
 [74,] 11.08111816 82.00000 47.54097 116.45903 68.71206  95.28794
 [75,] 11.23222957 82.00000 46.18985 117.81015 68.19104  95.80896
 [76,] 11.38334098 82.00000 45.02128 118.97872 67.74043  96.25957
 [77,] 11.53445239 82.00000 44.03869 119.96131 67.36152  96.63848
 [78,] 11.68556380 82.00000 43.24375 120.75625 67.05498  96.94502
 [79,] 11.83667521 82.00000 42.63673 121.36327 66.82091  97.17909
 [80,] 11.98778662 82.00000 42.21657 121.78343 66.65888  97.34112
 [81,] 12.13889803 82.00000 41.98099 122.01901 66.56804  97.43196
 [82,] 12.29000944 82.00000 41.92645 122.07355 66.54701  97.45299
 [83,] 12.44112086 82.00000 42.04799 121.95201 66.59388  97.40612
 [84,] 12.59223227 82.00000 42.33904 121.66096 66.70611  97.29389
 [85,] 12.74334368 82.00000 42.79111 121.20889 66.88044  97.11956
 [86,] 12.89445509 82.00000 43.39333 120.60667 67.11267  96.88733
 [87,] 13.04556650 82.00000 44.13191 119.86809 67.39747  96.60253
 [88,] 13.19667791 82.00000 44.98933 119.01067 67.72811  96.27189
 [89,] 13.34778932 82.00000 45.94338 118.05662 68.09600  95.90400
 [90,] 13.49890073 82.00000 46.96588 117.03412 68.49029  95.50971
 [91,] 13.65001214 82.00000 48.02122 115.97878 68.89725  95.10275
 [92,] 13.80112355 82.00000 49.06469 114.93531 69.29963  94.70037
 [93,] 13.95223496 82.00000 50.04089 113.95911 69.67607  94.32393
 [94,] 14.10334638 82.00000 50.88271 113.11729 70.00069  93.99931
 [95,] 14.25445779 82.00000 51.51186 112.48814 70.24330  93.75670
 [96,] 14.40556920 82.00000 51.84191 112.15809 70.37057  93.62943
 [97,] 14.55668061 82.00000 51.78491 112.21509 70.34859  93.65141
 [98,] 14.70779202 82.00000 51.26113 112.73887 70.14661  93.85339
 [99,] 14.85890343 82.00000 50.20957 113.79043 69.74111  94.25889
[100,] 15.01001484 82.00000 48.59544 115.40456 69.11868  94.88132
> (pa25   <- predict(a25,  interval = "both"))
                 z      fit    cb.lo     cb.up    ci.lo     ci.up
  [1,]  0.04998516 99.61896 71.09477 128.14315 82.67778 116.56014
  [2,]  0.20109657 98.47936 72.82560 124.13312 83.24300 113.71572
  [3,]  0.35220798 97.35829 73.63361 121.08298 83.26765 111.44893
  [4,]  0.50331939 96.25575 73.84849 118.66301 82.94756 109.56394
  [5,]  0.65443080 95.17173 73.76577 116.57769 82.45824 107.88523
  [6,]  0.80554221 94.10624 73.56650 114.64599 81.90721 106.30528
  [7,]  0.95665362 93.05929 73.26229 112.85628 81.30139 104.81718
  [8,]  1.10776504 92.03085 72.65556 111.40614 80.52342 103.53829
  [9,]  1.25887645 91.02095 71.60648 110.43542 79.49025 102.55165
 [10,]  1.40998786 90.02957 70.44130 109.61785 78.39564 101.66351
 [11,]  1.56109927 89.05672 69.42245 108.69100 77.39547 100.71798
 [12,]  1.71221068 88.10240 68.66968 107.53512 76.56086  99.64395
 [13,]  1.86332209 87.16661 68.16915 106.16408 75.88358  98.44965
 [14,]  2.01443350 86.24935 67.76475 104.73394 75.27092  97.22778
 [15,]  2.16554491 85.35061 67.14027 103.56095 74.53507  96.16615
 [16,]  2.31665632 84.47040 66.19583 102.74498 73.61671  95.32409
 [17,]  2.46776773 83.60872 65.27323 101.94421 72.71884  94.49860
 [18,]  2.61887914 82.76557 64.60266 100.92848 71.97819  93.55294
 [19,]  2.76999056 81.94094 64.26088  99.62100 71.44034  92.44154
 [20,]  2.92110197 81.13484 64.17227  98.09742 71.06037  91.20931
 [21,]  3.07221338 80.34727 64.07600  96.61855 70.68338  90.01116
 [22,]  3.22332479 79.57823 63.69729  95.45917 70.14617  89.01029
 [23,]  3.37443620 78.82772 63.18237  94.47306 69.53558  88.11985
 [24,]  3.52554761 78.09573 62.52299  93.66847 68.84672  87.34474
 [25,]  3.67665902 77.38227 61.43468  93.32987 67.91063  86.85392
 [26,]  3.82777043 76.68734 59.86999  93.50470 66.69913  86.67556
 [27,]  3.97888184 76.01094 58.14300  93.87888 65.39875  86.62313
 [28,]  4.12999325 75.35307 56.47916  94.22697 64.14341  86.56272
 [29,]  4.28110467 74.71372 55.01281  94.41462 63.01289  86.41454
 [30,]  4.43221608 74.09290 53.81417  94.37163 62.04889  86.13691
 [31,]  4.58332749 73.49061 52.90837  94.07284 61.26634  85.71488
 [32,]  4.73443890 72.90685 52.28474  93.52895 60.65890  85.15479
 [33,]  4.88555031 72.34161 51.89810  92.78512 60.19973  84.48349
 [34,]  5.03666172 71.79490 51.66412  91.92568 59.83877  83.75104
 [35,]  5.18777313 71.26672 51.45089  91.08256 59.49764  83.03581
 [36,]  5.33888454 70.75707 51.05385  90.46029 59.05487  82.45927
 [37,]  5.48999595 70.26595 50.31772  90.21418 58.41823  82.11366
 [38,]  5.64110736 69.79335 49.34679  90.23991 57.64966  81.93704
 [39,]  5.79221877 69.33928 48.25453  90.42403 56.81656  81.86200
 [40,]  5.94333019 68.90374 47.13530  90.67219 55.97496  81.83253
 [41,]  6.09444160 68.48673 46.06273  90.91073 55.16859  81.80486
 [42,]  6.24555301 68.08824 45.09223  91.08425 54.43038  81.74611
 [43,]  6.39666442 67.70829 44.26465  91.15193 53.78457  81.63201
 [44,]  6.54777583 67.34686 43.60963  91.08408 53.24877  81.44495
 [45,]  6.69888724 67.00396 43.14841  90.85950 52.83559  81.17232
 [46,]  6.84999865 66.67958 42.89570  90.46347 52.55378  80.80539
 [47,]  7.00111006 66.37374 42.86099  89.88649 52.40897  80.33850
 [48,]  7.15222147 66.08642 43.04930  89.12354 52.40414  79.76870
 [49,]  7.30333288 65.81763 43.46129  88.17397 52.53968  79.09558
 [50,]  7.45444429 65.56737 44.09290  87.04184 52.81318  78.32156
 [51,]  7.60555571 65.33564 44.93411  85.73716 53.21870  77.45257
 [52,]  7.75666712 65.12243 45.96662  84.27824 53.74535  76.49951
 [53,]  7.90777853 64.92775 47.15943  82.69607 54.37473  75.48077
 [54,]  8.05888994 64.75160 48.46122  81.04198 55.07637  74.42683
 [55,]  8.21000135 64.59398 49.78707  79.40088 55.79981  73.38814
 [56,]  8.36111276 64.45488 50.99804  77.91173 56.46255  72.44721
 [57,]  8.51222417 64.33432 51.87914  76.78949 56.93690  71.73173
 [58,]  8.66333558 64.23228 52.17569  76.28887 57.07159  71.39297
 [59,]  8.81444699 64.14877 52.01997  76.27756 56.94519  71.35234
 [60,]  8.96555840 64.08378 51.84402  76.32354 56.81431  71.35326
 [61,]  9.11666981 64.03733 51.86698  76.20767 56.80908  71.26558
 [62,]  9.26778123 64.00940 52.09265  75.92615 56.93177  71.08704
 [63,]  9.41889264 64.00000 52.30332  75.69669 57.05307  70.94693
 [64,]  9.57000405 64.00000 52.00581  75.99419 56.87637  71.12363
 [65,]  9.72111546 64.00000 51.09945  76.90055 56.33807  71.66193
 [66,]  9.87222687 64.00000 49.79742  78.20258 55.56476  72.43524
 [67,] 10.02333828 64.00000 48.28016  79.71984 54.66363  73.33637
 [68,] 10.17444969 64.00000 46.67416  81.32584 53.70978  74.29022
 [69,] 10.32556110 64.00000 45.06158  82.93842 52.75203  75.24797
 [70,] 10.47667251 64.00000 43.49485  84.50515 51.82152  76.17848
 [71,] 10.62778392 64.00000 42.00766  85.99234 50.93824  77.06176
 [72,] 10.77889533 64.00000 40.62200  87.37800 50.11526  77.88474
 [73,] 10.93000675 64.00000 39.35246  88.64754 49.36126  78.63874
 [74,] 11.08111816 64.00000 38.20881  89.79119 48.68201  79.31799
 [75,] 11.23222957 64.00000 37.19755  90.80245 48.08140  79.91860
 [76,] 11.38334098 64.00000 36.32292  91.67708 47.56194  80.43806
 [77,] 11.53445239 64.00000 35.58749  92.41251 47.12515  80.87485
 [78,] 11.68556380 64.00000 34.99252  93.00748 46.77178  81.22822
 [79,] 11.83667521 64.00000 34.53818  93.46182 46.50194  81.49806
 [80,] 11.98778662 64.00000 34.22371  93.77629 46.31517  81.68483
 [81,] 12.13889803 64.00000 34.04739  93.95261 46.21045  81.78955
 [82,] 12.29000944 64.00000 34.00657  93.99343 46.18621  81.81379
 [83,] 12.44112086 64.00000 34.09754  93.90246 46.24024  81.75976
 [84,] 12.59223227 64.00000 34.31538  93.68462 46.36962  81.63038
 [85,] 12.74334368 64.00000 34.65373  93.34627 46.57057  81.42943
 [86,] 12.89445509 64.00000 35.10447  92.89553 46.83828  81.16172
 [87,] 13.04556650 64.00000 35.65727  92.34273 47.16660  80.83340
 [88,] 13.19667791 64.00000 36.29902  91.70098 47.54774  80.45226
 [89,] 13.34778932 64.00000 37.01308  90.98692 47.97184  80.02816
 [90,] 13.49890073 64.00000 37.77838  90.22162 48.42637  79.57363
 [91,] 13.65001214 64.00000 38.56826  89.43174 48.89550  79.10450
 [92,] 13.80112355 64.00000 39.34926  88.65074 49.35935  78.64065
 [93,] 13.95223496 64.00000 40.07990  87.92010 49.79330  78.20670
 [94,] 14.10334638 64.00000 40.70997  87.29003 50.16751  77.83249
 [95,] 14.25445779 64.00000 41.18086  86.81914 50.44718  77.55282
 [96,] 14.40556920 64.00000 41.42789  86.57211 50.59390  77.40610
 [97,] 14.55668061 64.00000 41.38523  86.61477 50.56856  77.43144
 [98,] 14.70779202 64.00000 40.99320  87.00680 50.33573  77.66427
 [99,] 14.85890343 64.00000 40.20615  87.79385 49.86828  78.13172
[100,] 15.01001484 64.00000 38.99804  89.00196 49.15076  78.84924
> (pa75   <- predict(a75,  interval = "both"))
                 z fit cb.lo cb.up ci.lo ci.up
  [1,]  0.04998516 100   100   100   100   100
  [2,]  0.20109657 100   100   100   100   100
  [3,]  0.35220798 100   100   100   100   100
  [4,]  0.50331939 100   100   100   100   100
  [5,]  0.65443080 100   100   100   100   100
  [6,]  0.80554221 100   100   100   100   100
  [7,]  0.95665362 100   100   100   100   100
  [8,]  1.10776504 100   100   100   100   100
  [9,]  1.25887645 100   100   100   100   100
 [10,]  1.40998786 100   100   100   100   100
 [11,]  1.56109927 100   100   100   100   100
 [12,]  1.71221068 100   100   100   100   100
 [13,]  1.86332209 100   100   100   100   100
 [14,]  2.01443350 100   100   100   100   100
 [15,]  2.16554491 100   100   100   100   100
 [16,]  2.31665632 100   100   100   100   100
 [17,]  2.46776773 100   100   100   100   100
 [18,]  2.61887914 100   100   100   100   100
 [19,]  2.76999056 100   100   100   100   100
 [20,]  2.92110197 100   100   100   100   100
 [21,]  3.07221338 100   100   100   100   100
 [22,]  3.22332479 100   100   100   100   100
 [23,]  3.37443620 100   100   100   100   100
 [24,]  3.52554761 100   100   100   100   100
 [25,]  3.67665902 100   100   100   100   100
 [26,]  3.82777043 100   100   100   100   100
 [27,]  3.97888184 100   100   100   100   100
 [28,]  4.12999325 100   100   100   100   100
 [29,]  4.28110467 100   100   100   100   100
 [30,]  4.43221608 100   100   100   100   100
 [31,]  4.58332749 100   100   100   100   100
 [32,]  4.73443890 100   100   100   100   100
 [33,]  4.88555031 100   100   100   100   100
 [34,]  5.03666172 100   100   100   100   100
 [35,]  5.18777313 100   100   100   100   100
 [36,]  5.33888454 100   100   100   100   100
 [37,]  5.48999595 100   100   100   100   100
 [38,]  5.64110736 100   100   100   100   100
 [39,]  5.79221877 100   100   100   100   100
 [40,]  5.94333019 100   100   100   100   100
 [41,]  6.09444160 100   100   100   100   100
 [42,]  6.24555301 100   100   100   100   100
 [43,]  6.39666442 100   100   100   100   100
 [44,]  6.54777583 100   100   100   100   100
 [45,]  6.69888724 100   100   100   100   100
 [46,]  6.84999865 100   100   100   100   100
 [47,]  7.00111006 100   100   100   100   100
 [48,]  7.15222147 100   100   100   100   100
 [49,]  7.30333288 100   100   100   100   100
 [50,]  7.45444429 100   100   100   100   100
 [51,]  7.60555571 100   100   100   100   100
 [52,]  7.75666712 100   100   100   100   100
 [53,]  7.90777853 100   100   100   100   100
 [54,]  8.05888994 100   100   100   100   100
 [55,]  8.21000135 100   100   100   100   100
 [56,]  8.36111276 100   100   100   100   100
 [57,]  8.51222417 100   100   100   100   100
 [58,]  8.66333558 100   100   100   100   100
 [59,]  8.81444699 100   100   100   100   100
 [60,]  8.96555840 100   100   100   100   100
 [61,]  9.11666981 100   100   100   100   100
 [62,]  9.26778123 100   100   100   100   100
 [63,]  9.41889264 100   100   100   100   100
 [64,]  9.57000405 100   100   100   100   100
 [65,]  9.72111546 100   100   100   100   100
 [66,]  9.87222687 100   100   100   100   100
 [67,] 10.02333828 100   100   100   100   100
 [68,] 10.17444969 100   100   100   100   100
 [69,] 10.32556110 100   100   100   100   100
 [70,] 10.47667251 100   100   100   100   100
 [71,] 10.62778392 100   100   100   100   100
 [72,] 10.77889533 100   100   100   100   100
 [73,] 10.93000675 100   100   100   100   100
 [74,] 11.08111816 100   100   100   100   100
 [75,] 11.23222957 100   100   100   100   100
 [76,] 11.38334098 100   100   100   100   100
 [77,] 11.53445239 100   100   100   100   100
 [78,] 11.68556380 100   100   100   100   100
 [79,] 11.83667521 100   100   100   100   100
 [80,] 11.98778662 100   100   100   100   100
 [81,] 12.13889803 100   100   100   100   100
 [82,] 12.29000944 100   100   100   100   100
 [83,] 12.44112086 100   100   100   100   100
 [84,] 12.59223227 100   100   100   100   100
 [85,] 12.74334368 100   100   100   100   100
 [86,] 12.89445509 100   100   100   100   100
 [87,] 13.04556650 100   100   100   100   100
 [88,] 13.19667791 100   100   100   100   100
 [89,] 13.34778932 100   100   100   100   100
 [90,] 13.49890073 100   100   100   100   100
 [91,] 13.65001214 100   100   100   100   100
 [92,] 13.80112355 100   100   100   100   100
 [93,] 13.95223496 100   100   100   100   100
 [94,] 14.10334638 100   100   100   100   100
 [95,] 14.25445779 100   100   100   100   100
 [96,] 14.40556920 100   100   100   100   100
 [97,] 14.55668061 100   100   100   100   100
 [98,] 14.70779202 100   100   100   100   100
 [99,] 14.85890343 100   100   100   100   100
[100,] 15.01001484 100   100   100   100   100
> ## Generate Figure 2b (p.22, online version)
> plot(age,fci,xlim = c(0,15),ylim = c(0,100),xlab = "AGE",ylab = "FCI")
> lines(pa50)
> lines(pa50.1,lty = 2)
> lines(pa25, col = 2)
> lines(pa75, col = 3)
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  3.07 0.21 3.28 NA NA 
> 
> proc.time()
   user  system elapsed 
   3.07    0.21    3.28