g11sa:: Contingency Table Analysis (NAG Toolbox)

Given a set of

p

dichotomous variables

\tilde{x} = {(x_{1}, x_{2}, \dots, x_{p})}^{'}

, where

^{'}

denotes vector or matrix transpose, the objective is to investigate whether the association between them can be adequately explained by a latent variable model of the form (see Bartholomew (1980) and Bartholomew (1987))

G \{π_{i} (θ)\} = α_{i 0} + α_{i 1} θ .

(1)

The

x_{i}

are called item responses and take the value

0

1

θ

denotes the latent variable assumed to have a standard Normal distribution over a population of individuals to be tested on

p

items. Call

π_{i} (θ) = P (x_{i} = 1 ∣ θ)

the item response function: it represents the probability that an individual with latent ability

θ

will produce a positive response (1) to item

i

α_{i 0}

and

α_{i 1}

are item parameters which can assume any real values. The set of parameters,

α_{i 1}

, for

i = 1, 2, \dots, p

, being coefficients of the unobserved variable

θ

, can be interpreted as ‘factor loadings’.

G

is a function selected by you as either

Φ^{- 1}

or logit, mapping the interval

(0, 1)

onto the whole real line. Data from a random sample of

n

individuals takes the form of the matrices

X

and

R

defined below:

X_{s \times p} = [\begin{array}{l} x_{11}^{} & x_{12} & \dots & x_{1 p} \\ x_{21} & x_{22} & \dots & x_{2 p} \\ ⋮ & ⋮ & ⋮ \\ x_{s 1} & x_{s 2} & \dots & x_{s p} \end{array}] = [\begin{array}{l} {\tilde{x}}_{1} \\ {\tilde{x}}_{2} \\ ⋮ \\ {\tilde{x}}_{s} \end{array}], R_{s \times 1} = [\begin{array}{l} r_{1}^{} \\ r_{2} \\ ⋮ \\ r_{s} \end{array}]

where

{\tilde{x}}_{l} = (x_{l 1}, x_{l 2}, \dots, x_{l p})

denotes the

l

th score pattern in the sample,

r_{l}

the frequency with which

{\tilde{x}}_{l}

occurs and

s

the number of different score patterns observed. (Thus

\sum_{l = 1}^{s} r_{l} = n

). It can be shown that the log-likelihood function is proportional to

\sum_{l = 1}^{s} r_{l} \log P_{l},

where

P_{l} = P (\tilde{x} = {\tilde{x}}_{l}) = \int_{- \infty}^{\infty} P (\tilde{x} = {\tilde{x}}_{l} ∣ θ) ϕ (θ) d θ

(2)

(

ϕ (θ)

being the probability density function of a standard Normal random variable).

P_{l}

denotes the unconditional probability of observing score pattern

{\tilde{x}}_{l}

. The integral in (2) is approximated using Gauss–Hermite quadrature. If we take

G (z) = logit z = \log (\frac{z}{1 - z})

in (1) and reparameterise as follows,

\begin{array}{lcl} α_{i} & = & α_{i 1}, \\ π_{i} & = & {logit}^{- 1} α_{i 0}, \end{array}

then (1) reduces to the logit model (see Bartholomew (1980))

π_{i} (θ) = \frac{π_{i}}{π_{i} + (1 - π_{i}) \exp (- α_{i} θ)} .

If we take

G (z) = Φ^{- 1} (z)

(where

Φ

is the cumulative distribution function of a standard Normal random variable) and reparameterise as follows,

\begin{array}{lcl} α_{i} & = & \frac{α_{i 1}}{\sqrt{(1 + α_{i 1}^{2})}} \\ γ_{i} & = & \frac{- α_{i 0}}{\sqrt{(1 + α_{i 1}^{2})}} \end{array},

then (1) reduces to the probit model (see Bock and Aitkin (1981))

π_{i} (θ) = ϕ (\frac{α_{i} θ - γ_{i}}{\sqrt{(1 - α_{i}^{2})}}) .

An E-M algorithm (see Bock and Aitkin (1981)) is used to maximize the log-likelihood function. The number of quadrature points used is set initially to

10

and once convergence is attained increased to

20

The theta score of an individual responding in score pattern

{\tilde{x}}_{l}

is computed as the posterior mean, i.e.,

E (θ ∣ {\tilde{x}}_{l})

. For the logit model the component score

X_{l} = \sum_{j = 1}^{p} α_{j} x_{l j}

is also calculated. (Note that in calculating the theta scores and measures of goodness-of-fit nag_contab_binary (g11sa) automatically reverses the coding on item

j

α_{j} < 0

; it is assumed in the model that a response at the one level is showing a higher measure of latent ability than a response at the zero level.)

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

Accuracy

Further Comments

Timing

Initial Estimates

Heywood Cases

As in normal factor analysis, Heywood cases can often occur, particularly when

p

is small and

n

not very big. To overcome this difficulty the maximum likelihood search function is terminated when the absolute value of one of the

α_{j 1}

exceeds

10.0

. You have the option of deciding whether to exit from nag_contab_binary (g11sa) (by setting

ifail = 0

on entry) or to permit nag_contab_binary (g11sa) to proceed onwards as if it had exited normally from the maximum likelihood search function (setting

ifail = - 1

on entry). The elements in a, c, alpha and pigam may still be good approximations to the maximum likelihood estimates. You are advised to inspect the elements g to see whether this is confirmed.

Goodness of Fit Statistic

First and Second Order Margins

Example

function g11sa_example


fprintf('g11sa example results\n\n');

n = int64(1000);
x = [false, false, false, false;
     true,  false, false, false;
     false, false, false, true;
     false, true,  false, false;
     true,  false, false, true;
     true,  true,  false, false;
     false, true,  false, true;
     false, false, true,  false;
     true,  true,  false, true;
     true,  false, true,  false;
     false, false, true,  true;
     false, true,  true,  false;
     true,  false, true,  true;
     true,  true,  true,  false;
     false, true,  true,  true;
     true,  true,  true,  true];
irl = [int64(154);  11;  42;  49;
                  2;  10;  27;  84;
                 10;  25;  75; 129;
                 30;  50; 181; 121];

% Initial values
a = [0.5;     0.5;     0.5;     0.5];
c = [0;       0;       0;       0];

% Parameters
gprob  = false;
cgetol = 0.0001;
chisqr = true;
iprint = int64(-1);

% Fit a latent variable model
[x, irl, a, c, niter, alpha, pigam, cm, g, expp, obs, exf, y, ...
  iob, rlogl, chi, idf, siglev, ifail] = ...
  g11sa( ...
         n, gprob, x, irl, a, c, cgetol, chisqr, 'iprint', iprint);

% Display results
fprintf('Log likelihood kernel on exit = %14.4e\n\n', rlogl);

fprintf('Maximum likelihood estimates of item parameters are as follows\n');
fprintf('--------------------------------------------------------------\n\n');
fprintf('%8s%12s%11s%16s%10s%13s\n\n', 'item j', 'alpha(j)', 's.e.', ...
        'alpha(j,0)', 'pi(j)', 's.e.');
ivar = [1:numel(a)]';
se   = diag(cm);
results = [ivar alpha se(1:2:end) c pigam se(2:2:end)];
fprintf('%5d%13.3f%13.3f%13.3f%13.3f%13.3f\n',results');
fprintf('\n\n');

mtitle = 'Expected percentage of cases producing positive responses';
[ifail] = x04ca( ...
                 'Lower', 'Non-unit', expp, mtitle);
fprintf('\n');

mtitle = 'Observed percentage of cases producing positive responses';
[ifail] = x04ca( ...
                 'Lower', 'Non-unit', obs, mtitle);
fprintf('\n\n');

fprintf(' Observed   Expected   Theta   Component   Raw    Score\n');
fprintf(' frequency  frequency  score   score       score  pattern\n\n');
cs = double(x)*alpha;
results = [ double(irl)  exf  y  cs double(iob) double(x)];
fprintf('%7d%13.3f%8.3f%10.3f%8d    %1d%1d%1d%1d\n',results');
fprintf('---------  ---------\n');
fprintf('%7d%13.3f\n\n',n,n);

fprintf('Likelihood ratio goodness of fit statistic = %10.3f\n', chi);
fprintf('                        Significance level = %10.3f\n', siglev);
fprintf('(Based on %4d degrees of freedom)\n',idf)

g11sa example results

Log likelihood kernel on exit =    -2.4039e+03

Maximum likelihood estimates of item parameters are as follows
--------------------------------------------------------------

  item j    alpha(j)       s.e.      alpha(j,0)     pi(j)         s.e.

    1        1.045        0.148       -1.276        0.218        0.017
    2        1.409        0.179        0.424        0.604        0.022
    3        2.659        0.525        1.615        0.834        0.036
    4        1.122        0.140       -0.062        0.485        0.020


 Expected percentage of cases producing positive responses
             1          2          3          4
 1     25.8963
 2     19.0888    57.6547
 3     22.4987    47.9571    69.4214
 4     16.4381    33.8672    40.5658    48.7712

 Observed percentage of cases producing positive responses
             1          2          3          4
 1     25.9000
 2     19.1000    57.7000
 3     22.6000    48.1000    69.5000
 4     16.3000    33.9000    40.7000    48.8000


 Observed   Expected   Theta   Component   Raw    Score
 frequency  frequency  score   score       score  pattern

    154      147.061  -1.273     0.000       0    0000
     11       13.444  -0.873     1.045       1    1000
     42       42.420  -0.846     1.122       1    0001
     49       54.818  -0.747     1.409       1    0100
      2        5.886  -0.494     2.167       2    1001
     10        8.410  -0.399     2.455       2    1100
     27       27.511  -0.374     2.531       2    0101
     84       92.062  -0.332     2.659       1    0010
     10        6.237  -0.019     3.577       3    1101
     25       21.847   0.027     3.705       2    1010
     75       73.835   0.055     3.781       2    0011
    129      123.766   0.162     4.069       2    0110
     30       26.899   0.466     4.826       3    1011
     50       50.881   0.591     5.114       3    1110
    181      179.564   0.626     5.190       3    0111
    121      125.360   1.144     6.236       4    1111
---------  ---------
   1000     1000.000

Likelihood ratio goodness of fit statistic =      9.027
                        Significance level =      0.251
(Based on    7 degrees of freedom)

(i)	If the estimated factor loading for the $j$ th item is negative then that item is re-coded, i.e., $0$ s and $1$ s (or true and false) in the $j$ th column of x are interchanged.
(ii)	The rows of x will be reordered so that the theta scores corresponding to rows of x are in increasing order of magnitude.

$cm (2 \times i - 1, 2 \times i - 1)$	=	standard error $(alpha (i))$
$cm (2 \times i, 2 \times i)$	=	standard error $(pigam (i))$
$cm (2 \times i, 2 \times i - 1)$	=	correlation $(pigam (i), alpha (i))$ ,

$cm (2 \times i - 1, 2 \times j - 1)$	=	correlation $(alpha (i), alpha (j))$
$cm (2 \times i, 2 \times j)$	=	correlation $(pigam (i), pigam (j))$
$cm (2 \times i - 1, 2 \times j)$	=	correlation $(alpha (i), pigam (j))$
$cm (2 \times i, 2 \times j - 1)$	=	correlation $(alpha (j), pigam (i))$ ,

$idf = s_{0} - 2 \times p$	if $s_{0} < 2^{p}$ ;
$idf = s_{0} - 2 \times p - 1$	if $s_{0} = 2^{p}$ ,

On entry,	$ip < 3$ ,
or	$n < 7$ ,
or	$ns \leq 2 \times ip$ ,
or	$ns > n$ ,
or	$ns > 2^{ip}$ ,
or	two or more rows of x are identical,
or	$ldx < ns$ ,
or	$\sum_{l = 1}^{ns} irl (l) \neq n$ ,
or	at least one of $irl (l) < 0$ , for $l = 1, 2, \dots, ns$ ,
or	$maxit < 1$ ,
or	$ishow < 0$ ,
or	$ishow > 7$ ,
or	$ldcm < ip + ip$ ,
or	$ldexpp < ip$ ,
or	$lw < 4 \times ip \times (ip + 16)$ .

NAG Toolbox: nag_contab_binary (g11sa)

▸▿ Contents

Purpose

Syntax

Description