Demonstration for creating continuous norms with cNORM

Wolfgang Lenhard, Alexandra Lenhard


cNORM is a package for the R environment for statistical computing that aims at generating continuous test norms in psychometrics and biometrics and to analyze the model fit. It is based on the approach of A. Lenhard et al. (2016, 2019).

The method stems from psychometric test construction and was developed to create continuous norms for age or grade in performance assessment (e. g. vocabulary development, A. Lenhard, Lenhard, Segerer & Suggate, 2015; reading and writing development, W. Lenhard, Lenhard & Schneider, 2017). It can however be applied wherever test data like psychological (e. g. intelligence), physiological (e. g. weight) or other measures are dependent on continuous (e.g., age) or discrete (e.g., sex or test mode) explanatory variables. It has been applied to biometric data (fetal growth in dependence of gestation week) and macro economic data as well.

The package estimates percentile curves in dependence of the explanatory variable (e. g. schooling duration, age …) via Taylor polynomials, thus offering several advantages:

The rationale of the approach is to rank the results in the different age cohorts (or use a sliding window in case the data is distributed over a large age interval) and thus to determine the observed norm scores (= location). Afterwards, powers of the age specific location and of the age are computed, as well as all linear interactions. Finally, the data is fitted by a hyperplane via multiple regression and the most relevant terms are identified:

In a nutshell: Establishing continuous norming models

The ‘cnorm’ method combines most of the steps in one go. The example in a nutshell already suffices for establishing norm scores. It conducts the ranking, the computation of powers and the modeling. A detailed explanation of the distinct steps follows afterwards.


# We will use the internal dataset elfe on reading comprehension
# that includes raw scores and grade as a grouping variable.

# Most easy example: conventional norming without age or group.
# Just get the manifest norm scores and as well modelled via polynomial regression.
model <- cnorm(raw = elfe$raw)

# The method 'cnorm' prepares the data and establishes a model for continuous
# norming. Please specify a raw score and a grouping vector.
model <- cnorm(raw = elfe$raw, group = elfe$group)

# Fine tune model selection by checking model fit and by
# visual inspection, e. g. with the following checks
plot(model, "series", start=4, end = 8)   # series of percentile plots
plot(model, "subset")                     # plot information function$data)                      # run repeated cross validation

# Select final model
model <- cnorm(raw = elfe$raw, group = elfe$group, terms = 3)

# generate norm score table with 90% confiden interval and
# a reliability of .94 for several grade levels
normTable(c(3, 3.2, 3.4, 3.6), model, CI = .9, reliability = .94)

# To set powers of age and location independently, please specify k as the
# power parameter for location and t for age, e. g. for a cubic age trend
# with the ppvt demo data:
model.ppvt <- cnorm(raw = ppvt$raw, group = ppvt$group, k = 4, t = 3)

In the following, the single steps in detail:

1. Data Preparation and modelling

If a sufficiently large and representative sample has been established (missings will be excluded casewise), then the data must first be imported. It is advisable to start with a simply structured data object of type data.frame or numeric vectors containing raw score and grouping or the continuous explanatory variable. This explanatory variable in psychometric performance tests is usually age. We therefore refer to this variable as ‘age’ in the following. In fact, however, the explanatory variable is not necessarily age. A training or schooling duration or other explanatory variables can also be included in the modeling. However, it must be an interval-scaled (or, as the case may be, dichotomous) variable. Finally, a grouping variable is required to divide the explanatory variable into smaller standardization groups (e.g. grades or age groups). The method is relatively robust against changes in the granularity of the group subdivision. For example, the result of the standardization only marginally depends on whether one chooses half-year or full-year gradations (see A. Lenhard, Lenhard, Suggate & Segerer, 2016). The more the variable to be measured co-varies with the explanatory variable (e. g. a fast development over age in an intelligence test), the more groups should be formed beforehand to capture the trajectories adequately. By standard, we assign the variable name “group” to the grouping variable.

Conventional norming

In case, you are not interested in continuously modelling norm scores over age, cNORM can be used for conventional norming as well. It returns the manifest norm scores and percentiles but as well establishes a polynomial regression model that closes all the missings in the norm table, smoothes error variance and can be used for cautiously extrapolating more extreme raw scores.

# Apply conventional norming. You simply have to provide a raw score vector.
# Additionally, you can specify ranking order (parameter 'descend'), the 
# degree of the polynomial (parameter 'k') or the number of terms in the model
# via the parameter 'terms'.

model <- cnorm(raw=elfe$raw, terms=4)

Recoding Continuous Variables

No, let’s start with a detailed explanation of continuous norming. If, when using cNORM, you initially only have the continuous age variable available, it is advisable to as well recode it into a discrete grouping variable. Please take care, that the grouping variable aligns with the continuous variable. The values in the grouping variable should represent the group means. The following code could be helpful (another possibility is the’rankBySlidingWindow’ function described below):

# Creates a grouping variable for a fictitious age variable
# for children age 2 to 18 for the ppvt demo dataset. That way, 
# the age variable is recoded into a discrete group variable 
# with 12 distinct groups. The resulting vector specifies the 
# mean age of each group.

data$group <- getGroups(ppvt$age, 12)

Of course, it is also possible to use a data set for which standard scores already exist for individual age groups. Please pay attention, that the grouping variable corresponds to the group mean in case, you use a continuous age variable later on.

Sample Data

For demonstration purposes, cNORM includes a cleaned data set from a German test standardization (ELFE 1-6, W. Lenhard & Schneider, 2006, subtest sentence comprehension) that will be used for demonstrating the method. Another large (but unrepresentative) data set for demonstration purposes stems from the adaption of a vocabulary test to the German language (PPVT-4, A. Lenhard, Lenhard, Segerer & Suggate, 2015). For biometric modeling, it includes a large CDC dataset (N > 45,000) for growth curves from age 2 to 25 (weight, height, BMI; CDC, 2012) and for macro economical and sociological data the data on mortality and life expectancy at birth from 1960 to 2017 from the World Bank. You can retrieve information on the data by typing ?elfe, ?ppvt, ?CDC, ?life or ?mortality on the R console.

# Displays the first lines of the data of the example dataset 'elfe'
#>   personID group raw
#> 1      197     2   0
#> 2      295     2   0
#> 3     1261     2   0
#> 4     1317     2   0
#> 5     1331     2   0
#> 6      239     2   1

Data Exploration

As you can see, there is no age variable in the data set ‘elfe’, only a person ID, a raw score and a grouping variable. In this case, the grouping variable also serves as a continuous explanatory variable, since children were only examined at the very beginning and in the exact middle of the school year during the test standardization. For example, the value 2.0 means that the children were at the beginning of the second school year, the value 2.5 means that the children were examined in the middle of the second school year. In the ‘elfe’ data set there are seven groups with 200 cases each, i.e. a total of 1400 cases.

# Display some descriptive results by group
by(elfe$raw, elfe$group, summary)
#> elfe$group: 2
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00    4.00    7.00    7.32   10.00   23.00 
#> ------------------------------------------------------------ 
#> elfe$group: 2.5
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00    7.00   11.00   10.88   15.00   25.00 
#> ------------------------------------------------------------ 
#> elfe$group: 3
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00   11.00   14.00   14.51   19.00   28.00 
#> ------------------------------------------------------------ 
#> elfe$group: 3.5
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>     4.0    12.0    16.0    15.9    19.0    28.0 
#> ------------------------------------------------------------ 
#> elfe$group: 4
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    6.00   16.75   20.00   19.72   23.00   28.00 
#> ------------------------------------------------------------ 
#> elfe$group: 4.5
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>     1.0    19.0    22.0    21.3    25.0    28.0 
#> ------------------------------------------------------------ 
#> elfe$group: 5
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    2.00   19.00   23.00   22.27   26.00   28.00


The next step is to rank each person in each group using the rankByGroup function. It is already done by the cnorm function, so rankByGroup and subsequently computePowers is not necessary, if you use the comprehensive cnorm() function.

The function rankByGroups and as well rankBySlidingWindow return percentiles and also perform a normal-rank transformation in which T-Scores (M = 50, SD = 10) are returned by default. In principle, our mathematical method also works without normal rank transformation, i.e., the method could theoretically also be carried out with the percentiles. This is useful, for example, if you want to enter a variable that deviates extremely from the normal distribution or follows a completely different distribution. For most psychological or physical scales, however, the distributions are still sufficiently similar to the normal distribution even with strong bottom and ceiling effects. In these cases, the normal-rank transformation usually increases the model fit and facilitates the further processing of the data. In addition to T-Scores, the standard scores can also be expressed as z- or IQ-Scores or specified as desired. For bindings, RankIt (default), Blom, van der Warden, Tukey, Levenbach and Filliben, Yu & Huang is available.

# Determine percentiles by group
normData <- rankByGroup(elfe, group = elfe$group) 

To change the ranking method, please specify a method index with method = x (x = method index; see ?rankByGroup). The standard score can be specified as T-Score, IQ-Score, z-Score or by means of a double vector of M and SD, e.g. scale = c(10, 3) for Wechsler subtest scaled scores. The grouping variable can be deactivated by setting group = FALSE. The normal-rank transformation is then applied to the entire sample.

Please note that there is a second function for determining the rank, which works without discrete grouping variables. The rank of each individual subject is then estimated based on the continuous explanatory variables using a sliding window. The width of this window can be specified individually. In the case of a continuous age variable, the specification width = 0.5 means, for example, that the width of the window is half a year. As a consequence, the rank of a test persons is based on all participants who are no more than 3 months younger or older than the test person in question, i. e., the group comprises a total of 6 months.

# Percentile generation by a sliding window of the size 'width'
normData2 <- rankBySlidingWindow(data = elfe, age = elfe$group, raw = elfe$raw, width = 0.5)

Please note that the ‘rankBySlidingWindow’ function only makes sense if the age variable is actually continuous. In the ‘elfe’ data set the variable ‘group’ serves as continuous explanatory variable as well as discrete grouping variable. Therefore, with the function ‘rankBySlidingWindow’ we obtain the same standard scores as with the function ‘rankByGroup’ in this specific case.

In order to compensate for imbalances in the representativeness of the dataset, weights can be applied via the ‘weights’ parameter. In that case, weighted percentiles are computed and the weights are as well used in the regression modeling. Please have a look at the ‘WeightedRegression’ vignette on how to establish the weights.

Both ranking functions (‘rankBySlidingWindow’ and ‘rankByGroup’) add two additional columns, namely ‘percentile’ and ‘normValue’. In addition, descriptive information about each group is added, namely n, m, md and sd. Descriptive results are only necessary under certain circumstances. The creation of these variables can be deactivated via the parameter ‘descriptives’.

Computing powers and interactions

At this point, where many test developers already stop standardization, the actual modeling process begins. A function is determined which expresses the raw score as a function of the latent person parameter l and the explanatory variable. In the following, we will refer to the latter variable as ‘a’. In the ‘elfe’ example, we use the discrete variable ‘group’ for a. If there is an additional continuous age variable, it should be used instead as ‘a’ because of its higher precision.

To retrieve the mathematical model, all powers of the variables ‘l’ and ‘a’ up to a certain exponent k must be computed. Subsequently, all interactions between these powers must also be calculated by simple multiplication. As a rule of thumb, k > 5 leads to over-adjustment. In general, k = 4 or even k = 3 will already be sufficient to model human performance data with adequate precision. Please use the following function for the calculation:

# Calculation of powers and interactions up to k = 4
normData <- computePowers(normData, k = 4, norm = "normValue", age = "group") 

The data set now has 24 new variables ( \(2*k + k^{2}\) ), namely L1, L2, L3, L4 (powers of the norm value), A1, A2, A3, A4 (powers of the grouping variable) and the linear combinations L1A1, L2A1L4A3, L4A4.

It is as well possible to set the power parameter for location independently from age. If you want to model the age specific distributions with the power of 5, but a cubic age trajectory is sufficient, please additionally specify ‘t’ for age:

# Calculation of powers and interactions up to the fifth power of location and a
# cubic age trajectory (t = 3)
normData <- computePowers(normData, k = 5, t = 3, norm = "normValue", age = "group") 


We now want to find a regression model that models the original data as closely as possible with as few predictors as possible. We however want to smooth out noise from the original norm data, which can be due to the random sampling process or violations of representativeness. This is done through the ‘bestModel’ function. You can use this function in two different ways: If you specify \(R_{adjusted}^{2}\), then the regression function will be selected that meets this requirement with the smallest number of predictors. You can however also specify a fixed number of predictors. Then the model is selected that achieves the highest \(R_{adjusted}^{2}\) with this specification. To select the best model, cNORM uses the ‘regsubset’ function from the ‘leaps’ package. As we do not know beforehand, how well the data can be modeled, we start with the default values (k = 4 and \(R_{adjusted}^{2}\) = .99):

# If you only need the model, than use
# model <- bestModel(normData)
# Or just the convenience method that does everything at once
model <- cnorm(raw=elfe$raw, group=elfe$group)
#> Powers of location: k = 5
#> Powers of age:      t = 3
#> Multiple R2 between raw score and explanatory variable: R2 = 0.5129
#> Final solution: 3 terms
#> R-Square Adj. = 0.990838
#> Final regression model: raw ~ L3 + L1A1 + L3A3
#> Regression function: raw ~ -11.39665566 + (2.08886395e-05*L3) + (0.1649386974*L1A1) + (-5.892663055e-07*L3A3)
#> Raw Score RMSE = 0.68348
#> Use 'printSubset(model)' to get detailed information on the different solutions, 'plotPercentiles(model) to display percentile plot, plotSubset(model)' to inspect model fit.

Fine! The determined model already exceeds the predefined threshold of \(R_{adjusted}^{2}\) = .99 with only three predictors (plus intercept). The ‘bestModel’ function as well returns the coefficients and the complete regression formula, which - as was specified - captures more than 99% of the variance in the data set.

If you want to have a look at the selection procedure, all the information is available in ‘model$subsets’. The variable selection process per step is listed in ‘outmat’ and ‘which’. There, you can find the \(R^{2}\), \(R_{adjusted}^{2}\), \(C_p\) and \(BIC\):

#>        R2adj       BIC          CP       RSS      RMSE    DeltaR2adj
#> 1  0.9197525 -3518.209 14746.83265 5736.1349 2.0241638            NA
#> 2  0.9800576 -5461.138  2614.81241 1424.4767 1.0087038  6.030503e-02
#> 3  0.9908376 -6543.733   448.50221  653.9973 0.6834771  1.078002e-02
#> 4  0.9914116 -6628.074   333.95429  612.5836 0.6614830  5.740480e-04
#> 5  0.9917828 -6683.686   260.26527  585.6885 0.6467990  3.711756e-04
#> 6  0.9920204 -6718.526   213.45532  568.3446 0.6371502  2.376106e-04
#> 7  0.9922445 -6752.171   169.41665  551.9854 0.6279134  2.241163e-04
#> 8  0.9925535 -6802.851   108.45630  529.6133 0.6150571  3.089811e-04
#> 9  0.9926317 -6811.395    93.74201  523.6741 0.6115987  7.820923e-05
#> 10 0.9926689 -6812.246    87.24539  520.6549 0.6098331  3.720644e-05
#> 11 0.9927589 -6823.287    70.23231  513.8989 0.6058635  8.991488e-05
#> 12 0.9927927 -6823.616    64.42967  511.1263 0.6042270  3.387453e-05
#> 13 0.9928281 -6824.277    58.33021  508.2483 0.6025234  3.541165e-05
#> 14 0.9928384 -6820.050    57.25115  507.1542 0.6018746  1.027172e-05
#> 15 0.9928841 -6822.769    49.12978  503.5577 0.5997367  4.564866e-05
#> 16 0.9928940 -6818.497    48.12485  502.4899 0.5991005  9.954472e-06
#> 17 0.9929119 -6815.791    45.54800  500.8636 0.5981302  1.787355e-05
#> 18 0.9929725 -6821.589    34.47790  496.2193 0.5953506  6.063981e-05
#> 19 0.9930215 -6825.145    25.74684  492.4062 0.5930588  4.894852e-05
#> 20 0.9930324 -6821.116    24.56708  491.2763 0.5923780  1.096405e-05
#> 21 0.9930473 -6817.875    22.61991  489.8737 0.5915318  1.485024e-05
#> 22 0.9930457 -6811.317    23.94451  489.6337 0.5913868 -1.640457e-06
#> 23 0.9930504 -6806.050    24.00000  488.9428 0.5909694  4.766823e-06
#>               F            p nr
#> 1            NA           NA  1
#> 2  4228.4907653 0.000000e+00  2
#> 3  1644.6386423 0.000000e+00  3
#> 4    94.3089615 0.000000e+00  4
#> 5    64.0130593 2.553513e-15  5
#> 6    42.5096286 9.810242e-11  6
#> 7    41.2547575 1.829583e-10  7
#> 8    58.7590888 3.330669e-14  8
#> 9    15.7645352 7.540381e-05  9
#> 10    8.0544971 4.605090e-03 10
#> 11   18.2475156 2.071841e-05 11
#> 12    7.5236691 6.167729e-03 12
#> 13    7.8484272 5.157215e-03 13
#> 14    2.9879117 8.411115e-02 14
#> 15    9.8847574 1.701950e-03 15
#> 16    2.9387868 8.669971e-02 16
#> 17    4.4874050 3.432353e-02 17
#> 18   12.9252315 3.355968e-04 18
#> 19   10.6865636 1.105753e-03 19
#> 20    3.1715487 7.515144e-02 20
#> 21    3.9453954 4.719763e-02 21
#> 22    0.6749444 4.114753e-01 22
#> 23    1.9445051 1.634053e-01 23

Furthermore, information about the change of Radjusted and other information criteria (Mallow’s Cp or BIC) depending on the number of predictors (with fixed k) can also be graphically inspected. Please use the following command to do this:

plot(model, "subset", type = 0) 

The figure displays Radjusted2 as a function of the number of predictors by default. Alternatively, you can also plot log-transformed Mallow’s \(C_p\) (type = 1) and \(BIC\) (type = 2) as a function of \(R_{adjusted}^{2}\) or RMSE (type = 3) as a function of the number of terms.

plot(model, "subset", type = 1) 

The figure shows that the default value of \(R_{adjusted}^{2}\) = .99 is already achieved with only three predictors. The inclusion of further predictors only leads to small increases of \(R_{adjusted}^{2}\) or to small decreases of Mallow’s \(C_p\). Where the dots are close together, the inclusion of further predictors is of little use. To avoid over-fitting, a model with as few predictors as possible should therefore be selected from this area.

The model with three predictors seems to be suitable. Nevertheless, the model found in this way must still be tested for plausibility using the means described in Model Validation. Above all, it is necessary to determine the limits of model validity. If a model turns out to be suboptimal after this model check, \(R_{adjusted}^{2}\), the number of predictors or, if necessary, k should be chosen differently. What is more, you can use the cross validation function to get an impression on the quality of the norm sample and modeling process. The function determines RMSE for the raw score of the training and validation data (80% and 20% drawn from the data set), \(R^{2}\) for the norm scores, crossfit and norm score \(R^{2}\). Crossfit values below indicate an underfit and values greater 1 an overfit, with values between .9 and 1.1 being optimal:

# do a cross validation with 2 repetitions, restrict number of max terms to 10$data, max=10, repetitions = 2)

In this specific case, it might not be worthwhile to use more terms than 4, because \(R^{2}\) in the cross validated data set does not increase anymore.You should however as well inspect to percentile curves visually or use the ‘checkConsistency’ function on the final model to avoid intersecting percentile curves and use a model with more terms in that case.

2. Model Validation

From a mathematical point of view, the regression function represents a so-called hyperplane in three-dimensional space. If \(R^2\) is sufficiently high (e.g. \(R^2 > .99\)), this plane usually models the manifest data over wide ranges of the standardization sample very well. However, a Taylor polynomial, as used here, usually has a finite radius of convergence. This means that there are age or performance ranges for which the regression function no longer provides plausible values. With high R2, these limits of model validity are only reached at the outer edges of the age or performance range of the standardization sample or even beyond. Please note that such model limits occur not only because the method is not omnipotent, but also because the underlying test scales have only a limited validity range within which they can reliably map a latent ability to a meaningful numerical test score. In other words, the limits of model validity often show up at those points where the test has too strong floor or ceiling effects or where the standardization sample is too diluted.

Of course, norm tables and normal scores should generally only be issued within the validity range of the model and the test. It is therefore essential to determine the limits of model validity when applying cNORM (or any other procedure used to model normal scores). For this purpose, cNORM mainly provides graphical methods, which we present to you on this page. At this point, however, we would like to point out to the mathematically experienced users that it is also possible to approach the topic analytically. Since the regression equation is a polynomial of the nth degree that is very easy to handle from a mathematical point of view, it can be subjected to a conventional curve sketching. This makes it very easy to determine, for example, where extremes, turning points, saddle points, etc. occur or where the gradient has implausible values.

Below you will find three functions for checking the model fit graphically and making the limits of the model visible:


The following figure shows how well the model generally fits the manifest data:

# Plots the fitted and the manifest percentiles
# modeling already  displays the plot; you can call it
# directly with plot(results) as well
plot(model, "percentiles")

As the figure shows, the predicted percentiles run smoothly across all levels of the explanatory variable and are in good agreement with the original data. Small fluctuations between the individual groups are eliminated. It is important to ensure that the percentile lines do not intersect, since this would mean that different values of the latent person variable are assigned to one single raw score. The mapping of latent person variables to raw scores would no longer be biunique (=bijective) at this point, e.g. it would not be possible to distinguish between these different values of the latent variable by means of the test score. As already described above, intersecting percentiles predominantly occur when the regression model is extended to age or performance ranges that do not or only rarely occur in the standardization sample, or when the test shows strong floor or ceiling effects.

If you are not sure yet, which model to choose, you can display a series of plots:

# Displays a series of plots of the fitted and the manifest percentiles
# Here the maximum number of terms is fixed (optional)
plot(model, "series", end = 10)


In the next figure, the fitted and the manifest data are compared separately for each (age) group:

plot(model, "raw")

The adjustment is particularly good if all points are as close as possible to the bisecting line. However, it must be noted that deviations in the extremely upper, but particularly in the extremely lower performance range often occur because the manifest data in these areas are also associated with low reliability.


This function does the same as plotRaw, but instead uses norm values. Please note, that it might take some time depending on the sample size:

plot(model, "norm")

The plot can be split by group as well:

plot(model, "norm", group = "group")    # specifies the grouping variable


To check whether the mapping between latent person variables and test scores is biunique, the regression function can be searched numerically within each group for bijectivity violations using the ‘checkConsistency’ function. In addition, it is also possible to plot the first partial derivative of the regression function to l and search for negative values. This can be done in the following way:

plot(model, "derivative", minAge=1, maxAge=6, minNorm=20, maxNorm=80)
#> Horizontal and vertical extrapolation detected. Be careful using age groups and extreme norm scores outside the original sample.
#> The original data for the regression model spanned from age 2 to 5, with a norm score range from 21.93 to 78.07. The raw scores range from 0 to 28. Coefficients from the 1 order derivative function:
#>            L2            A1          L2A3 
#>  6.266592e-05  1.649387e-01 -1.767799e-06

# if parameters on age an norm are not specified, cnorm plots within
# the ranges of the current dataset

In this figure, we have extended both the age and the performance range beyond the limits of the standardization sample in order to better represent and check the limits of the model validity. (Please remember that the age variable in this norming sample comprises the values 2 to 5 and that 200 children per age group were examined.) As you can see, the first partial derivative of the regression function to l is only negative in the upper age and performance range. This does not mean that the modeling has failed, but that the test scale loses its ability to differentiate in this measuring range.

When, at the end of the modeling process, norm tables are generated, the identified limits of model validity must be respected. Or to put it in other words: Normal scores should only be issued for the valid ranges of the model.

3. Generating Norm Tables

In addition to the pure modeling functions, cNORM also contains functions for generating norm tables, retrieving the normal score for a specific raw score and vice versa or for the visualization of norm curves. These functions are described below.


The first function ‘getNormCurve’ returns the fitted raw scores for a certain normal score (e.g., T = 50) across different age groups. The parameter ‘step’ specifies the distance between two age groups. If no further specifications are made, the output is limited to actually occurring raw values and age groups.

getNormCurve(50, model, minAge = 2, maxAge = 5, step = 0.25, minRaw = 0, maxRaw = 28)


Plots the fitted raw scores for pre-specified normal scores (e.g., T = 30, 40, 50, 60, 70) across age.

plot(model, "curves", normList = c(30, 40, 50, 60, 70), minAge = 2, maxAge = 5, step = 0.1, minRaw = 0, maxRaw = 28)


The ‘predictNorm’ function returns the normal score for a specific raw score (e.g., raw = 15) and a specific age (e.g., a = 4.7). The normal scores can be limited to a minimum and maximum value in order to take into account the limits of model validity.

predictNorm(15, 4.7, model, minNorm = 25, maxNorm = 75)
#> [1] 36.59881


The ‘predictRaw’ function returns the predicted raw score for a specific normal score (e.g., T = 55) and a specific age (e.g., a = 4.5).

predictRaw(55, 4.5, model)
#> [1] 23.9672

… or a matrix, if you like …

predictRaw(c(45, 50, 55), c(2.5, 3, 3.5, 4, 4.5), model)
#>          2.5        3      3.5        4      4.5
#> 45  8.223411 11.32373 14.18241 16.75919 19.01378
#> 50 10.680851 13.96646 16.92060 19.48803 21.61352
#> 55 13.225900 16.64652 19.62596 22.09070 23.96720


The ‘normTable’ function returns the corresponding raw scores for a specific age (e.g., a = 3) or vector of ages and a pre-specified series of normal scores. The parameter ‘step’ specifies the distance between two normal scores.

# Generate norm table for grade 3
normTable(3, model, minRaw = 0, maxRaw = 28, minNorm=30, maxNorm=70, step = 1)

# You as well can generate multiple tables at once, by providing a vector for the age groups,
# here for grade 3, 3.5, 4 and 4.5
normTables -> normTable(c(3.5, 4.5), model, minRaw = 0, maxRaw = 28, 
              minNorm=30, maxNorm=70, step = 1)

# In case, you want to use the data for producing paper and pencil tests, you might
# want to export the data. Here, this is demonstratet with an export to Excel:
# Install psych package
install.packages("openxlsx", dependencies = TRUE)
write.xlsx(normTables, file="normtables.xlsx")

Norm tables, in which the raw score or the range of raw scores is given for a certain normal score, are usually needed if one has several test scales, all of which are to be converted into the same type of normal scale. Please note that as a test designer you cannot use this function directly to generate norm tables. Instead you have to convert the table to a suitable form first. Remember that raw scores are usually integers. Therefore, the normal score series should not start at an integer value, but e.g. at 30.5. If a step size of 1 is selected, normal score intervals with integers as interval centers are generated. In the example shown, the T-score interval [30.5; 31.5] contains only one single integer raw score, namely 4, which would therefore be assigned to the T-score 31. A whole range of raw scores (or no raw score at all) can thus be assigned to a particular integer normal score.


The function ‘rawTable’ is similar to ‘normTable’, but reverses the assignment: The normal scores are assigned to a pre-specified series of raw scores at a certain age. This requires an inversion of the regression function, which is determined numerically.

# Generate raw table for grade 3.5
rawTable(3.5, model, minRaw = 0, maxRaw = 28, minNorm = 25, maxNorm = 75, step = 1)

# Generate raw table for grade 3.5, 3.6, 3.7
rawTable(c(3.5, 3.6, 3.7), model, minRaw = 0, maxRaw = 28, minNorm = 25, maxNorm = 75, step = 1)

You need these kind of tables if you want to determine the exact percentile or the exact normal score for all occurring raw scores. With higher precision and smaller increments, this function becomes computationally intensive.


Please visit for further examples with the inbuilt data sets, especially with respect to continuous explanatory variables and non standard application cases.