centiles.pred: Creating predictive centiles values

Description Usage Arguments Value Warning Note Author(s) References See Also Examples

View source: R/centile-pred.R

Description

This function creates predictive centiles curves for new x-values given a GAMLSS fitted model. The function has three options: i) for given new x-values and given percentage centiles calculates a matrix containing the centiles values for y, ii) for given new x-values and standard normalized centile values calculates a matrix containing the centiles values for y, iii) for given new x-values and new y-values calculates the z-scores. A restriction of the function is that it applies to models with only one explanatory variable.

Usage

1
2
3
4
5
6
centiles.pred(obj, type = c("centiles", "z-scores", "standard-centiles"), 
             xname = NULL, xvalues = NULL, power = NULL, yval = NULL, 
             cent = c(0.4, 2, 10, 25, 50, 75, 90, 98, 99.6), 
             dev = c(-4, -3, -2, -1, 0, 1, 2, 3, 4), calibration = FALSE,
             plot = FALSE, legend = TRUE, 
             ...)

Arguments

obj

a fitted gamlss object from fitting a gamlss continuous distribution

type

the default, "centiles", gets the centiles values given in the option cent. type="standard-centiles" gets the standard centiles given in the dev. type="z-scores" gets the z-scores for given y and x new values

xname

the name of the unique explanatory variable (it has to be the same as in the original fitted model)

xvalues

the new values for the explanatory variable where the prediction will take place

power

if power transformation is needed (but read the note below)

yval

the response values for a given x required for the calculation of "z-scores"

cent

a vector with elements the % centile values for which the centile curves have to be evaluated

dev

a vector with elements the standard normalized values for which the centile curves have to be evaluated in the option type="standard-centiles"

calibration

whether to calibrate the "centiles", the default is calibrate=FALSE

plot

whether to plot the "centiles" or the "standard-centiles", the default is plot=FALSE

legend

whether a legend is required in the plot or not, the default is legent=TRUE

...

for extra arguments

Value

a vector (for option type="z-scores") or a matrix for options type="centiles" or type="standard-centiles" containing the appropriate values

Warning

See example below of how to use the function when power transformation is used for the x-variables

Note

The power option should be only used if the model

Author(s)

Mikis Stasinopoulos , d.stasinopoulos@londonmet.ac.uk, based on ideas of Elaine Borghie from the World Health Organization

References

Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, pp 507-554.

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC. An older version can be found in https://www.gamlss.com/.

Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, https://www.jstatsoft.org/v23/i07/.

Stasinopoulos D. M., Rigby R.A., Heller G., Voudouris V., and De Bastiani F., (2017) Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC.

(see also https://www.gamlss.com/).

See Also

gamlss, centiles, centiles.split

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
## bring the data and fit the model
data(abdom)
a<-gamlss(y~pb(x),sigma.fo=~pb(x), data=abdom, family=BCT)
## plot the centiles
centiles(a,xvar=abdom$x)
##-----------------------------------------------------------------------------
## the first use of the function centiles.pred()
## to calculate the centiles at new x values
##-----------------------------------------------------------------------------
newx<-seq(12,40,2)
mat <- centiles.pred(a, xname="x", xvalues=newx )
mat
## now plot the centile curves  
 mat <- centiles.pred(a, xname="x",xvalues=newx, plot=TRUE )
##-----------------------------------------------------------------------------
## the second use of the function centiles.pred()
## to calculate (nornalised) standard-centiles for new x
## values using the fitted model
##-----------------------------------------------------------------------------
newx <- seq(12,40,2)
mat <- centiles.pred(a, xname="x",xvalues=newx, type="standard-centiles" )
mat
## now plot the standard centiles  
mat <- centiles.pred(a, xname="x",xvalues=newx, type="standard-centiles",
       plot = TRUE )
##-----------------------------------------------------------------------------
## the third use of the function centiles.pred()
##  if we have new x and y values what are their z-scores?
##-----------------------------------------------------------------------------
# create new y and x values and plot them in the previous plot
newx <- c(20,21.2,23,20.9,24.2,24.1,25)
newy <- c(130,121,123,125,140,145,150)
for(i in 1:7) points(newx[i],newy[i],col="blue")
## now calculate their z-scores
znewx <- centiles.pred(a, xname="x",xvalues=newx,yval=newy, type="z-scores" )
znewx
## Not run: 
##-----------------------------------------------------------------------------
## What we do if the x variables is transformed?
##----------------------------------------------------------------------------
##  case 1 : transformed x-variable within the formula
##----------------------------------------------------------------------------
## fit model
aa <- gamlss(y~pb(x^0.5),sigma.fo=~pb(x^0.5), data=abdom, family=BCT)
## centiles is working in this case
centiles(aa, xvar=abdom$x, legend = FALSE)
## get predict for values of x at 12, 14, ..., 40
mat <- centiles.pred(aa, xname="x", xvalues=seq(12,40,2), plot=TRUE )
mat
# plot all prediction points
xx <- rep(mat[,1],9)
yy <- unlist(mat[,2:10])
points(xx,yy,col="red")
##----------------------------------------------------------------------------
##  case 2 : the x-variable is previously transformed 
##----------------------------------------------------------------------------
nx <- abdom$x^0.5
aa <- gamlss(y~pb(nx),sigma.fo=~pb(nx), data=abdom, family=BCT)
centiles(aa, xvar=abdom$x)
# equivalent to fitting
newd<-data.frame( abdom, nx=abdom$x^0.5)
aa1 <- gamlss(y~pb(nx),sigma.fo=~pb(nx), family=BCT, data=newd)
centiles(aa1, xvar=abdom$x)
# getting the centiles at x equal to 12, 14, ...40
mat <-  centiles.pred(aa, xname="nx", xvalues=seq(12,40,2), power=0.5, 
         data=newd, plot=TRUE)
# plot all prediction points         
xxx <- rep(mat[,1],9)
yyy <- unlist(mat[,2:10])
points(xxx,yyy,col="red")
# the idea is that if the transformed x-variable is used in the fit
# the power argument has to used in centiles.pred()

## End(Not run)

Example output

Loading required package: splines
Loading required package: gamlss.data

Attaching package: 'gamlss.data'

The following object is masked from 'package:datasets':

    sleep

Loading required package: gamlss.dist
Loading required package: MASS
Loading required package: nlme
Loading required package: parallel
 **********   GAMLSS Version 5.1-3  ********** 
For more on GAMLSS look at http://www.gamlss.org/
Type gamlssNews() to see new features/changes/bug fixes.

GAMLSS-RS iteration 1: Global Deviance = 4771.925 
GAMLSS-RS iteration 2: Global Deviance = 4771.039 
GAMLSS-RS iteration 3: Global Deviance = 4770.999 
GAMLSS-RS iteration 4: Global Deviance = 4770.994 
GAMLSS-RS iteration 5: Global Deviance = 4770.993 
% of cases below  0.4 centile is  0.3278689 
% of cases below  2 centile is  2.459016 
% of cases below  10 centile is  8.688525 
% of cases below  25 centile is  26.22951 
% of cases below  50 centile is  50.16393 
% of cases below  75 centile is  73.77049 
% of cases below  90 centile is  90 
% of cases below  98 centile is  98.03279 
% of cases below  99.6 centile is  99.67213 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
    x      C0.4        C2       C10       C25       C50       C75       C90
1  12  43.85374  47.01984  50.74981  53.54954  56.68358  60.02206  63.39547
2  14  64.99847  69.26230  74.24949  77.96894  82.10964  86.49528  90.90189
3  16  86.34623  91.48407  97.45352 101.87910 106.78079 111.94504 117.10700
4  18 107.36397 113.15337 119.83864 124.76784 130.20180 135.89910 141.56677
5  20 127.74961 134.02414 141.23116 146.51998 152.32682 158.38963 164.39613
6  22 147.47956 154.16417 161.80999 167.39985 173.51762 179.88399 186.17077
7  24 167.03055 174.15663 182.28328 188.20903 194.67984 201.39796 208.01698
8  26 186.43408 194.11651 202.86341 209.23224 216.17832 223.38070 230.46799
9  28 205.25984 213.65320 223.20618 230.15974 237.74154 245.60093 253.33262
10 30 223.05929 232.30271 242.82961 250.49626 258.85943 267.53296 276.06958
11 32 239.94648 250.14305 261.76914 270.24518 279.49955 289.10620 298.56977
12 34 255.83969 267.02269 279.79078 289.11063 299.29674 309.88185 320.32017
13 36 270.63078 282.78191 296.67371 306.82572 317.93249 329.48630 340.89150
14 38 284.50070 297.58358 312.55868 323.51411 335.51087 348.00232 360.34458
15 40 297.87102 311.82018 327.80182 339.50337 352.32631 365.68783 378.89934
         C98     C99.6
1   68.59603  73.80262
2   97.64928 104.35232
3  124.96154 132.70885
4  150.14145 158.54397
5  173.43863 182.24991
6  195.59846 204.74447
7  217.91579 227.48904
8  241.05126 251.26909
9  264.87444 276.01361
10 288.82014 301.13376
11 312.72020 326.40256
12 335.94760 351.07956
13 357.98730 374.56407
14 378.86564 396.84720
15 398.74199 418.02565
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
    x        -4        -3        -2        -1         0         1         2
1  12  35.15895  41.84842  47.29014  52.05426  56.68358  61.77326  68.18800
2  14  53.12637  62.28239  69.62499  75.98492  82.10964  88.78588  97.12182
3  16  71.85517  83.05595  91.91962  99.52115 106.78079 114.63156 124.34960
4  18  90.83984 103.63825 113.64261 122.14434 130.20180 138.85213 149.47547
5  20 109.65487 123.69455 134.55293 143.70767 152.32682 161.52221 172.73819
6  22 128.04345 143.14495 154.72632 164.42963 173.51762 183.16524 194.86972
7  24 146.19052 162.39874 174.75499 185.06194 194.67984 204.85446 217.15176
8  26 163.89552 181.43418 194.76105 205.85078 216.17832 227.08281 240.23506
9  28 180.61855 199.79570 214.35725 226.46805 237.74154 249.63991 263.98447
10 30 195.95489 217.04471 233.07831 246.42555 258.85943 271.99194 287.83667
11 32 210.11595 233.31792 250.99914 265.74380 279.49955 294.04831 311.62811
12 34 223.20982 248.57772 267.96226 284.15998 299.29674 315.33167 334.74071
13 36 235.26757 262.74845 283.80351 301.43179 317.93249 335.43951 356.66614
14 38 246.51544 276.02211 298.68419 317.69210 335.51087 354.44325 377.43348
15 40 257.44461 288.83780 312.99423 333.28385 352.32631 372.58122 397.20691
           3         4
1   77.56787  93.78078
2  109.16934 129.64914
3  138.24413 161.50485
4  164.51569 189.34583
5  188.48368 214.16921
6  211.19192 237.56934
7  234.22074 261.62408
8  258.44418 287.57313
9  283.83334 315.56069
10 309.78243 344.90894
11 336.02216 375.16925
12 361.73053 405.17351
13 386.24510 433.99563
14 409.53114 461.48714
15 431.63893 487.49014
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
[1] -2.445938 -3.983711 -4.692154 -3.552482 -4.392778 -4.108256 -4.250223
GAMLSS-RS iteration 1: Global Deviance = 4773.346 
GAMLSS-RS iteration 2: Global Deviance = 4771.828 
GAMLSS-RS iteration 3: Global Deviance = 4771.785 
GAMLSS-RS iteration 4: Global Deviance = 4771.779 
GAMLSS-RS iteration 5: Global Deviance = 4771.778 
% of cases below  0.4 centile is  0.3278689 
% of cases below  2 centile is  2.459016 
% of cases below  10 centile is  8.52459 
% of cases below  25 centile is  26.39344 
% of cases below  50 centile is  50.16393 
% of cases below  75 centile is  74.42623 
% of cases below  90 centile is  90.16393 
% of cases below  98 centile is  98.03279 
% of cases below  99.6 centile is  99.67213 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
    x      C0.4        C2       C10       C25       C50       C75       C90
1  12  42.69481  45.82570  49.50971  52.27247  55.36394  58.65747  61.98763
2  14  65.03564  69.32899  74.34107  78.07346  82.22485  86.62004  91.03685
3  16  86.58085  91.73915  97.71904 102.14454 107.04073 112.19600 117.34864
4  18 107.48667 113.27310 119.93927 124.84526 130.24726 135.90717 141.53688
5  20 127.80932 134.07786 141.26163 146.52397 152.29533 158.31743 164.28324
6  22 147.60275 154.30719 161.95995 167.54603 173.65399 180.00748 186.28235
7  24 167.16346 174.34787 182.52642 188.48204 194.98080 201.72660 208.37518
8  26 186.39931 194.17256 203.00897 209.43561 216.44087 223.70443 230.85560
9  28 204.99689 213.49241 223.14711 230.16709 237.81746 245.74814 253.55438
10 30 222.68119 232.01463 242.62711 250.34706 258.76357 267.49201 276.08692
11 32 239.44648 249.69704 261.36383 269.85820 279.12592 288.74453 298.22317
12 34 255.29150 266.49058 279.25185 288.55270 298.70931 309.26007 319.66661
13 36 270.24447 282.39181 296.25011 306.36122 317.41260 328.90355 340.24778
14 38 284.42645 297.50897 312.45140 323.36466 335.30319 347.72777 360.00456
15 40 298.02957 312.02057 328.01732 339.71145 352.51432 365.84931 379.03625
         C98     C99.6
1   67.12847  72.28531
2   97.80482 104.53665
3  125.19280 132.93726
4  150.05758 158.41442
5  173.26897 182.03337
6  195.69898 204.84538
7  218.32797 227.96797
8  241.54692 251.88704
9  265.22200 276.50297
10 288.93947 301.37285
11 312.41007 326.14844
12 335.25901 350.37687
13 357.26379 373.78249
14 378.43890 396.35608
15 398.85627 418.14132
GAMLSS-RS iteration 1: Global Deviance = 4773.346 
GAMLSS-RS iteration 2: Global Deviance = 4771.828 
GAMLSS-RS iteration 3: Global Deviance = 4771.785 
GAMLSS-RS iteration 4: Global Deviance = 4771.779 
GAMLSS-RS iteration 5: Global Deviance = 4771.778 
% of cases below  0.4 centile is  0.3278689 
% of cases below  2 centile is  2.459016 
% of cases below  10 centile is  8.52459 
% of cases below  25 centile is  26.39344 
% of cases below  50 centile is  50.16393 
% of cases below  75 centile is  74.42623 
% of cases below  90 centile is  90.16393 
% of cases below  98 centile is  98.03279 
% of cases below  99.6 centile is  99.67213 
GAMLSS-RS iteration 1: Global Deviance = 4773.346 
GAMLSS-RS iteration 2: Global Deviance = 4771.828 
GAMLSS-RS iteration 3: Global Deviance = 4771.785 
GAMLSS-RS iteration 4: Global Deviance = 4771.779 
GAMLSS-RS iteration 5: Global Deviance = 4771.778 
% of cases below  0.4 centile is  0.3278689 
% of cases below  2 centile is  2.459016 
% of cases below  10 centile is  8.52459 
% of cases below  25 centile is  26.39344 
% of cases below  50 centile is  50.16393 
% of cases below  75 centile is  74.42623 
% of cases below  90 centile is  90.16393 
% of cases below  98 centile is  98.03279 
% of cases below  99.6 centile is  99.67213 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
new prediction 
New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
Warning messages:
1: In plot.window(...) : "data" is not a graphical parameter
2: In plot.xy(xy, type, ...) : "data" is not a graphical parameter
3: In axis(side = side, at = at, labels = labels, ...) :
  "data" is not a graphical parameter
4: In axis(side = side, at = at, labels = labels, ...) :
  "data" is not a graphical parameter
5: In box(...) : "data" is not a graphical parameter
6: In title(...) : "data" is not a graphical parameter

gamlss documentation built on March 31, 2021, 5:10 p.m.