NAG Library Routine Document
g22ycf (lm_design_matrix)
1
Purpose
g22ycf generates a design matrix from a data matrix and model description.
2
Specification
Fortran Interface
Subroutine g22ycf ( 
hform, hddesc, dat, lddat, sddat, hxdesc, x, ldx, sdx, mx, ifail) 
Integer, Intent (In)  ::  lddat, sddat, ldx, sdx  Integer, Intent (Inout)  ::  ifail  Integer, Intent (Out)  ::  mx  Real (Kind=nag_wp), Intent (In)  ::  dat(lddat,sddat)  Real (Kind=nag_wp), Intent (Inout)  ::  x(ldx,sdx)  Type (c_ptr), Intent (In)  ::  hform, hddesc  Type (c_ptr), Intent (Inout)  ::  hxdesc 

C Header Interface
#include <nagmk26.h>
void 
g22ycf_ (void **hform, void **hddesc, const double dat[], const Integer *lddat, const Integer *sddat, void **hxdesc, double x[], const Integer *ldx, const Integer *sdx, Integer *mx, Integer *ifail) 

3
Description
g22ycf generates a design matrix from a data matrix and a model description. Design matrices encapsulate the observed values of the independent variables and the required model in a form that can be used by many of the model fitting routines available in the NAG Library, for example those in
Chapter G02.
3.1
Notation
Let $D$ denote a data matrix with $n$ observations on ${m}_{d}$ independent variables, denoted by ${V}_{j}$, for $\mathit{j}=1,2,\dots ,{m}_{d}$. If ${V}_{j}$ is a categorical variable, let ${L}_{j}$ denote the number of levels associated with it. If ${V}_{j}$ is a binary, ordinal or continuous variable, let ${L}_{j}=1$.
Let ${V}_{ji}$ denote the $i$th value of ${V}_{j}$.
Let $\mathcal{M}$ denote a model made up of one or more terms, denoted by ${T}_{i}$. Each term consists of either a main effect or an interaction and hence can be described using one or more variable names ${V}_{j}$ and the interaction operator ‘$.$’. The operator ‘$+$’ is used to denote the addition of a term to the model. Therefore, $\mathcal{M}={T}_{1}+{T}_{2}+{T}_{3}={V}_{1}+{V}_{2}+{V}_{1}.{V}_{2}$ denotes a model with three terms, the first two terms being the main effects for variables ${V}_{1}$ and ${V}_{2}$ and the last term the interaction between them. For simplicity we reorder the terms of the model by the number of variables in them, so main effects come first, then twoway interactions, then threeway interactions etc. By default it is assumed that the model $\mathcal{M}$ contains a mean effect (or intercept term), if the mean effect is excluded, this will be denoted by ‘$1$’, so $\mathcal{M}={T}_{1}$ is a model with one term and a mean effect and $\mathcal{M}={T}_{1}1$ is the same model with the mean effect dropped.
g22ycf generates an $n$ by ${m}_{x}$ design matrix, $X$, from $D$ and $\mathcal{M}$.
3.2
Dummy Variables
When constructing a design matrix, we cannot work directly with categorical variables. Categorical variables must first be recoded into dummy variables. A categorical variable
${V}_{j}$ requires
${L}_{j}$ dummy variables. Let
${\mathcal{D}}^{j}$ denote an
$n\times {L}_{j}$ matrix of dummy variables for
${V}_{j}$ defined as
where
${\mathcal{D}}_{l}^{j}$ is the
$l$th column of
${\mathcal{D}}^{j}$ and
${\mathcal{D}}_{li}^{j}$ is the
$i$th element of
${\mathcal{D}}_{l}^{j}$.
For a binary, ordinal or continuous variable, ${\mathcal{D}}_{1i}^{j}={V}_{ji}$.
3.3
Full Design Matrix
Given a model,
$\mathcal{M}$, and the matrices of dummy variables constructing the full design matrix
${X}_{F}$ is trivial. Each term is processed in order and
1. 
If term $i$ is a main effect, that is ${T}_{i}={V}_{j}$ for some $j$, ${\mathcal{D}}^{j}$ is copied into ${X}_{F}$. 
2. 
If term $i$ is a twoway interaction, that is ${T}_{i}={V}_{j}.{V}_{k}$, for some $j\ne k$, then
(i) 
Loop over ${l}_{j}=1,2,\dots {L}_{j}$. 
(ii) 
Loop over ${l}_{k}=1,2,\dots {L}_{k}$. 
(iii) 
Add a column to ${X}_{F}$ corresponding to the elementwise product of ${\mathcal{D}}_{{l}_{j}}^{j}$ and ${\mathcal{D}}_{{l}_{k}}^{k}$. 

3. 
Higher interaction terms are handled in a similar manner as the twoway interactions by adding columns constructed from multiplying all combinations of the columns of the corresponding $\mathcal{D}$s that correspond to the variables involved. In all cases, the variables towards the right hand side of a term are iterated over the quickest. 
3.4
Contrasts
Using the full design matrix ${X}_{F}$ in an analysis can result in an overparameterized model. This is due to ${X}_{F}$ often not being of full rank as the sum of all the dummy variables for a particular variable is a vector of ones. This source of overparameterization can be alleviated by using a design matrix $X$ where (some) dummy variables are replaced by contrasts. For a categorical variable ${V}_{j}$ the contrasts are a set of ${L}_{j}1$ functionally independent linear combinations of the dummy variables.
Whilst the choice of contrasts used in term ${T}_{i}$ will affect the individual model coefficients (parameters), it has no effect on the overall contribution of ${T}_{i}$.
For a given variable ${V}_{j}$, the contrasts can be represented by an ${L}_{j}$ by ${L}_{j}1$ matrix, ${C}_{j}$. The rows of ${C}_{j}$ correspond to a particular value of ${V}_{j}$ and the columns correspond to the values to use in the design matrix.
Six types of contrast are available in
g22ycf; two types of treatment contrasts, two types of sum contrasts, Helmert contrasts and polynomial contrasts. Unless specified otherwise, the contrasts used by
g22ycf are treatment contrasts relative to the first level. See the description of the optional parameter
Contrast in
g22yaf for ways of changing the contrasts used.
3.4.1
Treatment Contrasts
Treatment contrasts are taken relative to either the first or last level of the variable. For example, if
${L}_{j}=4$,
would be the contrast matrix for
${V}_{j}$ using treatment contrasts relative to the first level. The contrast matrix obtained when using treatment contrasts relative to the last level is similar, but the row of zeros appears at the bottom and all other rows are shifted up one.
Strictly speaking, the term contrast implies that each row in the contrast matrix sums to zero. That is not the case for treatment contrasts, however they are included as this coding is commonly used in practice.
3.4.2
Sum Contrasts
Sum contrasts are similar to treatment contrasts and again can be taken relative to the first or last level of the variable. Unlike treatment contrasts, sum contrasts effectively constrain the coefficients related to the variable to sum to zero. For example, if
${L}_{j}=4$,
would be the contrast matrix for
${V}_{j}$ using treatment contrasts relative to the last level. The contrast matrix obtained when using treatment contrasts relative to the first level is similar, but the row of
$1$s appears at the top and all other rows are shifted down one.
3.4.3
Helmert Contrasts
With Helmert contrasts level
$l$ of the variable is compared with the average effect of all previous levels. For example, if
${L}_{j}=4$,
would be the contrast matrix for
${V}_{j}$ using Helmert contrasts.
3.4.4
Polynomial Contrasts
With polynomial contrasts the entries in the columns of
${C}_{j}$ correspond in linear, quadratic, cubic, quartic, etc. terms to a hypothetical underlying numeric variable that takes equally spaced values at each level. For example, if
${L}_{j}=4$,
would be the contrast matrix for
${V}_{j}$ using polynomial contrasts.
3.4.5
When Contrasts Can Be Used
Depending on the specifics of the model,
$\mathcal{M}$, it may not be possible to always replace the
${L}_{j}$ dummy variables with
${L}_{j}1$ contrasts for all variables in all terms and retain the same model. A simple example of this is a data matrix,
$D$, with four observations and two variables which have two and three levels respectively. This data matrix might look something like:
For the sake of argument, assume that our model contains the main effect for each variable, but does not contain a mean effect (or intercept term). So using the notation established earlier,
$\mathcal{M}={V}_{1}+{V}_{2}1$. The full design matrix,
${X}_{F}$, for this data matrix and model would be
However, ${X}_{F}$ is not of full rank (and hence $\mathcal{M}$ is overparameterized) because the sum of the first two columns is a vector of ones as is the sum of the last three columns.
In order to alleviate this we might try constructing
${X}_{C}$ where the dummy variables have been replaced by contrasts. Assuming treatment contrasts, relative to the first level, we would have
However, using
${X}_{C}$ makes an implicit assumption that the expected value of the dependent variable (the quantity being modelled) is zero when
${V}_{1}=1$ and
${V}_{2}=1$. This assumption was not made when we used
${X}_{F}$ and hence the two design matrices are not equivalent. One solution would be to use dummy variables for
${V}_{1}$ and contrasts for
${V}_{2}$, which would result in a design matrix,
$X$ of
Using $X$ would give an equivalent model to using ${X}_{F}$.
The algorithm used by g22ycf to decide which variables, in which terms, can be coded as contrasts and which need to be coded as dummy variables is described below.
Suppose ${V}_{j}$ is any variable that appears in term ${T}_{i}$, let ${T}_{i\left(j\right)}$ denote the term obtained by dropping ${V}_{j}$ from ${T}_{i}$. For example, if ${T}_{3}={V}_{1}.{V}_{2}.{V}_{3}$, ${T}_{3\left(2\right)}={V}_{1}.{V}_{3}$. In this context, the empty term is taken to be the mean effect (or intercept term). We say that ${T}_{i\left(j\right)}$ appears in $\mathcal{M}$ if there exists a term ${T}_{k}$, $k<i$, that contains all of the variables appearing in ${T}_{i\left(j\right)}$. In most cases ${T}_{k}={T}_{i\left(j\right)}$, but this is not required. Note, as stated earlier, the terms in $\mathcal{M}$ are ordered by the number of variables in them.
A variable, ${V}_{j}$ in term ${T}_{i}$ is coded by contrasts if ${T}_{i\left(j\right)}$ appears in $\mathcal{M}$ and by dummy variables otherwise. It is therefore possible for variable ${V}_{j}$ to be coded by contrasts in some terms and dummy variables in others within the same $X$.
The above rule assumes the presence of a mean effect. If no such effect is present in the model, the main effect of the first categorical variable is coded by dummy variables to compensate. If no main effects appear in the model, the warning ${\mathbf{ifail}}={\mathbf{14}}$ is returned.
A longer description and informal proof that the resulting
$X$ is a suitable design matrix for the model of interest can be found in chapter two of
Chambers and Hastie (1992).
3.5
Mean Effect
The mean effect (or intercept term) is included in a design matrix by adding a column of ones as the first column of
$X$. However, many model fitting routines in the NAG Library handle the mean effect as a special case and do not require it to be explicitly added to the design matrix. Therefore, by default,
g22ycf does not explicitly add the mean effect to the design matrix. This behaviour can be changed via the optional parameter
Explicit Mean in
g22yaf.
4
References
Chambers J M and Hastie T J (1992) Statistical Models in S Wadsworth and Brooks/Cole Computer Science Series
5
Arguments
 1: $\mathbf{hform}$ – Type (c_ptr)Input

On entry: a G22 handle to the internal data structure containing a description of the model
$\mathcal{M}$ as returned in
hform by
g22yaf.
 2: $\mathbf{hddesc}$ – Type (c_ptr)Input

On entry: a G22 handle to the internal data structure containing a description of the data matrix,
$D$ as returned in
hddesc by
g22ybf.
 3: $\mathbf{dat}\left({\mathbf{lddat}},{\mathbf{sddat}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the data matrix,
$D$. By default
${D}_{ij}$, the
$\mathit{i}$th value for the
$\mathit{j}$th variable, for
$\mathit{i}=1,2,\dots ,n$ and
$\mathit{j}=1,2,\dots ,{m}_{d}$, should be supplied in
${\mathbf{dat}}\left(i,j\right)$.
If the optional parameter
Storage Order, described in
g22ybf, is set to
$\mathrm{VAROBS}$,
${D}_{ij}$ should be supplied in
${\mathbf{dat}}\left(j,i\right)$.
 4: $\mathbf{lddat}$ – IntegerInput

On entry: the first dimension of the array
dat as declared in the (sub)program from which
g22ycf is called.
Constraints:
 if the optional parameter Storage Order, described in g22ybf, is set to $\mathrm{VAROBS}$, ${\mathbf{lddat}}\ge {m}_{d}$;
 otherwise ${\mathbf{lddat}}\ge n$.
 5: $\mathbf{sddat}$ – IntegerInput

On entry: the second dimension of the array
dat as declared in the (sub)program from which
g22ycf is called.
Constraints:
 if the optional parameter Storage Order, described in g22ybf, is set to $\mathrm{VAROBS}$, ${\mathbf{sddat}}\ge n$;
 otherwise ${\mathbf{sddat}}\ge {m}_{d}$.
 6: $\mathbf{hxdesc}$ – Type (c_ptr)Input/Output

On entry: must be set to
c_null_ptr.
As an alternative an existing G22 handle may be supplied in which case this routine will destroy the supplied G22 handle as if
g22zaf had been called.
On exit: holds a G22 handle to the internal data structure containing a description of the design matrix,
$X$. You
must not change the G22 handle other than through the routines in
Chapter G22.
 7: $\mathbf{x}\left({\mathbf{ldx}},{\mathbf{sdx}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the design matrix,
$X$. By default
${X}_{ij}$, the
$\mathit{i}$th value for the
$\mathit{j}$th column, for
$\mathit{i}=1,2,\dots ,n$ and
$\mathit{j}=1,2,\dots ,{m}_{x}$, is returned in
${\mathbf{x}}\left(i,j\right)$ If the optional parameter
Storage Order, described in
g22yaf, is set to
$\mathrm{VAROBS}$,
${X}_{ij}$ is returned in
${\mathbf{x}}\left(j,i\right)$.
If
ldx or
sdx are too small to hold
x, the number of columns required to hold the design matrix is returned in
mx.
Under some conditions it is possible to use the data matrix in place of the design matrix. Specifically, if
$D$ has no categorical variables,
$\mathcal{M}$ has only main effects and either has no mean effect or the mean effect does not need to be explicitly added to the design matrix. If
ldx or
sdx are too small under such circumstances,
${\mathbf{ifail}}={\mathbf{71}}$ is returned and
hxdesc is set up in such a way as to allow
dat to be used as the design matrix.
 8: $\mathbf{ldx}$ – IntegerInput

On entry: the first dimension of the array
x as declared in the (sub)program from which
g22ycf is called.
Constraints:
 if the optional parameter Storage Order, described in g22yaf, is set to $\mathrm{VAROBS}$, ${\mathbf{ldx}}\ge {m}_{x}$;
 otherwise ${\mathbf{ldx}}\ge n$.
 9: $\mathbf{sdx}$ – IntegerInput

On entry: the second dimension of the array
x as declared in the (sub)program from which
g22ycf is called.
Constraints:
 if the optional parameter Storage Order, described in g22yaf, is set to $\mathrm{VAROBS}$, ${\mathbf{sdx}}\ge n$;
 otherwise ${\mathbf{sdx}}\ge {m}_{x}$.
 10: $\mathbf{mx}$ – IntegerOutput

On exit: the minimum number of columns required to hold the design matrix.
In most cases
${\mathbf{mx}}={m}_{x}$. The one exception is when
${\mathbf{ifail}}={\mathbf{71}}$, that is the size of
x was too small but the data matrix given in
dat can be used as the design matrix. In this case
mx holds the number of columns that would be required if only the relevant parts of
dat were copied into a new array.
 11: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=11$

hform has not been initialized or is corrupt.
 ${\mathbf{ifail}}=12$

hform is not a G22 handle as generated by
g22yaf.
 ${\mathbf{ifail}}=13$

A variable name used when creating
hform is not present in
hddesc.
Variable name:
$\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=14$

The model contains categorical variables, but no intercept or main effects terms have been requested.
Please check the design matrix returned matches the model you require.
 ${\mathbf{ifail}}=21$

hddesc has not been initialized or is corrupt.
 ${\mathbf{ifail}}=22$

hddesc is not a G22 handle as generated by
g22ybf.
 ${\mathbf{ifail}}=31$

On entry, column
$j$ of the data matrix,
$D$, is not consistent with information supplied in
hddesc,
$j=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=41$

On entry, $n=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{lddat}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lddat}}\ge n$.
 ${\mathbf{ifail}}=42$

On entry, ${m}_{d}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{lddat}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lddat}}\ge {m}_{d}$.
 ${\mathbf{ifail}}=51$

On entry, ${m}_{d}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{sddat}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{sddat}}\ge {m}_{d}$.
 ${\mathbf{ifail}}=52$

On entry, $n=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{sddat}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{sddat}}\ge n$.
 ${\mathbf{ifail}}=61$

On entry,
hxdesc is not
c_null_ptr or a recognised G22 handle.
 ${\mathbf{ifail}}=71$

On entry, the size of
x is too small to hold the design matrix.
dat can be used instead.
 ${\mathbf{ifail}}=81$

On entry, $n=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{ldx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ldx}}\ge n$.
 ${\mathbf{ifail}}=82$

On entry, ${m}_{x}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{ldx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ldx}}\ge {m}_{x}$.
 ${\mathbf{ifail}}=91$

On entry, ${m}_{x}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{sdx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{sdx}}\ge {m}_{x}$.
 ${\mathbf{ifail}}=92$

On entry, $n=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{sdx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{sdx}}\ge n$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
Not applicable.
8
Parallelism and Performance
g22ycf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g22ycf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
g22ydf can be used to obtain labels for the columns of the design matrix
$X$.
Many of the analysis routines that require a design matrix to be supplied allow submodels to be defined through the use of a vector of ones or zeros indicating whether a column of
$X$ should be included or excluded from the analyses (see for example
isx in
g02daf or
g02gaf). This allows nested models to be fit without having to reconstruct the design matrix for each analysis.
g22ydf offers a mechanism for constructing these vectors using submodels specified using
g22yaf.
10
Example
This example creates and outputs two design matrices for a simple linear regression model. The first design matrix uses sum contrasts for all variables and the second uses a combination of polynomial and Helmert contrasts. Column labels are generated using
g22ydf.
See also the examples for
g22yaf,
g22ybf and
g22ydf.
10.1
Program Text
Program Text (g22ycfe.f90)
10.2
Program Data
Program Data (g22ycfe.d)
10.3
Program Results
Program Results (g22ycfe.r)
11
Optional Parameters
As well as the optional parameters common to all G22 handles described in
g22zmf and
g22znf, a number of additional optional parameters can be specified for a G22 handle holding the description of a design matrix as returned by
g22ycf in
hxdesc.
The value of an optional parameter can be queried using
g22znf.
The remainder of this section can be skipped if you wish to use the default values for all optional parameters.
The following is a list of the optional parameters available. A full description of each optional parameter is provided in
Section 11.1.
11.1
Description of the Optional Parameters
For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
 a parameter value,
where the letters $a$, $i$ and $r$ denote options that take character, integer and real values respectively;
Keywords and character values are case and white space insensitive.
This optional parameter returns a verbose formula string describing the model, $\mathcal{M}$, used to create the design matrix. This formula will only contain variable names, the operators ‘$+$’ and ‘$.$’ and any contrast identifiers present.
Min Number of Columns  $i$  Read Only 
This optional parameter returns the minimum number of columns required to hold the design matrix,
$X$. In most cases
${\mathbf{Min\; Number\; of\; Columns}}={\mathbf{Number\; of\; Columns}}$. The one exception is when
${\mathbf{ifail}}={\mathbf{71}}$, that is the size of
x was too small but the data matrix given in
dat can be used as the design matrix. In this case,
${\mathbf{Number\; of\; Columns}}={m}_{x}={m}_{d}$ and
${\mathbf{Min\; Number\; of\; Columns}}$ holds the number of columns that would be required if only the relevant parts of
dat were copied into a new array.
Number of Columns  $i$  Read Only 
This optional parameter returns ${m}_{x}$, the number of columns in the design matrix.
Number of Observations  $i$  Read Only 
This optional parameter returns $n$, the number of observations in the design matrix.
Storage Order  $a$  Read Only 
This optional parameter returns how the design matrix,
$X$, is stored in
x.
If ${\mathbf{Storage\; Order}}=\mathrm{OBSVAR}$, ${X}_{ij}$, the value for the $j$th variable of the $i$th observation of the design matrix is stored in ${\mathbf{x}}\left(i,j\right)$.
If ${\mathbf{Storage\; Order}}=\mathrm{VAROBS}$, ${X}_{ij}$, the value for the $j$th variable of the $i$th observation of the design matrix is stored in ${\mathbf{x}}\left(j,i\right)$.
It should be noted that
Storage Order is not writeable. If you wish to change the storage order of the design matrix you need to change
Storage Order in
hform as described in
Section 11 in
g22yaf prior to calling
g22ycf.