roof.Rout-N-save
R : Copyright 2002, The R Development Core Team
Version 1.5.0 Under development (unstable) (2002-04-17)
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.
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.
> library(cobs)
> data(USArmyRoofs)
> attach(USArmyRoofs)#-> "age" and "fci"
>
> if(!interactive()) postscript("roof.ps", horizontal = TRUE)
>
> ## 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)),
+ lstart = 1e3, trace = 2)# trace > 1 : more tracing
You are using 10 knots instead of the default number of 20 knots.
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) Increase `factor' to 2 or 3;
(c) Set lambda==0 to exclude the penalty term.
dcrqL1lt(): N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12)
LAMBDA will be updated; (tau,lam)= (0.5, -1):
after setup() & newpen(): ifl = 0 (a|cg)mag = (1000, 58.254)
search used icyc = 35 iter.; ifl =1; NLT: LAMBDA = 507.887
after setup() & newpen(): ifl = 0 (a|cg)mag = (507.887, 58.254)
search used icyc = 19 iter.; ifl =1; NLT: LAMBDA = 498.878
after setup() & newpen(): ifl = 0 (a|cg)mag = (498.878, 58.254)
search used icyc = 0 iter.; ifl =1; NLT: LAMBDA = 485.665
after setup() & newpen(): ifl = 0 (a|cg)mag = (485.665, 58.254)
search used icyc = 3 iter.; ifl =1; NLT: LAMBDA = 446.085
after setup() & newpen(): ifl = 0 (a|cg)mag = (446.085, 58.254)
search used icyc = 4 iter.; ifl =1; NLT: LAMBDA = 239.003
after setup() & newpen(): ifl = 0 (a|cg)mag = (239.003, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 223.339
after setup() & newpen(): ifl = 0 (a|cg)mag = (223.339, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 218.431
after setup() & newpen(): ifl = 0 (a|cg)mag = (218.431, 58.254)
search used icyc = 8 iter.; ifl =1; NLT: LAMBDA = 209.02
after setup() & newpen(): ifl = 0 (a|cg)mag = (209.02, 58.254)
search used icyc = 12 iter.; ifl =1; NLT: LAMBDA = 189.155
after setup() & newpen(): ifl = 0 (a|cg)mag = (189.155, 58.254)
search used icyc = 4 iter.; ifl =1; NLT: LAMBDA = 19.5685
after setup() & newpen(): ifl = 0 (a|cg)mag = (19.5685, 58.254)
search used icyc = 3 iter.; ifl =1; NLT: LAMBDA = 14.4345
after setup() & newpen(): ifl = 0 (a|cg)mag = (14.4345, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 14.4078
after setup() & newpen(): ifl = 0 (a|cg)mag = (14.4078, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 14.0291
after setup() & newpen(): ifl = 0 (a|cg)mag = (14.0291, 58.254)
search used icyc = 1 iter.; ifl =1; NLT: LAMBDA = 9.21145
after setup() & newpen(): ifl = 0 (a|cg)mag = (9.21145, 58.254)
search used icyc = 1 iter.; ifl =1; NLT: LAMBDA = 7.84163
after setup() & newpen(): ifl = 0 (a|cg)mag = (7.84163, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 4.3684
after setup() & newpen(): ifl = 0 (a|cg)mag = (4.3684, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 2.44346
after setup() & newpen(): ifl = 0 (a|cg)mag = (2.44346, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 1.73756
after setup() & newpen(): ifl = 0 (a|cg)mag = (1.73756, 58.254)
search used icyc = 2 iter.; ifl =1; NLT: LAMBDA = 1.59067
after setup() & newpen(): ifl = 0 (a|cg)mag = (1.59067, 58.254)
search used icyc = 1 iter.; ifl =1; NLT: LAMBDA = 1.22009
> 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, lstart = 1000)
{tau=0.5}-quantile; dimensionality of fit: 4 (4)
knots[1 .. 10]: 0.049999, 1.060000, 2.110000, ... , 15.010001
lambda = 19.58811, selected via SIC, out of 21 ones.
with 1 pointwise constraints
coef :
[1] 99.8569171 98.4177163 95.9739856 94.1164468 93.0604245 91.7548114
[7] 86.8671634 86.1622032 86.0000000 86.0000000 86.0000000 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 $sic :
> plot(a50$pp.lambda, a50$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$sic), 6)
[,1] [,2]
[1,] 1000.00000 2.78723
[2,] 508.39600 2.78432
[3,] 499.37800 2.78432
[4,] 486.15100 2.77828
[5,] 446.53200 2.77828
[6,] 239.24200 2.77766
[7,] 223.56300 2.79138
[8,] 218.65000 2.79124
[9,] 209.22900 2.78075
[10,] 189.34500 2.78075
[11,] 19.58810 2.77626
[12,] 14.44890 2.78858
[13,] 14.42220 2.78856
[14,] 14.04310 2.78836
[15,] 9.22067 2.80359
[16,] 7.84948 2.82000
[17,] 4.37277 2.83408
[18,] 2.44591 2.83265
[19,] 1.73930 2.84903
[20,] 1.59226 2.86466
[21,] 1.22009 2.86466
>
> ## 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)), lstart = 1e3)
You are using 10 knots instead of the default number of 20 knots.
N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12) _fixed_ (tau,lam)= (0.5, 1):
> 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)), lstart = 1000)
{tau=0.5}-quantile; dimensionality of fit: 3 (3)
knots[1 .. 10]: 0.049999, 1.060000, 2.110000, ... , 15.010001
lambda = 239.2424
with 1 pointwise constraints
coef :
[1] 99.9071005 98.9703611 97.1887743 95.6052663 94.5143866 92.9483429
[7] 89.9471296 88.5873474 88.2398212 86.0322399 86.0322399 0.1239925
R^2 = -13.39% ; empirical tau (over all): 82/153 = 0.5359477 (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)))
You are using 10 knots instead of the default number of 20 knots.
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) Increase `factor' to 2 or 3;
(c) Set lambda==0 to exclude the penalty term.
N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12)
LAMBDA will be updated; (tau,lam)= (0.25, -1):
> 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: 4 (4)
knots[1 .. 10]: 0.049999, 1.060000, 2.110000, ... , 15.010001
lambda = 57.62234, selected via SIC, out of 22 ones.
with 1 pointwise constraints
coef :
[1] 99.618856 95.779467 88.792730 82.920768 79.115907 73.879833 65.488742
[8] 64.278470 64.000000 64.000000 64.000000 0.811392
empirical tau (over all): 44/153 = 0.2875817 (target tau : 0.25)
>
> ## Again plot $sic against $pp.lambda
> plot(a25$pp.lambda,a25$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)))
You are using 10 knots instead of the default number of 20 knots.
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) Increase `factor' to 2 or 3;
(c) Set lambda==0 to exclude the penalty term.
N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12)
LAMBDA will be updated; (tau,lam)= (0.75, -1):
WARNING! The algorithm has not converged after 3060 iterations.
Increase the `maxiter' counter and restart cobs with both
`coef' and `knots' set to the values at the last iteration.
> summary(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: 10 (10)
knots[1 .. 10]: 0.049999, 1.060000, 2.110000, ... , 15.010001
lambda = 7872, selected via SIC, out of 2 ones.
with 1 pointwise constraints
coef :
[1] 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02
[6] 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02
[11] 1.000000e+02 1.776357e-13
empirical tau (over all): 124/153 = 0.8104575 (target tau : 0.75)
> ## We rerun cobs with a larger lstart at 10^8 as suggested by warning of cobs
> a75 <- cobs(age, fci, constraint = "decrease",lambda = -1, nknots = 10,
+ tau = .75,lstart = 10^8,pointwise = rbind(c(0, 0, 100)))
You are using 10 knots instead of the default number of 20 knots.
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) Increase `factor' to 2 or 3;
(c) Set lambda==0 to exclude the penalty term.
N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12)
LAMBDA will be updated; (tau,lam)= (0.75, -1):
WARNING! Since the optimal lambda chosen by SIC reached the smoothest
possible fit allowed by `lstart', you might want to rerun cobs with
a larger `lstart' value to see if it makes a difference if you haven't
done so.
> ## still lstart...
> a75 <- cobs(age, fci, constraint = "decrease",lambda = -1, nknots = 10,
+ tau = .75,lstart = 10^10,pointwise = rbind(c(0, 0, 100)))
You are using 10 knots instead of the default number of 20 knots.
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) Increase `factor' to 2 or 3;
(c) Set lambda==0 to exclude the penalty term.
N[rq,L1]= (153, 1); ([ei]qc= 1,28, vars=12)
LAMBDA will be updated; (tau,lam)= (0.75, -1):
WARNING! Since the optimal lambda chosen by SIC reached the smoothest
possible fit allowed by `lstart', you might want to rerun cobs with
a larger `lstart' value to see if it makes a difference if you haven't
done so.
> summary(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)), lstart = 10^10)
{tau=0.75}-quantile; dimensionality of fit: 2 (2)
knots[1 .. 10]: 0.049999, 1.060000, 2.110000, ... , 15.010001
lambda = 1e+10, selected via SIC, out of 7 ones.
with 1 pointwise constraints
coef :
[1] 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02
[6] 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02 1.000000e+02
[11] 1.000000e+02 1.300979e-29
empirical tau (over all): 100/153 = 0.6535948 (target tau : 0.75)
>
> ## Again we plot $sic against $pp.lambda
> plot(a75$pp.lam,a75$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.04999912 99.85692 71.53034 128.18349 81.68492 118.02892
[2,] 0.20111025 99.43166 74.04410 124.81922 83.14509 115.71823
[3,] 0.35222138 99.01720 75.70451 122.32988 84.06170 113.97270
[4,] 0.50333251 98.61353 76.79243 120.43463 84.61490 112.61215
[5,] 0.65444363 98.22065 77.56410 118.87719 84.96911 111.47219
[6,] 0.80555476 97.83856 78.18084 117.49628 85.22779 110.44933
[7,] 0.95666589 97.46726 78.65431 116.28022 85.39842 109.53611
[8,] 1.10777702 97.10676 78.80342 115.41010 85.36484 108.84868
[9,] 1.25888815 96.75705 78.49749 115.01662 85.04321 108.47089
[10,] 1.40999928 96.41813 78.04503 114.79124 84.63146 108.20481
[11,] 1.56111041 96.09001 77.70011 114.47990 84.29256 107.88745
[12,] 1.71222154 95.77267 77.58307 113.96227 84.10372 107.44162
[13,] 1.86333267 95.46613 77.68726 113.24500 84.06067 106.87159
[14,] 2.01444379 95.17038 77.87110 112.46966 84.07258 106.26817
[15,] 2.16555492 94.88542 77.84165 111.92918 83.95154 105.81930
[16,] 2.31666605 94.61125 77.50632 111.71618 83.63814 105.58437
[17,] 2.46777718 94.34788 77.18511 111.51065 83.33766 105.35810
[18,] 2.61888831 94.09530 77.09354 111.09705 83.18837 105.00223
[19,] 2.76999944 93.85351 77.30350 110.40351 83.23639 104.47063
[20,] 2.92111057 93.62251 77.74426 109.50075 83.43633 103.80868
[21,] 3.07222170 93.40230 78.17174 108.63286 83.63163 103.17298
[22,] 3.22333283 93.19289 78.32845 108.05733 83.65708 102.72869
[23,] 3.37444396 92.99427 78.35101 107.63753 83.60035 102.38818
[24,] 3.52555508 92.80644 78.23168 107.38119 83.45647 102.15640
[25,] 3.67666621 92.62498 77.69984 107.55013 83.05024 102.19973
[26,] 3.82777734 92.43414 76.69530 108.17297 82.33740 102.53088
[27,] 3.97888847 92.23250 75.51063 108.95437 81.50512 102.95987
[28,] 4.12999960 92.02007 74.35685 109.68329 80.68880 103.35134
[29,] 4.28111073 91.79685 73.35972 110.23398 79.96910 103.62459
[30,] 4.43222186 91.56284 72.58495 110.54072 79.38819 103.73748
[31,] 4.58333299 91.31803 72.05610 110.57996 78.96116 103.67490
[32,] 4.73444412 91.06243 71.76318 110.36169 78.68162 103.44325
[33,] 4.88555524 90.79604 71.66391 109.92817 78.52244 103.06964
[34,] 5.03666637 90.51886 71.67939 109.35833 78.43301 102.60471
[35,] 5.18777750 90.23089 71.68616 108.77561 78.33411 102.12766
[36,] 5.33888863 89.93585 71.49653 108.37518 78.10670 101.76501
[37,] 5.48999976 89.64979 70.98117 108.31840 77.67354 101.62603
[38,] 5.64111089 89.37451 70.23953 108.50948 77.09908 101.64993
[39,] 5.79222202 89.11002 69.37781 108.84224 76.45146 101.76859
[40,] 5.94333315 88.85633 68.48429 109.22837 75.78730 101.92536
[41,] 6.09444428 88.61343 67.62789 109.59897 75.15083 102.07602
[42,] 6.24555540 88.38132 66.86047 109.90216 74.57531 102.18732
[43,] 6.39666653 88.16000 66.22025 110.09976 74.08526 102.23475
[44,] 6.54777766 87.94948 65.73498 110.16398 73.69848 102.20048
[45,] 6.69888879 87.74975 65.42452 110.07497 73.42771 102.07178
[46,] 6.84999992 87.56080 65.30264 109.81897 73.28180 101.83981
[47,] 7.00111105 87.38266 65.37824 109.38707 73.26643 101.49888
[48,] 7.15222218 87.21530 65.65600 108.77460 73.38462 101.04597
[49,] 7.30333331 87.05873 66.13655 107.98092 73.63678 100.48069
[50,] 7.45444444 86.91296 66.81607 107.00986 74.02045 99.80548
[51,] 7.60555556 86.77798 67.68521 105.87076 74.52963 99.02634
[52,] 7.75666669 86.65379 68.72681 104.58078 75.15332 98.15427
[53,] 7.90777782 86.54040 69.91190 103.16890 75.87292 97.20787
[54,] 8.05888895 86.43779 71.19243 101.68316 76.65762 96.21797
[55,] 8.21000008 86.34598 72.48892 100.20304 77.45643 95.23553
[56,] 8.36111121 86.26496 73.67136 98.85857 78.18594 94.34398
[57,] 8.51222234 86.19473 74.53854 97.85092 78.71708 93.67239
[58,] 8.66333347 86.13530 74.85213 97.41846 78.89695 93.37365
[59,] 8.81444460 86.08665 74.73592 97.43739 78.80495 93.36835
[60,] 8.96555572 86.04880 74.59422 97.50339 78.70048 93.39712
[61,] 9.11666685 86.02174 74.63212 97.41137 78.71510 93.32839
[62,] 9.26777798 86.00548 74.85317 97.15778 78.85107 93.15988
[63,] 9.41888911 86.00000 75.05365 96.94635 78.97772 93.02228
[64,] 9.57000024 86.00000 74.77524 97.22476 78.79912 93.20088
[65,] 9.72111137 86.00000 73.92705 98.07295 78.25498 93.74502
[66,] 9.87222250 86.00000 72.70855 99.29145 77.47330 94.52670
[67,] 10.02333363 86.00000 71.28863 100.71137 76.56239 95.43761
[68,] 10.17444476 86.00000 69.78565 102.21435 75.59820 96.40180
[69,] 10.32555588 86.00000 68.27651 103.72349 74.63007 97.36993
[70,] 10.47666701 86.00000 66.81029 105.18971 73.68946 98.31054
[71,] 10.62777814 86.00000 65.41850 106.58150 72.79660 99.20340
[72,] 10.77888927 86.00000 64.12173 107.87827 71.96470 100.03530
[73,] 10.93000040 86.00000 62.93363 109.06637 71.20251 100.79749
[74,] 11.08111153 86.00000 61.86334 110.13666 70.51590 101.48410
[75,] 11.23222266 86.00000 60.91695 111.08305 69.90878 102.09122
[76,] 11.38333379 86.00000 60.09842 111.90158 69.38368 102.61632
[77,] 11.53444492 86.00000 59.41017 112.58983 68.94215 103.05785
[78,] 11.68555604 86.00000 58.85335 113.14665 68.58494 103.41506
[79,] 11.83666717 86.00000 58.42816 113.57184 68.31217 103.68783
[80,] 11.98777830 86.00000 58.13385 113.86615 68.12337 103.87663
[81,] 12.13888943 86.00000 57.96884 114.03116 68.01751 103.98249
[82,] 12.29000056 86.00000 57.93063 114.06937 67.99300 104.00700
[83,] 12.44111169 86.00000 58.01575 113.98425 68.04761 103.95239
[84,] 12.59222282 86.00000 58.21961 113.78039 68.17839 103.82161
[85,] 12.74333395 86.00000 58.53625 113.46375 68.38152 103.61848
[86,] 12.89444508 86.00000 58.95807 113.04193 68.65212 103.34788
[87,] 13.04555621 86.00000 59.47540 112.52460 68.98400 103.01600
[88,] 13.19666733 86.00000 60.07597 111.92403 69.36927 102.63073
[89,] 13.34777846 86.00000 60.74422 111.25578 69.79797 102.20203
[90,] 13.49888959 86.00000 61.46042 110.53958 70.25743 101.74257
[91,] 13.65000072 86.00000 62.19963 109.80037 70.73164 101.26836
[92,] 13.80111185 86.00000 62.93053 109.06947 71.20053 100.79947
[93,] 13.95222298 86.00000 63.61431 108.38569 71.63918 100.36082
[94,] 14.10333411 86.00000 64.20397 107.79603 72.01746 99.98254
[95,] 14.25444524 86.00000 64.64467 107.35533 72.30017 99.69983
[96,] 14.40555637 86.00000 64.87587 107.12413 72.44850 99.55150
[97,] 14.55666749 86.00000 64.83597 107.16403 72.42290 99.57710
[98,] 14.70777862 86.00000 64.46912 107.53088 72.18756 99.81244
[99,] 14.85888975 86.00000 63.73260 108.26740 71.71506 100.28494
[100,] 15.01000088 86.00000 62.60202 109.39798 70.98978 101.01022
> (pa50.1 <- predict(a50.1,interval = "both"))
z fit cb.lo cb.up ci.lo ci.up
[1,] 0.04999912 99.90710 74.07733 125.73687 81.65008 118.16412
[2,] 0.20111025 99.62822 76.47840 122.77803 83.26545 115.99098
[3,] 0.35222138 99.35216 78.09434 120.60999 84.32669 114.37763
[4,] 0.50333251 99.07894 79.18123 118.97665 85.01483 113.14305
[5,] 0.65444363 98.80855 79.97274 117.64436 85.49501 112.12209
[6,] 0.80555476 98.54099 80.61597 116.46601 85.87122 111.21076
[7,] 0.95666589 98.27626 81.12155 115.43098 86.15096 110.40157
[8,] 1.10777702 98.01437 81.32435 114.70439 86.21752 109.81122
[9,] 1.25888815 97.75530 81.10520 114.40541 85.98666 109.52394
[10,] 1.40999928 97.49907 80.74543 114.25270 85.65725 109.34088
[11,] 1.56111041 97.24566 80.47671 114.01461 85.39303 109.09830
[12,] 1.71222154 96.99509 80.40878 113.58140 85.27155 108.71864
[13,] 1.86333267 96.74735 80.53557 112.95913 85.28853 108.20617
[14,] 2.01444379 96.50244 80.72798 112.27691 85.35273 107.65216
[15,] 2.16555492 96.26037 80.71890 111.80184 85.27534 107.24540
[16,] 2.31666605 96.02112 80.42388 111.61836 84.99667 107.04557
[17,] 2.46777718 95.78471 80.13472 111.43469 84.72297 106.84644
[18,] 2.61888831 95.55112 80.04796 111.05429 84.59316 106.50908
[19,] 2.76999944 95.32037 80.22914 110.41160 84.65358 105.98716
[20,] 2.92111057 95.09245 80.61377 109.57113 84.85862 105.32628
[21,] 3.07222170 94.86736 80.97927 108.75545 85.05097 104.68374
[22,] 3.22333283 94.64510 81.09086 108.19934 85.06469 104.22552
[23,] 3.37444396 94.42567 81.07312 107.77823 84.98781 103.86353
[24,] 3.52555508 94.20908 80.91899 107.49916 84.81537 103.60278
[25,] 3.67666621 93.99531 80.38572 107.60491 84.37577 103.61485
[26,] 3.82777734 93.78438 79.43282 108.13594 83.64040 103.92836
[27,] 3.97888847 93.57628 78.32834 108.82422 82.79872 104.35384
[28,] 4.12999960 93.37101 77.26469 109.47733 81.98673 104.75529
[29,] 4.28111073 93.16857 76.35655 109.98059 81.28549 105.05165
[30,] 4.43222186 92.96896 75.66386 110.27407 80.73736 105.20057
[31,] 4.58333299 92.77219 75.20807 110.33630 80.35751 105.18686
[32,] 4.73444412 92.57824 74.98009 110.17639 80.13950 105.01698
[33,] 4.88555524 92.38713 74.94137 109.83289 80.05610 104.71815
[34,] 5.03666637 92.19884 75.01995 109.37774 80.05645 104.34124
[35,] 5.18777750 92.01339 75.10326 108.92352 80.06096 103.96582
[36,] 5.33888863 91.83077 75.01675 108.64480 79.94627 103.71527
[37,] 5.48999976 91.65099 74.62788 108.67409 79.61871 103.68326
[38,] 5.64111089 91.47403 74.02568 108.92238 79.14117 103.80688
[39,] 5.79222202 91.29990 73.30695 109.29285 78.58212 104.01769
[40,] 5.94333315 91.12861 72.55223 109.70499 77.99844 104.25878
[41,] 6.09444428 90.96014 71.82434 110.09594 77.43457 104.48572
[42,] 6.24555540 90.79451 71.17059 110.41844 76.92392 104.66511
[43,] 6.39666653 90.63171 70.62580 110.63762 76.49112 104.77230
[44,] 6.54777766 90.47174 70.21530 110.72818 76.15407 104.78941
[45,] 6.69888879 90.31460 69.95720 110.67201 75.92557 104.70364
[46,] 6.84999992 90.16030 69.86404 110.45655 75.81449 104.50611
[47,] 7.00111105 90.00882 69.94395 110.07369 75.82656 104.19109
[48,] 7.15222218 89.86018 70.20119 109.51917 75.96480 103.75556
[49,] 7.30333331 89.71437 70.63633 108.79240 76.22962 103.19911
[50,] 7.45444444 89.57138 71.24590 107.89687 76.61855 102.52421
[51,] 7.60555556 89.43123 72.02136 106.84111 77.12558 101.73689
[52,] 7.75666669 89.29392 72.94708 105.64075 77.73963 100.84820
[53,] 7.90777782 89.15943 73.99662 104.32223 78.44204 99.87681
[54,] 8.05888895 89.02777 75.12618 102.92936 79.20184 98.85370
[55,] 8.21000008 88.89895 76.26329 101.53460 79.96781 97.83009
[56,] 8.36111121 88.77295 77.28939 100.25652 80.65613 96.88977
[57,] 8.51222234 88.64979 78.02102 99.27856 81.13716 96.16243
[58,] 8.66333347 88.52946 78.24083 98.81809 81.25724 95.80168
[59,] 8.81444460 88.41196 78.06172 98.76220 81.09619 95.72773
[60,] 8.96555572 88.29729 77.85235 98.74223 80.91459 95.67999
[61,] 9.11666685 88.18546 77.79975 98.57116 80.84462 95.52629
[62,] 9.26777798 88.07645 77.90715 98.24575 80.88858 95.26432
[63,] 9.41888911 87.97028 77.98877 97.95178 80.91514 95.02541
[64,] 9.57000024 87.86693 77.63156 98.10230 80.63236 95.10150
[65,] 9.72111137 87.76642 76.75762 98.77522 79.98517 95.54767
[66,] 9.87222250 87.66874 75.54884 99.78864 79.10215 96.23533
[67,] 10.02333363 87.57389 74.15923 100.98855 78.09213 97.05565
[68,] 10.17444476 87.48187 72.69671 102.26704 77.03141 97.93233
[69,] 10.32555588 87.39269 71.23141 103.55396 75.96956 98.81581
[70,] 10.47666701 87.30633 69.80807 104.80459 74.93820 99.67446
[71,] 10.62777814 87.22281 68.45543 105.99018 73.95764 100.48797
[72,] 10.77888927 87.14212 67.19227 107.09196 73.04116 101.24308
[73,] 10.93000040 87.06425 66.03103 108.09747 72.19754 101.93097
[74,] 11.08111153 86.98922 64.98005 108.99840 71.43269 102.54576
[75,] 11.23222266 86.91703 64.04488 109.78917 70.75052 103.08353
[76,] 11.38333379 86.84766 63.22914 110.46618 70.15360 103.54172
[77,] 11.53444492 86.78112 62.53501 111.02724 69.64347 103.91878
[78,] 11.68555604 86.71742 61.96357 111.47127 69.22089 104.21395
[79,] 11.83666717 86.65654 61.51498 111.79811 68.88597 104.42712
[80,] 11.98777830 86.59850 61.18857 112.00844 68.63824 104.55877
[81,] 12.13888943 86.54329 60.98289 112.10369 68.47667 104.60991
[82,] 12.29000056 86.49091 60.89567 112.08616 68.39967 104.58216
[83,] 12.44111169 86.44136 60.92374 111.95899 68.40498 104.47775
[84,] 12.59222282 86.39465 61.06292 111.72638 68.48966 104.29964
[85,] 12.74333395 86.35076 61.30776 111.39376 68.64985 104.05167
[86,] 12.89444508 86.30971 61.65135 110.96807 68.88067 103.73875
[87,] 13.04555621 86.27149 62.08485 110.45812 69.17588 103.36710
[88,] 13.19666733 86.23609 62.59710 109.87509 69.52756 102.94463
[89,] 13.34777846 86.20353 63.17389 109.23318 69.92570 102.48136
[90,] 13.49888959 86.17381 63.79723 108.55038 70.35758 101.99003
[91,] 13.65000072 86.14691 64.44439 107.84943 70.80712 101.48670
[92,] 13.80111185 86.12284 65.08679 107.15889 71.25413 100.99155
[93,] 13.95222298 86.10161 65.68907 106.51415 71.67360 100.52961
[94,] 14.10333411 86.08320 66.20835 105.95806 72.03525 100.13116
[95,] 14.25444524 86.06763 66.59463 105.54063 72.30371 99.83155
[96,] 14.40555637 86.05489 66.79272 105.31707 72.43999 99.66979
[97,] 14.55666749 86.04498 66.74642 105.34354 72.40436 99.68560
[98,] 14.70777862 86.03790 66.40483 105.67098 72.16084 99.91496
[99,] 14.85888975 86.03366 65.72898 106.33834 71.68189 100.38542
[100,] 15.01000088 86.03224 64.69664 107.36784 70.95180 101.11268
> (pa25 <- predict(a25, interval = "both"))
z fit cb.lo cb.up ci.lo ci.up
[1,] 0.04999912 99.61886 71.59776 127.63995 81.64282 117.59489
[2,] 0.20111025 98.47926 73.36548 123.59304 82.36832 114.59020
[3,] 0.35222138 97.35819 74.29692 120.41947 82.56397 112.15241
[4,] 0.50333251 96.25565 74.66988 117.84143 82.40799 110.10331
[5,] 0.65444363 95.17164 74.73786 115.60542 82.06301 108.28027
[6,] 0.80555476 94.10616 74.66043 113.55188 81.63138 106.58094
[7,] 0.95666589 93.05920 74.44913 111.66927 81.12051 104.99789
[8,] 1.10777702 92.03077 73.92482 110.13673 80.41548 103.64606
[9,] 1.25888815 91.02087 72.95822 109.08352 79.43336 102.60838
[10,] 1.40999928 90.02950 71.85453 108.20447 78.36993 101.68906
[11,] 1.56111041 89.05665 70.86508 107.24823 77.38643 100.72687
[12,] 1.71222154 88.10234 70.10889 106.09578 76.55922 99.64545
[13,] 1.86333267 87.16655 69.57941 104.75369 75.88408 98.44901
[14,] 2.01444379 86.24929 69.13656 103.36201 75.27117 97.22740
[15,] 2.16555492 85.35055 68.49059 102.21051 74.53459 96.16652
[16,] 2.31666605 84.47034 67.54988 101.39081 73.61556 95.32513
[17,] 2.46777718 83.60867 66.63098 100.58635 72.71718 94.50015
[18,] 2.61888831 82.76552 65.94711 99.58393 71.97621 93.55482
[19,] 2.76999944 81.94089 65.56937 98.31242 71.43827 92.44352
[20,] 2.92111057 81.13480 65.42779 96.84181 71.05847 91.21113
[21,] 3.07222170 80.34723 65.28092 95.41354 70.68192 90.01254
[22,] 3.22333283 79.57819 64.87405 94.28233 70.14522 89.01116
[23,] 3.37444396 78.82768 64.34233 93.31302 69.53507 88.12029
[24,] 3.52555508 78.09570 63.67812 92.51327 68.84656 87.34483
[25,] 3.67666621 77.38224 62.61805 92.14643 67.91075 86.85373
[26,] 3.82777734 76.68731 61.11821 92.25642 66.69945 86.67517
[27,] 3.97888847 76.01091 59.46938 92.55245 65.39922 86.62260
[28,] 4.12999960 75.35304 57.88030 92.82577 64.14397 86.56211
[29,] 4.28111073 74.71369 56.47539 92.95199 63.01350 86.41389
[30,] 4.43222186 74.09288 55.31966 92.86610 62.04952 86.13623
[31,] 4.58333299 73.49059 54.43638 92.54479 61.26698 85.71420
[32,] 4.73444412 72.90683 53.81569 91.99796 60.65953 85.15413
[33,] 4.88555524 72.34159 53.41578 91.26740 60.20035 84.48283
[34,] 5.03666637 71.79489 53.15859 90.43119 59.83937 83.75040
[35,] 5.18777750 71.26671 52.92197 89.61145 59.49823 83.03518
[36,] 5.33888863 70.75706 52.51658 88.99753 59.05547 82.45865
[37,] 5.48999976 70.26593 51.79864 88.73323 58.41884 82.11303
[38,] 5.64111089 69.79334 50.86472 88.72196 57.65029 81.93639
[39,] 5.79222202 69.33927 49.81985 88.85869 56.81722 81.86133
[40,] 5.94333315 68.90373 48.75139 89.05608 55.97564 81.83182
[41,] 6.09444428 68.48672 47.72750 89.24595 55.16931 81.80413
[42,] 6.24555540 68.08824 46.79948 89.37700 54.43112 81.74536
[43,] 6.39666653 67.70828 46.00513 89.41144 53.78532 81.63124
[44,] 6.54777766 67.34685 45.37192 89.32179 53.24954 81.44417
[45,] 6.69888879 67.00395 44.91948 89.08843 52.83637 81.17153
[46,] 6.84999992 66.67958 44.66145 88.69771 52.55456 80.80460
[47,] 7.00111105 66.37374 44.60662 88.14085 52.40974 80.33773
[48,] 7.15222218 66.08642 44.75962 87.41322 52.40490 79.76794
[49,] 7.30333331 65.81763 45.12107 86.51419 52.54042 79.09484
[50,] 7.45444444 65.56737 45.68720 85.44753 52.81389 78.32085
[51,] 7.60555556 65.33564 46.44876 84.22251 53.21937 77.45190
[52,] 7.75666669 65.12243 47.38877 82.85609 53.74597 76.49889
[53,] 7.90777782 64.92775 48.47858 81.37693 54.37531 75.48019
[54,] 8.05888895 64.75160 49.67064 79.83256 55.07690 74.42631
[55,] 8.21000008 64.59398 50.88635 78.30161 55.80029 73.38767
[56,] 8.36111121 64.45488 51.99709 76.91268 56.46299 72.44678
[57,] 8.51222234 64.33432 52.80383 75.86481 56.93731 71.73133
[58,] 8.66333347 64.23228 53.07079 75.39377 57.07199 71.39257
[59,] 8.81444460 64.14877 52.92044 75.37710 56.94559 71.35194
[60,] 8.96555572 64.08378 52.75273 75.41484 56.81471 71.35286
[61,] 9.11666685 64.03733 52.77053 75.30413 56.80948 71.26518
[62,] 9.26777798 64.00940 52.97736 75.04144 56.93215 71.08665
[63,] 9.41888911 64.00000 53.17170 74.82830 57.05345 70.94655
[64,] 9.57000024 64.00000 52.89629 75.10371 56.87677 71.12323
[65,] 9.72111137 64.00000 52.05724 75.94276 56.33851 71.66149
[66,] 9.87222250 64.00000 50.85189 77.14811 55.56525 72.43475
[67,] 10.02333363 64.00000 49.44728 78.55272 54.66417 73.33583
[68,] 10.17444476 64.00000 47.96050 80.03950 53.71038 74.28962
[69,] 10.32555588 64.00000 46.46765 81.53235 52.75268 75.24732
[70,] 10.47666701 64.00000 45.01724 82.98276 51.82222 76.17778
[71,] 10.62777814 64.00000 43.64046 84.35954 50.93899 77.06101
[72,] 10.77888927 64.00000 42.35767 85.64233 50.11606 77.88394
[73,] 10.93000040 64.00000 41.18238 86.81762 49.36209 78.63791
[74,] 11.08111153 64.00000 40.12363 87.87637 48.68288 79.31712
[75,] 11.23222266 64.00000 39.18745 88.81255 48.08231 79.91769
[76,] 11.38333379 64.00000 38.37775 89.62225 47.56287 80.43713
[77,] 11.53444492 64.00000 37.69691 90.30309 47.12610 80.87390
[78,] 11.68555604 64.00000 37.14610 90.85390 46.77275 81.22725
[79,] 11.83666717 64.00000 36.72549 91.27451 46.50292 81.49708
[80,] 11.98777830 64.00000 36.43436 91.56564 46.31615 81.68385
[81,] 12.13888943 64.00000 36.27113 91.72887 46.21144 81.78856
[82,] 12.29000056 64.00000 36.23333 91.76667 46.18719 81.81281
[83,] 12.44111169 64.00000 36.31754 91.68246 46.24121 81.75879
[84,] 12.59222282 64.00000 36.51920 91.48080 46.37058 81.62942
[85,] 12.74333395 64.00000 36.83242 91.16758 46.57152 81.42848
[86,] 12.89444508 64.00000 37.24969 90.75031 46.83920 81.16080
[87,] 13.04555621 64.00000 37.76144 90.23856 47.16750 80.83250
[88,] 13.19666733 64.00000 38.35554 89.64446 47.54862 80.45138
[89,] 13.34777846 64.00000 39.01658 88.98342 47.97269 80.02731
[90,] 13.49888959 64.00000 39.72506 88.27494 48.42719 79.57281
[91,] 13.65000072 64.00000 40.45630 87.54370 48.89630 79.10370
[92,] 13.80111185 64.00000 41.17931 86.82069 49.36012 78.63988
[93,] 13.95222298 64.00000 41.85572 86.14428 49.79405 78.20595
[94,] 14.10333411 64.00000 42.43902 85.56098 50.16825 77.83175
[95,] 14.25444524 64.00000 42.87497 85.12503 50.44792 77.55208
[96,] 14.40555637 64.00000 43.10368 84.89632 50.59464 77.40536
[97,] 14.55666749 64.00000 43.06421 84.93579 50.56932 77.43068
[98,] 14.70777862 64.00000 42.70131 85.29869 50.33651 77.66349
[99,] 14.85888975 64.00000 41.97273 86.02727 49.86911 78.13089
[100,] 15.01000088 64.00000 40.85435 87.14565 49.15165 78.84835
> (pa75 <- predict(a75, interval = "both"))
z fit cb.lo cb.up ci.lo ci.up
[1,] 0.04999912 100 100 100 100 100
[2,] 0.20111025 100 100 100 100 100
[3,] 0.35222138 100 100 100 100 100
[4,] 0.50333251 100 100 100 100 100
[5,] 0.65444363 100 100 100 100 100
[6,] 0.80555476 100 100 100 100 100
[7,] 0.95666589 100 100 100 100 100
[8,] 1.10777702 100 100 100 100 100
[9,] 1.25888815 100 100 100 100 100
[10,] 1.40999928 100 100 100 100 100
[11,] 1.56111041 100 100 100 100 100
[12,] 1.71222154 100 100 100 100 100
[13,] 1.86333267 100 100 100 100 100
[14,] 2.01444379 100 100 100 100 100
[15,] 2.16555492 100 100 100 100 100
[16,] 2.31666605 100 100 100 100 100
[17,] 2.46777718 100 100 100 100 100
[18,] 2.61888831 100 100 100 100 100
[19,] 2.76999944 100 100 100 100 100
[20,] 2.92111057 100 100 100 100 100
[21,] 3.07222170 100 100 100 100 100
[22,] 3.22333283 100 100 100 100 100
[23,] 3.37444396 100 100 100 100 100
[24,] 3.52555508 100 100 100 100 100
[25,] 3.67666621 100 100 100 100 100
[26,] 3.82777734 100 100 100 100 100
[27,] 3.97888847 100 100 100 100 100
[28,] 4.12999960 100 100 100 100 100
[29,] 4.28111073 100 100 100 100 100
[30,] 4.43222186 100 100 100 100 100
[31,] 4.58333299 100 100 100 100 100
[32,] 4.73444412 100 100 100 100 100
[33,] 4.88555524 100 100 100 100 100
[34,] 5.03666637 100 100 100 100 100
[35,] 5.18777750 100 100 100 100 100
[36,] 5.33888863 100 100 100 100 100
[37,] 5.48999976 100 100 100 100 100
[38,] 5.64111089 100 100 100 100 100
[39,] 5.79222202 100 100 100 100 100
[40,] 5.94333315 100 100 100 100 100
[41,] 6.09444428 100 100 100 100 100
[42,] 6.24555540 100 100 100 100 100
[43,] 6.39666653 100 100 100 100 100
[44,] 6.54777766 100 100 100 100 100
[45,] 6.69888879 100 100 100 100 100
[46,] 6.84999992 100 100 100 100 100
[47,] 7.00111105 100 100 100 100 100
[48,] 7.15222218 100 100 100 100 100
[49,] 7.30333331 100 100 100 100 100
[50,] 7.45444444 100 100 100 100 100
[51,] 7.60555556 100 100 100 100 100
[52,] 7.75666669 100 100 100 100 100
[53,] 7.90777782 100 100 100 100 100
[54,] 8.05888895 100 100 100 100 100
[55,] 8.21000008 100 100 100 100 100
[56,] 8.36111121 100 100 100 100 100
[57,] 8.51222234 100 100 100 100 100
[58,] 8.66333347 100 100 100 100 100
[59,] 8.81444460 100 100 100 100 100
[60,] 8.96555572 100 100 100 100 100
[61,] 9.11666685 100 100 100 100 100
[62,] 9.26777798 100 100 100 100 100
[63,] 9.41888911 100 100 100 100 100
[64,] 9.57000024 100 100 100 100 100
[65,] 9.72111137 100 100 100 100 100
[66,] 9.87222250 100 100 100 100 100
[67,] 10.02333363 100 100 100 100 100
[68,] 10.17444476 100 100 100 100 100
[69,] 10.32555588 100 100 100 100 100
[70,] 10.47666701 100 100 100 100 100
[71,] 10.62777814 100 100 100 100 100
[72,] 10.77888927 100 100 100 100 100
[73,] 10.93000040 100 100 100 100 100
[74,] 11.08111153 100 100 100 100 100
[75,] 11.23222266 100 100 100 100 100
[76,] 11.38333379 100 100 100 100 100
[77,] 11.53444492 100 100 100 100 100
[78,] 11.68555604 100 100 100 100 100
[79,] 11.83666717 100 100 100 100 100
[80,] 11.98777830 100 100 100 100 100
[81,] 12.13888943 100 100 100 100 100
[82,] 12.29000056 100 100 100 100 100
[83,] 12.44111169 100 100 100 100 100
[84,] 12.59222282 100 100 100 100 100
[85,] 12.74333395 100 100 100 100 100
[86,] 12.89444508 100 100 100 100 100
[87,] 13.04555621 100 100 100 100 100
[88,] 13.19666733 100 100 100 100 100
[89,] 13.34777846 100 100 100 100 100
[90,] 13.49888959 100 100 100 100 100
[91,] 13.65000072 100 100 100 100 100
[92,] 13.80111185 100 100 100 100 100
[93,] 13.95222298 100 100 100 100 100
[94,] 14.10333411 100 100 100 100 100
[95,] 14.25444524 100 100 100 100 100
[96,] 14.40555637 100 100 100 100 100
[97,] 14.55666749 100 100 100 100 100
[98,] 14.70777862 100 100 100 100 100
[99,] 14.85888975 100 100 100 100 100
[100,] 15.01000088 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)
>
