RHPSS
Computes the matrix of sums of squares and crossproducts for the multivariate general linear hypothesis H BU = G given the coefficient estimates
and the R matrix.
Required Arguments
H — NH by NCOEF matrix H with each row corresponding to a row in the hypothesis and containing the constants that specify an estimable linear combination of the regression coefficients. (Input)
B — NCOEF by NDEP matrix
containing a least-squares solution for the regression coefficients. (Input)
G — Matrix containing the null hypothesis values. (Input)
If NU = 0, then G is NH by NDEP; otherwise, G is NH by NU.
R — NCOEF by NCOEF upper triangular matrix containing the R matrix. (Input)
The R matrix can come from a regression fit based on a QR decomposition of the matrix of regressors or based on a Cholesky factorization RTR of the matrix of sums of squares and crossproducts of the regressors. Elements to the right of a diagonal element of R that is zero must also be zero. A zero row indicates a nonfull rank model. For an R matrix that comes from a regression fit with linear equality restrictions on the parameters, each row of R corresponding to a restriction must have a corresponding diagonal element that is negative. The remaining rows of R must have positive diagonal elements. Only the upper triangle of R is referenced.
SCPH — Matrix containing sums of squares and crossproducts attributable to the hypothesis. (Output)
If NU = 0, SCPH is a NDEP by NDEP matrix, otherwise, SCPH is a NU by NU matrix.
Optional Arguments
NH — Number of rows in the hypothesis. (Input)
Default: NH = size (H,1).
NCOEF — Number of regression coefficients in the model. (Input)
Default: NCOEF = size (H,2).
LDH — Leading dimension of H exactly as specified in the dimension statement of the calling program. (Input)
Default: LDH = size (H,1).
NDEP — Number of dependent (response) variables. (Input)
Default: NDEP = size (B,2).
LDB — Leading dimension of B exactly as specified in the dimension statement in the calling program. (Input)
Default: LDB = size (B,1).
NU — U matrix option. (Input)
For positive NU, NU is the number of linear combinations of the dependent variables to be considered. If NU = 0, the hypothesis is HB = G, i.e., U is automatically taken to be the identity. NU must be less than or equal to NDEP.
Default: NU = 0.
U — NDEP by NU matrix U in test H BU = G. (Input, if NU is positive)
If NU = 0, U is not referenced and can be a 1 x 1 array.
Default: U is a 1x 1 array.
LDU — Leading dimension of U exactly as specified in the dimension statement in the calling program. (Input)
Default: LDU = size (U, 1,).
LDG — Leading dimension of G exactly as specified in the dimension statement in the calling program. (Input)
Default: LDG = size (G,1).
LDR — Leading dimension of R exactly as specified in the dimension statement in the calling program. (Input)
Default: LDR = size (R,1).
DFH — Degrees of freedom for SCPH. (Output)
DFH equals the rank of H.
LDSCPH — Leading dimension of SCPH exactly as specified in the dimension statement in the calling program. (Input)
Default: LDSCPH = size (SCPH,1).
FORTRAN 90 Interface
Generic: CALL RHPSS (H, B, G, R, SCPH [, …])
Specific: The specific interface names are S_RHPSS and D_RHPSS.
FORTRAN 77 Interface
Single: CALL RHPSS (NH, NCOEF, H, LDH, NDEP, B, LDB, NU, U, LDU, G, LDG, R, LDR, DFH, SCPH, LDSCPH)
Double: The double precision name is DRHPSS.
Description
Routine RHPSS computes the matrix of sums of squares and crossproducts for the general linear hypothesis H BU = G for the multivariate general linear model Y = XB + E with possible linear equality restrictions AB = Z. (See the chapter introduction for a description of the Multivariate General Linear Model.) Routines RGLM, RGIVN, RLEQU, and RCOV can be used to compute the fit of the general linear model prior to invoking RHPSS. The R matrix and 
from any of those routines are required for input to RHPSS.
The rows of H must be linear combinations of the rows of R, i.e., HB = G must be completely testable. If the hypothesis is not completely testable, Routine CESTI can be used to construct an equivalent completely testable hypothesis.
Computations are based on an algorithm discussed by Kennedy and Gentle (1980, page 317) that is extended by Sallas and Lionti (1988) for multivariate nonfull rank models with possible linear equality restrictions. The algorithm is as follows:
1. Form
2. Find C as the solution of RTC = HT using routine GIRTS (IMSL MATH/LIBRARY). If the equations are declared inconsistent within a computed tolerance, an error message with code 1 is issued that the hypothesis is not completely testable.
3. For all rows of R corresponding to restrictions, i.e., containing negative diagonal elements from a restricted least-squares fit using RLEQU, zero out the corresponding rows of C, i.e., form DC.
4. Decompose DC using Householder transformations and column pivoting to yield a square, upper triangular matrix T with diagonal elements of nonincreasing magnitude and permutation matrix P such that
where Q is an orthogonal matrix.
5. Determine the rank of T, say r. If t11 = 0, then r = 0. Otherwise, the rank of T is r if
where ɛ = 10.0 * AMACH(4). Then, zero out all rows of T below row r. Set the degrees of freedom for the hypothesis, output in DFH, to r.
6. Find V as a solution to T TV = PTW using routine GIRTS. If the equations are inconsistent, an error message with code 2 is issued that the hypothesis is inconsistent within a computed tolerance, i.e., the linear system
H BU = G
AB = Z
does not have a solution for B.
7. Form VTV, which is the required matrix of sum of squares and crossproducts output in SCPH.
In general, the two errors with code 1 and 2 are serious user errors that require the user to correct the hypothesis before any meaningful sums of squares from this routine can be computed. However, in some cases, the user may know the hypothesis is consistent and completely testable, but the checks in RHPSS are too tight. For this reason, RHPSS continues with the computations.
Routine RHPSS gives a matrix of sums of squares and crossproducts that could also be obtained from separate fittings of the two models
|
Y* = XB* + E* |
|
|
AB* = Z* |
(1) |
|
HB* = G |
|
and
|
Y* = XB* + E* |
|
|
AB* = Z* |
(2) |
where Y* = YU, B* = BU, E* = EU, and Z* = ZU. The error sum of squares and crossproduct matrix for (1) minus that for (2) is the matrix of sum of squares and crossproducts output in SCPH. Note that this approach avoids entirely the question of testability.
Comments
1. Workspace may be explicitly provided, if desired, by use of R2PSS/DR2PSS. The reference is:
CALL R2PSS (NCOEF, NH, H, LDH, NDEP, B, LDB, NU, U, LDU, G, LDG, R, LDR, DFH, SCPH, LDSCPH, IWK, WK)
The additional arguments are as follows:
IWK — Work vector of length NH.
WK — Work vector of length
NH * (NDEP + NCOEF + max(NCOEF, NH) + 3) + NU * NDEP ‑ 1.
2. Informational errors
|
Type |
Code |
Description |
|
3 |
1 |
The hypothesis is not completely testable. Each row of H must be in the space spanned by the rows of R. |
|
3 |
2 |
The hypothesis is inconsistent. The linear system HB U = G combined with any restrictions from a regression fit with linear equality restrictions must have a solution for B. |
3. 
where (CTDC)− is a generalized inverse of CTDC, C is a solution to RTC = HT, and D is a diagonal matrix with
Examples
Example 1
A two-way analysis-of-variance model is fitted to balanced data discussed by Snedecor and Cochran (1967, Table 12.5.1, page 347). The responses are the weight gains (in grams) of rats fed diets varying in two components-level of protein and source of protein. The model is
yijk = μ + α i + βj + γij + ɛ ijki = 1, 2; j = 1, 2, 3; k = 1, 2, …, 10
where
The model is fitted using routine RGLM. Next, the sum of squares for interaction
is computed using RHPSS. Finally, the F statistic is computed along with the associated p‑value using routine FDF (see Chapter 17, "Probability Distribution Functions and Inverses").
USE RHPSS_INT
USE RGLM_INT
USE UMACH_INT
USE FDF_INT
IMPLICIT NONE
INTEGER LDB, LDG, LDH, LDR, LDSCPE, LDSCPH, LDU, LDX, LINDEF, &
MAXB, NCOL, NDEP, NEF, NH, NROW, MAXCL, NCLVAR, J
PARAMETER (NDEP=1, LINDEF=4, MAXB=12, MAXCL=5, NCLVAR=2, NCOL=3, &
NEF=3, NH=2, NROW=60, LDB=MAXB, LDG=NH, LDH=NH, &
LDR=MAXB, LDSCPE=NDEP, LDSCPH=NDEP, LDX=NROW)
!
INTEGER INDCL(NCLVAR), INDDEP(NDEP), INDEF(LINDEF), INTCEP,&
IRANK, IRBEF(NEF+1), NCOEF, NOUT, NVEF(NEF)
REAL B(LDB,NDEP), DFE, DFH, F, G(LDG,NDEP), H(LDH,MAXB), &
PVALUE, R(LDR,MAXB), SCPE(LDSCPE,NDEP), &
SCPH(LDSCPH,NDEP),X(LDX,NCOL), XMAX(MAXB), &
XMIN(MAXB)
!
DATA X/73.0, 102.0, 118.0, 104.0, 81.0, 107.0, 100.0, 87.0, &
117.0, 111.0, 98.0, 74.0, 56.0, 111.0, 95.0, 88.0, 82.0, &
77.0, 86.0, 92.0, 94.0, 79.0, 96.0, 98.0, 102.0, 102.0, &
108.0, 91.0, 120.0, 105.0, 90.0, 76.0, 90.0, 64.0, 86.0, &
51.0, 72.0, 90.0, 95.0, 78.0, 107.0, 95.0, 97.0, 80.0, &
98.0, 74.0, 74.0, 67.0, 89.0, 58.0, 49.0, 82.0, 73.0, 86.0, &
81.0, 97.0, 106.0, 70.0, 61.0, 82.0, 30*1.0, 30*2.0, &
10*1.0, 10*2.0, 10*3.0, 10*1.0, 10*2.0, 10*3.0/
DATA INDCL/2, 3/, NVEF/1, 1, 2/, INDEF/2, 3, 2, 3/, INDDEP/1/
DATA (H(1,J),J=1,MAXB)/6*0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0/
DATA (H(2,J),J=1,MAXB)/6*0.0, 1.0, 0.0, -1.0, -1.0, 0.0, 1.0/
DATA G/2*0.0/
!
CALL RGLM (X, INDCL, NVEF, INDEF, NDEP, INDDEP, MAXCL, B, &
IRBEF=IRBEF, R=R, DFE=DFE, SCPE=SCPE)
!
NCOEF = IRBEF(NEF+1) - 1
CALL RHPSS (H, B, G, R, SCPH, DFH=DFH)
!
F = (SCPH(1,1)/DFH)/(SCPE(1,1)/DFE)
PVALUE = 1.0 - FDF(F,DFH,DFE)
CALL UMACH (2, NOUT)
WRITE (NOUT,*) 'Degrees of Sum of Prob. of'
WRITE (NOUT,*) ' Freedom Squares F-statistic Larger F'
WRITE (NOUT,99999) DFH, SCPH(1,1), F, PVALUE
99999 FORMAT (F8.1, 3X, 1F10.3, F11.3, 2X, F10.4)
END
Output
Degrees of Sum of Prob. of
Freedom Squares F-statistic Larger F
2.0 1178.135 2.746 0.0732
Example 2
The data for the second example are taken from Maindonald (1984, pages 203‑204). The data are saved in the matrix X. A multivariate regression model containing two dependent variables and three independent variables is fit using routine RGIVN. The sum of squares and crossproducts matrix is computed for the third independent variable in the model.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER INTCEP, LDB, LDG, LDH, LDR, LDSCPE, LDSCPH, LDX, &
NCOEF, NCOL, NDEP, NH, NIND, NROW, J, LDU
PARAMETER (INTCEP=1, LDU=1, NCOL=5, NDEP=2, NH=1, NIND=3, &
NROW=9, LDG=NH, LDH=NH, LDSCPE=NDEP, LDSCPH=NDEP, &
LDX=NROW, NCOEF=INTCEP+NIND, LDB=NCOEF, LDR=NCOEF)
!
INTEGER IDEP, IIND, INDDEP(1), INDIND(1),&
NOUT, NRMISS
REAL B(LDB,NDEP), D(NCOEF), DFE, DFH, G(LDG,NDEP), &
H(LDH,NCOEF), R(LDR,NCOEF), SCPE(LDSCPE,NDEP), &
SCPH(LDSCPH,NDEP), X(LDX,NCOL)
!
DATA (X(1,J),J=1,NCOL)/7.0, 5.0, 6.0, 7.0, 1.0/
DATA (X(2,J),J=1,NCOL)/2.0, -1.0, 6.0, -5.0, 4.0/
DATA (X(3,J),J=1,NCOL)/7.0, 3.0, 5.0, 6.0, 10.0/
DATA (X(4,J),J=1,NCOL)/-3.0, 1.0, 4.0, 5.0, 5.0/
DATA (X(5,J),J=1,NCOL)/2.0, -1.0, 0.0, 5.0, -2.0/
DATA (X(6,J),J=1,NCOL)/2.0, 1.0, 7.0, -2.0, 4.0/
DATA (X(7,J),J=1,NCOL)/-3.0, -1.0, 3.0, 0.0, -6.0/
DATA (X(8,J),J=1,NCOL)/2.0, 1.0, 1.0, 8.0, 2.0/
DATA (X(9,J),J=1,NCOL)/2.0, 1.0, 4.0, 3.0, 0.0/
DATA H/3*0.0, 1.0/, G/0.0, 0.0/
!
IIND = -NIND
IDEP = -NDEP
CALL RGIVN (X, IIND, INDIND, IDEP, INDDEP, B, R=R)
CALL RHPSS (H, B, G, R, SCPH, DFH=DFH)
CALL UMACH (2, NOUT)
WRITE (NOUT,*) 'DFH = ', DFH
CALL WRRRN ('SCPH', SCPH)
END
Output
DFH = 1.00000
SCPH
1 2
1 100.0 -40.0
2 -40.0 16.0
