Categorical Data Analysis for Social Scientists: Labs

• Load the following file which contains information on housing tenure and

related details:

use http://teaching.sociology.ul.ie/categorical/labs/hten

• Examine the contents, and create a binary variable indicating home ownership, e.g.

recode tenure 1/2=1 3/99=0, gen(owner)

• The Wald test (\(\frac{\hat \beta}{SE} \sim \mathcal{N}(0,1)\) or \((\frac{\hat \beta}{SE})^2 \sim \chi^2\)) tests each parameter estimate individually against the null of no effect

. logit owner age

Iteration 0:   log likelihood = -8462.2668  
Iteration 1:   log likelihood = -8396.1105  
Iteration 2:   log likelihood = -8395.8593  
Iteration 3:   log likelihood = -8395.8593  

Logistic regression                               Number of obs   =      15499
                                                  LR chi2(1)      =     132.81
                                                  Prob > chi2     =     0.0000
Log likelihood = -8395.8593                       Pseudo R2       =     0.0078

------------------------------------------------------------------------------
       owner |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   .0119102   .0010456    11.39   0.000     .0098608    .0139596
       _cons |   .6421907   .0495884    12.95   0.000     .5449992    .7393822
------------------------------------------------------------------------------

• Here the PE for age is clearly significant

• The likelihood ratio test comparing the model with versus without that parameter tests the same hypothesis and is likely more stable with small data sets:

logit owner // The null model with no explanatory variables
estimates store null 
logit owner age
lrtest age

• Note this code will not work, for a good reason: age has a missing value so the data set changes between the two models
• Try this approach instead

logit owner age 
estimates store model1
logit owner if e(sample)
lrtest model1

• The e(sample) is temporary information created by the previous logistic regression, indicating which cases were use to fit it. This strategy means the two models are fitted on the same data.

• The advantage for the LR test for single variables is greater in small data sets, but it is very important where we add groups of variables, either blocks that go together for theoretical reasons, or as groups of dummies indicating categorical variables

logit owner age i.sex i.mastat
est store model2
logit owner age i.sex if e(sample)
...
lrtest model2

Likelihood-ratio test                                  LR chi2(6)  =    928.68
(Assumption: . nested in model2)                       Prob > chi2 =    0.0000

• The LR \(\chi^2\) attributable to adding marital status to the model is 928.68. Since mastat has 7 values, this requires 6 parameter estimates, hence 6 DF – this is a very strong effect

• Armed with these tools, search for a persuasive model to predict house-ownership

• For individual data, we can't compare predicted and observed values in quite the same way. estat gof generates the Pearson goodness of fit statistic, and estat gof, group(10) the Hosmer-Lemeshow test: both are analagous to the LR test against the saturated model for grouped data, but are far from perfect.
• Use the example1 data set as you have changed and saved it, or the tenure data to experiment with this:

gen univ = educ==1
logit univ age
estat gof, group(10)

• Long and Freese's book comes with a suite of add-ons called SPost

ssc install spost

• This includes the fitstat command which outputs a large range of fit-related statistics: examine these for a few of the models you have fitted.

• Also explore the classification table: estat class for a few models
• Examine the ROC curves for a few models using lroc

• Consider a linear regression with a non-linear RHS:

use http://teaching.sociology.ul.ie/categorical/labs/example1
keep if rdoby>0
gen age = 2008 -rdoby
reg rfimn c.age##c.age
predict pi 
scatter pi age, msize(0.2)

• The fitted effect of age is nonlinear, rising to a peak in the late 40s, then falling
• The marginal effect of age is (by differentiation) \(\beta_1 + 2\beta_2\times x\)
• Find the marginal effect of age at its mean: this is the slope of the curve at the mean value of age. Is it a useful summary?
• Now find the average marginal effect of age at each observed value of age and summarise it

gen marg = _b[age] + 2*age*_b[c.age#c.age]
su marg

• Compare your results with the official margins command:

margins, dydx(age)

• Logistic regression is non-linear in the probability
• The marginal effect is the slope of the tangent: \(\beta p (1-p)\)

use http://teaching.sociology.ul.ie/categorical/labs/shuttledata, clear
logit fail temp
predict p
gen slope = _b[temp]*p*(1-p)
su slope
margins, dydx(temp)

• Things are a little more complex with multiple variables

use http://teaching.sociology.ul.ie/categorical/labs/hten
recode tenure 1/2=1 3/99=0, gen(owner)
logit owner i.sex age
predict p //

Load the data at http://teaching.sociology.ul.ie/so5032/bhpsqual.dta and examine the variables. The variable vote has four categories. Examine bivariate relationships between vote and some of the other variables, and then search for a multinomial logistic regession that makes sense:

mlogit vote ...

• It is better to have a base category that is easily interpretable, and the default here will use the biggest category ("Other/nationalist"). To avoid this, use the baseoutcome(1) option, which will force category one as the base.

• Use the likelihood-ratio test for each variable, since there are three times as many parameters as usual

• Use predict to generate predicted values (note you have to supply one variable name to hold the prediction for each category). How often is the most probable predicted category the same as the observed one?

predict v1-v4

• Fit a multinomial logit with qual as the dependent variable

mlogit qual c.age##c.age i.sex pahgs

• Note any patterns in the predictors: do you think the relationship between them and qual suggests ordinality

• Fit a stereotype logit:

slogit qual c.age##c.age i.sex pahgs

• This constrains the mlogit as if the dependent variable could be placed on a single continuum (ideally but not necessarily in its "natural" order)
• How well does it work?

• We can also try ordered logistic regression:

ologit qual c.age##c.age i.sex pahgs

• How do the results compare?
• Are you happy with the assumption of proportional odds?
• We can test this with the Brant test (part of spostado). However, the brant command is a little old and doesn't work with the latest syntax. We need an older format of the model specification:

gen age2 = age*age
xi: ologit qual age age2 i.sex pahgs
brant

• If the Brant test is significant, try generalised ordinal logit

ssc install gologit2
xi: gologit2 qual age age2 i.sex pahgs

• Compare the output with mlogit : it is equally imparsimonious but is structured differently

• seqlogit is also available from SSC:

ssc install seqlogit
seqlogit qual c.age##c.age i.sex pahgs, tree(1 : 2 3 4 , 2 : 3 4 , 3 : 4 )

• You can constrain effects to be equal across contrasts:

constraint 1 [_3_4v2=_4v3]: 2.sex
seqlogit qual c.age##c.age i.sex pahgs, ///
 tree(1 : 2 3 4 , 2 : 3 4 , 3 : 4 ) ///
 constraint (1)

• You can test whether a constraint is allowable using the lrtest command

• Fit each of the models again, and run fitstat after each one, noting the BIC statistics
• Which gives the highest value of BIC?