Computes model estimates and covariances in a fitted log-linear model.
NCLVAL Vector of length NCLVAR containing, in its i-th element, the number of levels or categories of the i-th classification variable. (Input)
NVEF Vector of length NEF containing the number of classification variables associated with each effect. (Input)
INDEF Vector of
length NVEF(1) +
+ NVEF(NEF) containing, in
consecutive positions, the indices of the variables that are included in each
effect. (Input)
The entries in INDEF are sequenced so
that the first NVEF(1) elements
contain the indices of the variables in effect 1, the next NVEF(2) elements of
INDEF contain
the indices of the variables in effect 2, etc. See Comment 4 for an
example.
FIT Vector of
length NCLVAL(1)
* NCLVAL(2) *
* NCLVAL(NCLVAR) containing the
model estimates of the cell counts. (Input)
See Comment 3 for the
ordering of the elements of FIT.
To obtain a first iteration approximation to the optimal parameter values, the
observed counts may be input in FIT,
in which case a least-squares model is fit. In all cases, values of zero in
FIT
are assumed to correspond to structural zeros in the table. See the Description
section for details.
NCOEF Number of regression coefficients in the model. (Output)
COEF NCOEF by 4 matrix containing the estimated coefficients and associated statistics. (Output)
Col. Statistic
1 Coefficient estimate
2 Estimated standard error of the estimated coefficient
3 Asymptotic normal score for testing that the coefficient is zero
4 p-value associated with the normal score in column 3 (two-sided alternative)
COV NCOEF by NCOEF covariance matrix of the estimated coefficients. (Output)
NCLVAR Number
of classification variables. (Input)
A variable specifying a
margin in the table is a classification variable. The first classification
variable is named A, the
second classification variable is named B, etc.
Default: NCLVAR = size (NCLVAL,1).
NEF Number of
effects in the model. (Input)
A marginal table is implied by
each effect in the model. Lower-order effects should not be included since their
inclusion is automatic in the hierarchical models fit here (e.g., do not include
effects A or B if
effect AB is in the model).
Default: NEF = size (NVEF,1).
TOL Tolerance
used in determining linear dependence in COV.
(Input)
TOL
= 100.0 * AMACH(4) is a common
choice. See the documentation for routine AMACH
in Reference
Material.
Default: TOL = 1.19e-5 for
single precision and 2.d 14 for double precision.
IPRINT Printing
option. (Input)
Default: IPRINT = 0.
IPRINT |
Action |
0 |
No printing is performed. |
1 |
Printing of COEF and COV is performed. |
2 |
COEF, COV, and FIT are printed. |
In the printing, A * B(2) denotes the second variable in the AB interaction effect.
LDCOEF Leading
dimension of COEF exactly as
specified in the dimension statement in the calling program.
(Input)
Default: LDCOEF = size (COEF,1).
LDCOV Leading
dimension of COV
exactly as specified in the dimension statement in the calling
program. (Input)
Default: LDCOV = size (COV,1).
Generic: CALL CTPAR (NCLVAL, NVEF, INDEF, FIT, NCOEF, COEF, COV [, ])
Specific: The specific interface names are S_CTPAR and D_CTPAR.
Single: CALL CTPAR (NCLVAR, NCLVAL, NEF, NVEF, INDEF, FIT, TOL, IPRINT, NCOEF, COEF, LDCOEF, COV, LDCOV)
Double: The double precision name is DCTPAR.
Routine CTPAR computes estimates of parameters and associated variances and covariances in hierarchical loglinear models. A weighted least-squares algorithm is used.
A hierarchical analysis of variance model is a factorial analysis of variance model in which a lower-order effect is included in a model whenever a higher-order effect containing it is in the model. Thus, if the effect ADF is in the model, then effects A, D, F, AD, AF, and DF are automatically in the model.
Input to CTPAR
may be either the expected table values for the given hierarchical model as
output, for example, by routine PRPFT, or the observed table values. When the
fitted values are input, the estimates computed are the maximum likelihood
estimates. When observed values are input, weighted least-squares estimates of
the parameters in the log-linear model are computed.
(Least-squares
estimates and maximum likelihood estimates can also be computed via routines
CTWLS and CTGLM, respectively.)
When an expected count (as input in FIT) is zero, the cell is taken to be a structural zero. Such cells are not included in the weighted least-squares analysis. Estimates corresponding to structural zeros are set to the missing value indicator (NaN). To avoid this (and to determine the total degrees of freedom for each effect), add a positive constant such as 0.5 to each of the observed cell counts of zero, the sampling zeros. When structural zeros are present in the data the estimates may be written as
where
are vectors, and ρ → ∞ Routine CTPAR
estimates the finite portion of the estimate, The infinite portion,
ensures that the fitted values for cells corresponding to
structural zeros are zero (sampling zeros are considered to be structural zeros
in CTPAR).
If there are no structural zeros
Let fi denote the i-th element of the vector FIT. The asymptotic variance-covariance matrix of the cell counts is estimated by a diagonal matrix S = diag(f) where diag(f) denotes the diagonal matrix in which sij = 0 for i ≠ j and sii = fi along the diagonal. If X denotes the design matrix for the hierarchical model (with rows in X corresponding to structural zeros omitted), and yi = logfi, then the weighted least-squares estimates are
and the estimated variance-covariance matrix is
(XT S−1X)−1
(see Grizzle, Starmer, and Koch [1969]).
If main effect A has, for example, four levels, then the design matrix X contains three dummy variables corresponding to this effect. Main effect dummy variables are generated as follows: For an observation fi corresponding to level j of the effect, if j < 3, then the j-th dummy variable is set to 1 with the remaining dummy variables set to 0. If j = 4, then all three dummy variables are set to −1. Dummy variables for interactions are generated as the product of the corresponding dummy variables in the usual manner with the smallest index in the specification of the interaction varying fastest. The indices of the classification variables for each effect are always sorted from smallest to largest when computing the columns of X.
1. Workspace may be explicitly provided, if desired, by use of C2PAR/DC2PAR. The reference is:
CALL C2PAR (NCLVAR, NCLVAL, NEF, NVEF, INDEF, FIT, TOL, IPRINT, NCOEF, COEF, LDCOEF, COV, LDCOV, IRANK, NCVEF, IXEF, IINDEF, IA, INDCL, CLVAL, REG, X, D, XMIN, XMAX, WK)
The additional arguments are as follows:
IRANK Rank of COV.
NCVEF Vector of length 2NCLVARNCLVAR − 1.
IXEF Vector of length NCLVAR * 2NCLVARNCLVAR−1.
IINDEF Vector of length NVEF(1) + L + NVEF(NEF).
IA Vector of length NCLVAR.
INDCL Vector of length NCLVAR.
CLVAL Vector of length NCLVAL(1) + L + NCLVAL(NCLVAR).
REG Vector of length NCOEF + 1.
X Vector of length NCLVAR.
D Vector of length NCOEF.
XMIN Vector of length NCOEF.
XMAX Vector of length NCOEF.
WK Vector of length NCOEF + 1 if IPRINT ≠ 2. Otherwise, its length is the maximum of NCOEF + 1 and the product of the two largest elements of NCLVAL.
2. Informational errors
Type Code
3 5 The label for one or more of the tables exceeds the buffer limit.
3 11 The label for one or more effects exceeds the buffer limit.
4 1 LDCOEF or LDCOV is less than NCOEF.
3. The cells
of the vector FIT are sequenced so
that the first variable cycles from 1 to NCLVAL(1) the slowest,
the second variable cycles from 1 to NCLVAL(2) the next
slowest, etc., up to the NCLVAR-th variable,
which cycles from 1 to NCLVAL(NCLVAR) the fastest.
Example: For NCLVAR = 3, NCLVAL(1) = 2, NCLVAL(2) = 3, and
NCLVAL(3) = 2,
the cells of table X(I, J, K)
are entered into FIT(1) through FIT(12) in the
following order: X(1, 1, 1), X(1, 1, 2), X(1, 2, 1), X(1, 2, 2), X(1, 3, 1), X(1, 3, 2), X(2, 1, 1), X(2, 1, 2), X(2, 2, 1), X(2, 2, 2),
X(2, 3, 1),
X(2, 3, 2).
4. INDEF is used to describe the marginal tables to be fit. For example, if NCLVAR = 3 and the first effect is to fit the marginal table for variables 1 and 3 and the second effect is to fit the marginal table for variable 2, then: NEF = 2, NVEF(1) = 2, and NVEF(2) = 1. Since the sum of the NVEF(I) is 3, then INDEF is a vector of length 3 with values: INDEF(1) = 1, INDEF(2) = 3, and INDEF(3) = 2.
The example illustrates the use of CTPAR in a simple four-way table in which the first three factors have two levels, and the fourth factor has three levels. The data, which is taken from Lee (1977), involve the brand preference in different situations.
USE PRPFT_INT
USE CTPAR_INT
IMPLICIT NONE
INTEGER IPRINT, LDCOEF, LDCOV, LTAB, NCLVAR
PARAMETER (IPRINT=2, LDCOEF=13, LDCOV=13, LTAB=24, NCLVAR=4)
!
INTEGER INDEF(6), NCLVAL(NCLVAR), NCOEF, NVEF(3)
REAL COEF(LDCOEF,4), COV(LDCOV,LDCOV), FIT(LTAB), &
TABLE(LTAB)
!
DATA TABLE/19, 57, 29, 63, 29, 49, 27, 53, 23, 47, 33, 66, 47, &
55, 23, 50, 24, 37, 42, 68, 43, 52, 30, 42/
DATA NVEF/2, 2, 2/, INDEF/2, 4, 1, 4, 2, 3/
DATA NCLVAL/3, 2, 2, 2/, FIT/24*1.0/
!
CALL PRPFT (NCLVAL, TABLE, NVEF, INDEF, FIT)
!
CALL CTPAR (NCLVAL, NVEF, INDEF, FIT, NCOEF, COEF, COV, IPRINT=IPRINT)
!
END
Variable Number of Levels
1
A
3
2
B
2
3
C
2
4
D
2
----------
Table 1: B = 1 C = 1
D (row) by
A (column)
1 2
3
1 19.52 23.65
26.09
2 47.85 46.99
42.89
----------
Table 2: B = 1 C = 2
D (row) by
A (column)
1 2
3
1 28.39 34.40
37.94
2 69.58 68.32
62.37
----------
Table 3: B = 2 C = 1
D (row) by
A (column)
1 2
3
1 30.85 37.37
41.23
2 57.52 56.48
51.56
----------
Table 4: B = 2 C = 2
D (row) by
A (column)
1 2
3
1 25.24 30.58
33.73
2 47.06 46.21
42.18
Coefficient Statistics
Standard Asymptotic
Coefficient Error Z-statistic P-value
1 intercept 3.6827 0.0333 110.66 0.0000
2 A(1) -0.0591 0.0475 -1.24 0.2341
3 A(2) 0.0278 0.0461 0.60 0.5562
4 B -0.0166 0.0331 -0.50 0.6242
5 C -0.0434 0.0319 -1.36 0.1943
6 D -0.2783 0.0329 -8.45 0.0000
7 A*D(1) -0.1016 0.0475 -2.14 0.0506
8 A*D(2) 0.0034 0.0461 0.07 0.9414
9 B*C -0.1438 0.0319 -4.51 0.0005
10 B*D -0.0684 0.0328 -2.09 0.0558
Asymptotic Coefficient Covariance
1 2 3 4 5
1 1.1076E-03 9.7132E-05 -3.5887E-05 4.3244E-05 4.3786E-05
2 2.2562E-03 -1.1408E-03 -3.4043E-11 2.6829E-11
3 2.1232E-03 2.5675E-11 -5.1643E-11
4 1.0968E-03 1.4480E-04
5 1.0146E-03
6 7 8 9 10
1 2.9815E-04 1.3065E-04 -1.6147E-05 1.4480E-04 7.6307E-05
2 1.3065E-04 7.2117E-04 -4.0976E-04 6.2343E-11 -1.0681E-11
3 -1.6147E-05 -4.0976E-04 5.7437E-04 -4.9217E-11 -2.3482E-11
4 7.6307E-05 1.2601E-11 -4.1730E-11 4.3786E-05 2.8917E-04
5 -1.4272E-11 -5.5301E-11 4.2801E-11 4.5231E-06 -4.6962E-11
6 1.0851E-03 9.7132E-05 -3.5887E-05 -4.9749E-11 3.0847E-05
7 2.2562E-03 -1.1408E-03 5.9300E-11 -1.0361E-10
8 2.1232E-03 -2.4481E-11 2.9160E-11
9 1.0146E-03 1.1201E-11
10 1.0743E-03
PHONE: 713.784.3131 FAX:713.781.9260 |