Topics in Categorical Data Analysis

1. Introduction

• Up to this point, we have concentrated on statistical models that analyze continuous variables.
   
• Categorical Data Analysis considers count or crosstabulated data, rather than continuous measurements.
   
• This document discusses various tools and techniques that are useful for categorical data analysis.
   
• Here are some of the discrete probability distributions that are useful for categorical data analysis:
   
  
  Bernoulli (returns 1 with probability p, 0 otherwise)
   
  Binomial (number of successes in n trials)
   
  Multinomial (returns one of many outcomes with specified probabilities)
   
  Negative Binomial (number of trials to obtain a success)
   
  Discrete Uniform (self explanatory)
   
  Poisson (number arrivals in a specified time or region)
   
  Hypergeometric (number of successes in n trials when drawing without replacement)
  
  
• We are specifically interested in generalized regression models where the dependent variable is a value obtained from a Bernoulli random variable and the independent variables determine the parameter for this random variable. Such models are called logistic regression models.

2. Bernoulli and Binomial Distributions

• A Bernoulli random variable outputs the value 1 with probability p and 0 with probability 1 - p. Here is its discrete probability distribution:
   
  

  Outcome xi	Probability pi
  0	s
  1	1-s

  
  
• The probability of success (outcome = 1) is s; the probability of failure (outcome = 0) is 1-s.
   
• This probability distribution can also be written as
   
  
  p(x) = s^x (1 - s)^1-x, where x = 0 (failure) or 1 (success).
  
  
• The expected value of a Bernoulli random variable is s:
   
  
  E(x) = ∑_i=0^k x_i p_i = 0(1 - s) + 1s = s.
  
  
• The variance of a Bernoulli random variable is p(1-p):
   
  
  Var(x) = E(x - E(x))² = ∑_i=0^k(x_i - p)² p_i = (0 - s)²(1 - s) + (1 - s)²s
   
  = s(1 - s)(1 - s + s) = s(1 - s)
  
  
• What are some situations that can be modeled with a Bernoulli distribution?
   
• A Binomial random variable is the sum of n independent Bernoulli outcomes.
   
• Here is the probability distribution function for a binomial random variable where s is the probability of success:
   
  
  p(x) = s^x (1 - s)^1-x n!/[x! (n - x)!]
   
  
• The expected value of Binomial random variable is np; the variance of a Binomial random variable is np(1 - p).
   
• Use the R functions pbinom and qbinom to obtain probability information about binomial random variables.
   
• Use the R function rbinom to generate binomial random numbers.
   
• Practice Problems   Use R for the following:
   
  
  1. If a multiple choice test has 20 questions, what is the probability of passing (score of 14/20 or higher) the test by guessing?
      
  2. What is the probability of obtaining a score of exactly 10?
      
  3. Simulate the scores of 30 students taking the test if they are all guessing.

3. Poisson Distribution

• A Poisson process is a random variable that counts the number of random arrivals that occur in a fixed time interval. The arrivals are independent of each other and the probability that two arrivals occur simultaneously is 0. The probability that an arrival occurs in a time interval is proportional to the length of the time interval.
   
• Here are some random variables that can be modeled by Poisson processes:
   
  
  The number of sunspots in a given month
   
  Arrivals at a freeway toolbooth in one minute
   
  Telephone arrivals at a switchboard in one minute
   
  Lobsters caught in a lobster trap in one week
   
  Radioactive particles detected by a Geiger counter in one second
   
  Requests for a document on a Web server in one minute
  
  
• A Poisson distribution function shows the probability of x arrivals in a time period if the mean number of arrivals in that time period is λ. Here is the Poisson distribution function
   
  
  p(x) = e^-λ λ^x / x!
  
  
• Practice Problems: Use R for the following:
   
  
  1. A web page on a server is requested on the average 10 times per hour. What is the probability that this page is accessed exactly 10 times in an hour?
      
  2. Simulate the number times that the webpage is accessed per hour for 10 different hours.

4. Maximum Likelihood Estimators

• The likelihood function of a parameter θ by ∏_i=1ⁿ p(x_i | θ), where p(x | θ) is the probability distribution function for a discrete distribution or the density function for a continuous distribution.
   
• In the traditional frequentist framework for statistical models, the parameter θ is fixed and unknown and must be estimated using the data.
   
• By contrast, in the Bayesian framework, the parameter is a random variable and the likelihood function is equivalent to the posterior distribution of θ, given the data and assuming a uniform prior distribution for θ.
   
• The maximum likelihood estimator (MLE) for a parameter in a probability distribution is that value that maximizes the likelihood function or the log-likelihood function.
   
• To find the MLE, it is often more convenient to minimize the log-likelihood function defined by
   
  
  l = log ∏_i=1ⁿ p(x_i | θ) = ∑_i=1ⁿ log p(x_i | θ)
  
  (l is a script letter el.)
   
• Here is a simple example of a maximum likelihood estimator for a discrete uniform distribution. Consider a deck of cards numbered 1, 2, ... , N. Draw a card from the deck without showing how large the deck is. On the basis of this single card, what is the maximum likelihood estimator of N?
   
• Practice Problems: Find the MLE for the parameters in these situations.
   
  
  1. The parameter s in the Bernoulli distribution. The log-likelihood function is
      
     
     l(s) = ∑_i=1ⁿ log [s ^ x_i (1 - s) ^ (1 - x_i)] = ∑_i=1ⁿ [x_i log s + (1 - x_i) log(1 - s)]
      
     Ans: s = ∑_i=1ⁿ x_i / n = number of successes / number of trials.
     
     
  2. The parameter λ in the Poisson distribution. The log-likelihood function is
      
     
     l(λ) = ∑_i=1ⁿ log(e^-λ λ ^ x_i / x_i!) = ∑_i=1ⁿ (-λ + x_i log λ - log x_i!)
      
     Ans: λ = ∑_i=1ⁿ x_i / n = x
     
     
  3. The parameter μ in the normal distribution if σ is known to be 1. The log-likelihood function is
      
     
     l(μ) = ∑_i=1ⁿ log[(1/√2π) exp[(x_i - μ)²)/2] = ∑_i=1ⁿ log((1/√2π) + ∑_i=1ⁿ (x_i - μ)²/2
      
     Ans: μ = x.
     
     
  Hint for Problems 1, 2, and 3: compute the derivative (d/dx)l, set it to zero, and solve for the parameter. Recall that the derivative of log(x) is 1/x.
   
• Here is an example plot of the log-likelihood function for the Bernoulli distribution, when the number of trials is 5 and the number of successes is 2.

5. Generalized Regression

• The binomial and Poisson is useful for predicting successes and arrivals in restricted situations. However, in real world situations, the success rate or arrival rate depends on the situation.
   
• For example, it may be interesting to know that the average probability of contracting the AIDS disease is p = 0.005. However, the actual probability of contracting AIDS for a specific person depends on many factors, such as his or her gender, location of residence, income level, availability of public health programs, etc.
   
• It is useful to use a regression equation to predict the success rate or arrival rate. Such a regression equation is called a generalized regression model:
   
  
  g(θ) = β₀ + ∑_j=1ⁿ β_j x_ij
  
  The function g is called the link function.
   
• Here are the standard link functions used for three types of generalized linear regression:
   
  

  Regression Type	Distribution	Link Function	Generalized Linear Regression Equation
  Ordinary	Normal	Identity	μ = β0 + ∑j=1n βj xi
  Logistic	Bernoulli	Logit	log(p/(1-p)) = β0 + ∑j=1n βj xi
  Poisson	Poisson	Log	log(λ) = β0 + ∑j=1n βj xi

  
  
• Maximum likelihood estimation is used to obtain the parameter estimates for generalized regression.

6. Odds Ratios

• An alternative to specifying the probability that an event occurs is to specify the odds ratio for that event.
   
• Las Vegas may post the odds that a boxer's odds of winning a fight are 7 to 1 or 1/7. This means the probability of success (he wins the fight) over the probability of failure (he loses) is 1/7.
   
• Odds of 1 to 1 or 1/1 means that P(success)/P(failure) = 1, or P(success) = P(failure).
   
• If ODDS > 1, P(success) > P(failure); if ODDS < 1, P(success) < P(failure).
   
• To calculate the odds ratio of an event from its probability p,
   
  
  ODDS = P(success) / P(failure) = p / (1 - p).
  
  
• To calculate the probability of an event from its odds ratio, solve for p in the preceding equation:
   
  
  p = ODDS / (1 + ODDS)

7. Logistic Regression

• Logistic regression is useful if the dependent variable is binary (0 or 1) with continuous or dummy independent variables.
   
• A logistic regression equation predicts the probability that y = 1 for given values of the dependent variables.
   
• More specifically the dependent variables predict the log of the odds ratio from the dependent variables.
   
• The odds ratio (ODDS) is the the ratio of probability of success to the probability of failure.
   
• Let 1 represent success and 0 represent failure:
   
  
  ODDS = p / (1 - p) = P(y = 1) / P(y = 0)
  
  
• Thus the logistic regression equation
   
  
  log(p/(1-p)) = β₀ + ∑_j=1ⁿ β_j x_j
  
  can be written as
   
  
  log(ODDS) = β₀ + ∑_j=1ⁿ β_j x_j
  
  
• Example 1: A logistic regression model is set up to predict the probability that a person gets a cold in a given three month period based on whether he or she has taken a supplement, for which the active ingredient is ascorbic acid:
   
  
  log(ODDS) = -0.94 - 0.45 x,
  
  where x is a dummy variable such that x = 1 when the subject takes the supplement and x = 0 when the subject does not take it. For a person with no supplement, the probability of getting a cold is
   
  
  log(ODDS) = -0.94 - 0.45 * 0 = -0.94.
  
  This means that ODDS = exp(-0.94) = 0.391. Furthermore ODDS = p / (1 - p), so p = ODDS / (ODDS + 1) = 0.391 / (0.391 + 1) = 0.281: the probability that a subject gets a cold with no supplement is 0.281. When x = 1 (subject takes supplement),
   
  
  log(ODDS) = -0.94 - 0.45 * 1 = -1.39,
  
  ODDS = exp(-1.39) = 0.249, and p = ODDS / (ODDS + 1) = 0.249 / (1 + 0.249) = 0.199.
   
• Example 2: Suppose we have the independent variables x₁ = 34.32, and x₂ =16.41, and also the estimated regression parameters β₀ = -0.0211, β₁ = 0.0345, and β₂ = 0.0421, substitute into the logit equation:
   
  
  logit(p) = β₀ + β₀ x₁ + β₂ x₂ = -0.211 + 0.0345 * 34.32 + 0.0421 * 16.41 = 2.0859
  
  Then solve for p:
   
  

  logit(p)	=	2.0859
  log[p/(1-p)]	=	2.0859
  p/(1-p)	=	exp(2.0859)
  p/(1-p)	=	8.0518
  p	=	8.0518(1-p)
  9.0518p	=	8.0518
  p	=	0.8895

  
  
• Look at the Discrim, Depress and Remission Examples.
   
• To solve for the estimated regression parameters, use the method of maximum likelihood combined with the Newton-Raphson method for finding the actual maximum: find the parameters β₀, β₁, ... , β_k that maximize the log-likelihood function
   
  
  l = log p(y₁, ... , y_n | independent variables' values, β₀, β₁, ... , β_k)

8. The χ-² Distribution

• A χ-² random variable with k degrees of freedom is obtained as the sum of k independent N(0, 1) random variables.
   
• Try this R experiment to show that the chisquared random variable x with 4 degrees of defined by
   
  

  x <- rnorm(10000)^2 + rnorm(10000)^2 + 
       rnorm(10000)^2 + rnorm(10000)^2
  

  agrees with the R chisq functions.
   
• If p is a probability between 0 and 1, then the return values of the function calls quantile(x, p) and qchisq(p, 4) should be approximately equal.
   
• The random vector x could also be generated by rchisq(10000, 4).

9. Traditional χ² Tests

1. χ² Goodness-of-fit Test
    
   
   • Problem 1: Test whether a coin is fair. Flip a coin 50 times recording the number of heads and tails. Perform a χ² test for equal proportions to test whether the number of heads equals the number of tails. We simulate the coin flips using R:
     

     x = rbinom(50, 1, 0.5)
     cat("Number of heads =", sum(x), "out of 50 coin flips.\n")
     

     The output is
     

     Number of heads = 21 out of 50 coin flips.
     

     This result can be summarized in a contingency table:
      
     

     Cell	1	2
     Result	Head	Tail
     Count	21	29

     
     The five steps of the χ² goodness-of-fit test:
      
     
     1. State the null and alternative hypotheses:
         
        
        P(outcome₁) = P(outcome₂)
        
        
     2. Compute the Pearson χ² test statistic. O_i is the observed frequency for cell i; E_i is the expected frequency of cell i assuming the null hypothesis.
         
        
        X² = ∑_i=1^C (O_i - E_i)² / E_i = (21 - 25)² / 25 + (21 - 25)² / 25 = 1.28.
        
        
     3. Find a 95% confidence interval for the test statistic. X² has a χ² distribution with c - 1 degrees of freedom. Here c is the number of cells = 2. Looking in the χ² table, we find the confidence interval as [0, 5.024].
         
     4. Accept the null hypothesis because 1.28 ∈ [0, 5.024].
         
     5. Use SAS or R to calculate the p-value as 0.2579. See the FairCoin1 or FairCoin2 examples to see the scripts for computing p.
     
     Conclusion: accept the null hypotheses that the coin (or random number generator) is fair.
      
   • Look at the FairCoin1, FairCoin2, and FairCoin3 Examples.
   
   
2. χ² Test for Independence
    
   
   • Recall that if A and B are independent events, then the multiplication rule holds:
      
     
     P(A and B) = P(A) P(B)
     
     
   • Problem 2: Consider the following contingency table of data. For the values of Cold, 1 = Cold, 0 = No Cold.
      
     

     Observed Values		Cold
     		1	0	Totals
     Treatment	Placebo	31	109	140
     	Supplement	17	122	129
     	Totals	48	231	279

     
     We want to test whether the probability of contracting a cold is independent of whether or not the subject takes the supplement. Now the proportion of subjects that contract a cold P(Cold=1) = 48 / 279 = 0.1720  and the proportion of subjects that take the supplement is P(Supplement) = 140 / 279 = 0.4982.  If we assume that Treatment and Cold are independent, P(Placebo and Cold=1) = 0.1710 * 0.4982 = 0.0857, so the expected number of subjects that take the placebo is 0.0858 * 279 = 23.91.  This gives us this table:
      
     

     Expected Values Assuming Independence		Cold
     		1	0	Totals
     Treatment	Placebo	24.09	115.91	140.00
     	Supplement	23.91	115.09	129.00
     	Totals	48.00	231.00	279.00

     
     
     We want to test whether Cold and Treatment are independent. The χ² test for this is very similar to the goodness-of-fit test in the previous section. Here are the five steps of this test:
      
     
     1. H₀: The effects Cold and Treatment are independent.
        H₁: The effects Cold and Treatment are not independent.
         
     2. Compute the test statistic.  As with the goodness-of-fit test, O_ij represents an observed value and E_ij represents an expected value assuming independence.
         
        
        X² = ∑_i=1^C ∑_j=1^R (O_ij - E_ij)² / E_ij
            = (31 - 24.09)² / 24.09 + (109 - 115.91)² / 115.91 +
               (17 - 23.91)² / 23.91 + (122 - 115.09)² / 115.09 = 4.81.
        
        
     3. Although there are four entries in the contingency table, only one of these entries is free to vary freely; the other three entries are determined by the marginal totals. Therefore, there is only one degree of freedom for this test. Look up the level α=0.05 critical value for 1 degree of freedom in the χ² table: 3.84. This means that the acceptence region for the test is [0, 3.84].
         
     4. Reject H₀ because 4.81 ∉ [0, 3.84].
         
     5. Let SAS or R compute the p-value: 0.028. (We know that p < 0.05 because we are rejecting H₀).
     
     
   • The Cold1 Example conducts the χ² test for independence.
      
   • If the number of subjects is large, the independent two-sample z-test gives the same result as the χ² test for independence. See the Cold2 Example for details.
      
   • See the Cold3 Example, which uses a logistic regression model for predicting the probability of contracting a cold.

10. The Maximum Likelihood χ² Statistic

• An alternative to the Pearson χ² statistic is the likelihood ratio χ² statistic G²:
   
  
  G² = 2 O_i ∑_i=1^C log(O_i / E_i)
  
  
• This test statistic is derived as the log of the following ratio of likelihood functions:
   
  
  G² = -2 log[max_{θs statisfy H0} L(θ₁, ... , θ_k) / max_{ψs are unrestricted} L(ψ₁, ... , ψ_k)]
  
  
• See this paper by Jesse Hoey for motivation and details.
   
• Example:   Compute the G-statistic for the Cold data in Section 8.
   
  
  G² = 2 O_i ∑_i=1^C log(O_i / E_i)
     = 2[ 31 log(31 / 24.09) + 109 log(109 / 115.91) +
            17 log(17 / 23.91) + 122 log(122 / 115.09) ] = 4.87.
  
  Find this value in the SAS output for the Cold1 Example.
   
• G² is used instead of the F statistic for logistic regression.
   
• The biggest advantage of the G² statistic over the Pearson χ² statistic is that it works better for small sample sizes.

11. Tests of Significance for Logistic Regression

• Because a logistic regression model is nonlinear, the R-squared and F statistics cannot be used to assess the adequacy of the model.
   
• For testing whether an estimated parameter β_j^ is significant in a logistic regression model, use the Wald χ² Test, which uses the test statistic
   
  
  z = (β_j^ - β) / SE(β_j^),
  
  It has an approximately N(0, 1) distribution, assuming the null hypothesis: H⁰: β_j = 0. Therefore, z² has an approximate χ²(1) distribution. For a binomial random variable x, E(x) = np and Var(x) = np(1-p); E(x) and Var(x) both depend on the probability of success p.  They are not independent as they are for the normal distribution.  This is also true for logistic regression where, additionally, the value of p changes, depending on the independent variables.
   
• Here are three options for checking the global null hypothesis (same idea as the overall F-test for regression) that all regression coefficients are equal to zero (H₀: β₁ = ... = β_p-1 = 0).
   
  
  1. The Wald χ² Test can also be used for testing that all regression coefficients are equal to zero. Instead of normalizing a single parameter β_j with its standard error SE(β_j), the estimated parameter vector (β₁ ... β₁)^T is normalized by its standard error matrix. Details are omitted.
      
  2. The Likelihood Ratio χ² Test is also used for testing whether an estimated parameter β_j^ is significant. Use
      
     
     G² = -2 log(L₁ / L₀) = -2 (log L₁ - log L₀)
     
     where log L₁ is the maximum log-likelihood function with all the parameters being tested included in the model, and log L₀ is the log-likelihood function with all parameters except β₀ removed from the model. See Section 10 for more details about the likelihood ratio G² statistic.
      
  3. The Efficient Score Test, also called Rao's Score Test or the Lagrange Multiplier Test in econometrics. Details are omitted.
   
• The Wald, Maximum Likelihood, and Efficient Score tests can also be used for the Global log likelihood ratio statistic is computed as
   
  
  LLR = -2 log(L₁ / L₀) = -2 (log L₁ - log L₀)
  
  where L₁ is the probability function for the full model and L₀ is the probability function for the null hypothesis model with the parameters being tested removed from the model. This statistic also has an approximately chi-squared distribution.
   
• Look at the Discrim, Depress and Remission Examples.

12. Deviances for Logistic Regression

• The deviances for logistic regression are a scaled version of the residuals for standard regression.
   
• To calculate the deviances, obtain the model object from the glm function. Also, assume that the vector a contains the actual values of the dependent variables. Then:
  

  p <- fitted(model, type="response")
  deviances = sqrt(-2 * 
     log(abs(p - (1 - a))) * sign(a - p) 
  

• The deviances are related to the maximum likelihood test for the global null hypothesis.
   
• Parameter estimates for logistic regression consists of finding parameters such that the sum of squares of deviances are minimized.

13. Rank Correlation Statistics

• To assess the accuracy of a logistic regression model, proc logistic provides four rank correlation statistics: c, Somer's D, Goodman-Kruskal Gamma, Kendell's Tau-a.
   
• To define these statistics, look at all pairs of responses r_i, r_j that differ in value: r_i = 0 and r_j = 1 or r_i = 1 and r_j = 0.
   
• If an observation with response 1 has a higher predicted probability, the pair is called concordant.
   
• If an observation with response 1 has a higher predicted probability, the pair is called discordant.
   
• If an observation is neither concordant or discordant, it is called a tie.
   
• Let N be the number of observations, t the number of pairs with different responses, n_c the number of concordant pairs, n_d the number of discordant pairs, and t - n_c - n_d the number of ties.
   
• Then define
   
  
  c = (n_c + 0.5(t - n_c - n_d)) / t
   
  Somer's D = (n_c - n_d) / t
   
  Goodman Kruskal γ = (n_c - n_d) / (n_c + n_d)
   
  Kendall's τ-a = (n_c - n_d) / (2N(N - 1))
  
  
• If each of these measures is large, there is good agreement between the observed responses and their predicted probabilities.