# Can multinomial regression do what I want?

I have a list of observations of females and their responses, categorized into 5 behaviors, to a potential threat. I'm wondering if, for each response type, whether the presence of an infant makes it more or less likely a female would perform that response. (E.g., we might hypothesize that females with infants are more likely to hide and less likely to go looking for food.)

The following was suggested:

1. Bin females by social group (i.e., each bin would only contain females who associated with each other).
2. Calculate the proportion of each response type observed in each group (all these, for each group, would add to 1).
3. Arcsine-transform the proportion data.
4. Perform a t-test, for each behavior, to see whether the proportions of females responding that way differs depending on the presence of an infant.

This sounded sketchy to me, and it seems that others agree. The consensus seems to be that, in such cases, it's better to use multinomial regression (in this case, BEHAVIOR ~ INFANT_STATUS) to determine whether infant presence has an effect. However, I was wondering whether I could use those results determine whether (and in what direction) the presence of an infant affects the probability of each response behavior. Also, would that be possible if I were to include additional independent variables?

Any advice you can give, on analytical design or actual coding, is much appreciated. I've been doing this in R, but I'm more comfortable in MATLAB if that works as well. Thanks in advance.

The description suggests an analysis with chi-square test of proportions would be appropriate. You can summarize your data with a contingency table of the following form

$$\begin{array}{cccccc} & \mbox{Behaviour 1} & \mbox{Behaviour 2} & \mbox{Behaviour 3} & \mbox{Behaviour 4} & \mbox{Behaviour 5} \\ \mbox{Child} & N_{11} & N_{12} & N_{13} & N_{14} & N_{15} \\ \mbox{No child} & N_{21} & N_{22} & N_{23} & N_{24} & N_{25} \\ \end{array}$$

You would of course look at marginal frequencies to determine where cell specific differences lie. You can alternately calculate proportions as you suggested in point 2, but you shouldn't base inference on those values.

It should be noted that, by virtue of having a binary predictor, you can use logistic regression to obtain an equivalent $\chi_2$ test of proportions using the following R code:

glm(child ~ factor(behavior), family=binomial)


Logistic regression allows us to treat predictors as outcomes. The reason why this works is that odds ratios are conditioning invariant:

$$\frac{\mbox{Odds(Child | Behaviour)}}{\mbox{Odds(No Child | Behaviour)}} = \frac{\mbox{Odds(Behavior | Child)}}{\mbox{Odds(No behavior | Child)}}$$

1. If you need to bin females by social group, then two options are multivariate logistic regression or Cochran Mantel Haenszel. Example of the first in R is glm(child ~ factor(behavior) + factor(group), family=binomial). Example of the latter in R is mantelhaen.test(table(behavior, child, group)). The choice of whether to do this should rely solely upon the scientific question.

2. Aggregating the behavior groups into proportions will not allow you to perform any statistical inference, but would be useful in an inspection of trends. You can determine which proportions are higher.

3-4 are suggestions for applying tests for continuous data in circumstances where tests for categorical data are available, established, well known, and superior.

• I initially thought chi-squared test on the counts (actually, Fisher's exact test since one count is 0). The result is a significant association between Behavior and Child. However, I don't know how to go from that result to saying that females with vs. without children differ in any of the individual behaviors. From your example, let's say that I can see that N11 > N21 when converted to proportion (where each row sum to 1). How can I say that's a significant difference when there are other proportion pairs N1* and N2* that could have accounted for the association instead? – Sam R Feb 1 '13 at 18:00
• I like the idea of using logistic regression; that's out-of-the-box but makes sense. My only question is how to interpret the output. summary(glm(Child~Behaviour)) will give p-values for each of 4 behaviors compared to the 5th (reference) behavior, as well as for the intercept. Does the p-value of the intercept represent whether there's a significant difference between females with vs. without children for the reference behavior? If so, should I just run a series of glm models, one with each behavior as reference? Also — thanks for your help :-). – Sam R Feb 1 '13 at 18:07
• If the tabular (Fisher's Exact) test is significant, it's important to report and interpret the p-value from that analysis. Subsequently, you can make offhand observations about which proportions differ, but remember this is NHST and you haven't proven any alternative hypothesis. You can perform a global hypothesis test of factor behavior using the lmtest package: lrtest(glm(child ~ factor(behavior), family=binomial)) Gives you the likelihood ratio test for the behavior variable compared to the null (only proportions by child). – AdamO Feb 1 '13 at 18:31

Another option, one that I like better: Use logistic regression to predict the probability of each behavior based on a number of independent variables, including whether the female is a mother. A significant result for the latter (p-value of coefficient < 0.05 after correction) means that motherhood affects the probability of performing the given behavior, all other things being equal.

In order to actually test whether mothers vs. non-mothers are differently likely to perform each behavior, my advisor suggested that I perform a likelihood ratio test for each behavior. Here are the log-likelihoods, keeping in mind the following conventions:

• N = number (by itself, total sample size; with one of the following subscripts, the number falling into the designated category)
• m = mothers
• n = non-mothers
• p = performed behavior
• d = did not perform behavior

$$L_{null} = N_{p}ln(\frac{N_p}{N}) + N_{d}ln(\frac{N_d}{N})$$ $$L_{alt} = N_{m,p}ln(\frac{N_{m,p}}{N_m}) + N_{m,d}ln(\frac{N_{m,d}}{N_m}) + N_{n,p}ln(\frac{N_{n,p}}{N_n}) + N_{n,d}ln(\frac{N_{n,d}}{N_n})$$

Then the dispersion parameter $D = −2*(L_{null}-L_{alt})$ is $\chi^2$-distributed. We then find the one-tailed $\chi^2$ p-value. Of course, we use some sort of correction (I like the Sidak) to adjust the critical p-value due to multiple tests.

Does anybody see anything wrong with this approach?