Solved – predict() – multinomial logistic regression

multinomial-distributionr

I come to you today because I face a huge problem that I cannot explain.

I have run a multinomial logistic regression (using the mlogit package) on behavioral data. I prepare the data by doing

    mlogit <- mlogit.data(Merge, choice = "Choice", shape = "long", alt.var = "Comp", 
                          drop.index = TRUE)

on my Merge data.

which gives me the following:

                Date     Time ActivityX ActivityY Temp Behavior Valley Age Month Year kid Individual Choice
    1.F   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26   TRUE
    1.R   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE
    1.M   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE
    1.RUN 01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE
    2.F   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26   TRUE
    2.R   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE
    2.M   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE
    2.RUN 01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE
    3.F   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE
    3.R   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE
    3.M   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26   TRUE
    3.RUN 01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE
    4.F   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE
    4.R   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26   TRUE
    4.M   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE
    4.RUN 01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE
    5.F   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE
    5.R   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26   TRUE
    5.M   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE
    5.RUN 01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE

then I ran my regression :

m1 <- mlogit(Choice ~ 1 |Temp + Valley + Age + kid + Month , mlogit)

and it gave me significant results :

                          Estimate  Std. Error  t-value  Pr(>|t|)    
    M:(intercept)      -4.2153e-01  5.7533e-02  -7.3268 2.358e-13 ***
    R:(intercept)       6.2325e-01  3.4958e-02  17.8284 < 2.2e-16 ***
    RUN:(intercept)    -1.2275e+01  4.0526e-01 -30.2895 < 2.2e-16 ***
    M:Temp              1.5371e-02  9.8680e-04  15.5764 < 2.2e-16 ***
    R:Temp             -3.9871e-02  6.7926e-04 -58.6975 < 2.2e-16 ***
    RUN:Temp           -4.4532e-02  6.8696e-03  -6.4825 9.023e-11 ***
    M:ValleyTrupchun   -3.6154e-01  1.6362e-02 -22.0968 < 2.2e-16 ***
    R:ValleyTrupchun   -4.0186e-02  9.7968e-03  -4.1020 4.096e-05 ***
    RUN:ValleyTrupchun  1.2895e+00  8.5357e-02  15.1066 < 2.2e-16 ***
    M:Age              -1.1026e-02  2.6902e-03  -4.0985 4.158e-05 ***
    R:Age               1.9465e-02  1.6479e-03  11.8119 < 2.2e-16 ***
    RUN:Age             5.5473e-02  1.6661e-02   3.3294 0.0008703 ***
    M:kidY              6.0686e-02  2.2638e-02   2.6807 0.0073460 ** 
    R:kidY             -4.1638e-01  1.2391e-02 -33.6024 < 2.2e-16 ***
    RUN:kidY            6.2311e-01  1.0410e-01   5.9854 2.158e-09 ***
    M:Month            -2.0466e-01  8.4448e-03 -24.2346 < 2.2e-16 ***
    R:Month             2.4148e-02  5.2317e-03   4.6157 3.917e-06 ***
    RUN:Month           9.8715e-01  5.6209e-02  17.5622 < 2.2e-16 ***

those results were in line with what I expected to find in literature so I was quite happy.

My next step was to plot my results and here is when I have some trouble.

First of all when I plot my original data and compare it with the result of my regression I find some huge differences. For example, when I plot the %of time spend in a behavior (M for moving, F for feeding, R for resting and Run for running, in my regression F is the reference) in function of age, I find that the older an individual gets, the more they will rest and the more they will move, but the estimates I got from my regression shows that they should rest more (when they get older) but move less. So to summarize, my graph on the original data shows the opposite as what I got from the regression.

I don't know if it is normal, in the sense that I don't know if I can compare my original data to the result of my regression in a way that my regression shows the probability from switching to a behavior from an other each time my variable grows of one unit.

So I wanted to use the predict() function but I don't know how to do that. I was hoping to get some help here.

Best Answer

When the estimates do not have the expected sign, the usual suspect is multicollinearity. Below is a passage from the Wikipedia page.

The usual interpretation of a regression coefficient is that it provides an estimate of the effect of a one unit change in an independent variable, $X_1$, holding the other variables constant. If is highly correlated with another independent variable, $X_2$, in the given data set, then we have a set of observations for which $X_1$ and $X_2$ have a particular linear stochastic relationship. We don't have a set of observations for which all changes in $X_1$ are independent of changes in $X_2$, so we have an imprecise estimate of the effect of independent changes in $X_1$.

You can find more about the adverse effects of multicollinearity, and strategies to fight it by reading this question and that question.

Related Solutions

Logistic – Why Exact Singularity When Adding Dummy Variable to Multinomial Logistic Using R mlogit?

I think I have discovered the source of the error. In my model, as well as in the example above, the dummy is an individual specific attribute, but I had included it in the model statement as an alternative specific attribute.

Properly specifying the dummy as an individual-specific attribute (i.e. after the '|') produces a sensible result:

summary(mlogit(choice ~ pf + cl + loc + wk + tod + seas | odd.dummy-1 , data=Electr)) # Dummy as individual-specific attribute

Coefficients :
             Estimate Std. Error  t-value Pr(>|t|)    
pf          -0.624054   0.023369 -26.7048   <2e-16 ***
cl          -0.108130   0.008281 -13.0577   <2e-16 ***
loc          1.442730   0.050706  28.4527   <2e-16 ***
wk           0.997726   0.044911  22.2157   <2e-16 ***
tod         -5.454041   0.185186 -29.4517   <2e-16 ***
seas        -5.831122   0.188030 -31.0116   <2e-16 ***
odd.dummy:2  0.057435   0.067904   0.8458   0.3976    
odd.dummy:3  0.064776   0.068510   0.9455   0.3444    
odd.dummy:4  0.061038   0.067753   0.9009   0.3676

I'll leave it to someone else to try to explain the mathematics of what was happening.

Linear Regression in R – How to Interpret Model Diagnostics

This is a long and rambling question, so you are getting a long and rambling answer. Apologies. Using the example from the ?lm() call,

ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2,10,20, labels=c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
summary(lm.D9)
#output#
Call:
lm(formula = weight ~ group)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0710 -0.4938  0.0685  0.2462  1.3690 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.0320     0.2202  22.850 9.55e-15 ***
groupTrt     -0.3710     0.3114  -1.191    0.249    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6964 on 18 degrees of freedom
Multiple R-squared: 0.07308,    Adjusted R-squared: 0.02158 
F-statistic: 1.419 on 1 and 18 DF,  p-value: 0.249

I don't entirely understand your confusion on the "coefficients." The table simply presents the OLS estimate of $\beta$, standard error of the estimate $SE(\beta)$, the "distance" that $\beta$ is from 0 on the Normal$(0, SE(\beta))$ distribution, and the probability of observing a $\beta$ that far away from 0. Forgive me for the basic statistics review; I can't tell if this is what you are asking for.

Proper OLS-estimated regression modeling (which is what the lm command runs) requires several assumptions, and these diagnostic plots are designed to test them.

The "Residuals vs Fitted" and "Scale-Location" charts are essentially the same, and show if there is a trend to the residuals. OLS models require that the residuals be "identically and independently distributed," that their distribution does not change substantially for different values of $x$. None of your charts is really satisfactory on this regard. If this assumption is not met, your $\beta$ estimates will still be good, but your $t$-statistics, and corresponding $p$-values, are garbage.

Another assumption is that the errors are approximately normally distributed, which is what the Q-Q plot allows you to see. Again, none of your plots really satisfies me in this regard. The consequences of this assumption not being met are the same as above ($\beta$'s good, $t$'s worthless).

The "outliers" principle is actually not an assumption of OLS regression. But if you have outliers in certain locations, your $\beta$ parameters will be unduly influenced by them. In this case, both your $\beta$ and $t$ measurements are useless. You can remove an influential observation from a data frame by identifying its row number and issuing the command

data <- data[-offending.row,]

Where offending.row is the number of the row you want to eliminate. The R diagnostic plots label the row numbers of potential outliers.

I don't know what kind of data you have, but you should be very careful about eliminating observations that you don't like. You should instead ask yourself how that observation became this way. If it is due to measurement error, by all means discard it. If not, then is this observation a part of the system you are trying to model? If so, you should keep it in and adapt for it in other ways.

I have two suggestions for your analysis. First, try to use GLS estimators. This method assigns weights to your observations to correct for heteroskedasticity, outliers, and some degree of non-normality. The R command for this is gls().

But it seems from your plots that your data are restricted in some ways. In particular Test-P seems like a variable that is either 1 or 0, or restricted to that range. For such a variable, you may want to look at binary logit or probit models, available with the command

glm(model, family=binomial(link="logit"))

If your data is censored at 0 but not on the upper end, a tobit model is what you want, tobit() from the AER package looks like the right thing (I've never run a tobit model, I have just looked at it theoretically).

Finally, predictions are done with the predict() function. If you want to perturb your data afterwards (to create a distribution of possible predictions), the best way I know of it to add a random number to the prediction. Using the example above,

#base prediction
pred.values <- predict(lm.D9)
# get standard error of residuals
SER <- (summary(lm.D9)$sigma)^2
#perturbations
pert <- rnorm(length(pred.values), mean=0, sd=SER)
SIMULATION.VALUES <- pred.values + pert

You can get multiple alternate simulations by repeating the last two steps. Good luck.

Best Answer

Related Solutions

Logistic – Why Exact Singularity When Adding Dummy Variable to Multinomial Logistic Using R mlogit?

Linear Regression in R – How to Interpret Model Diagnostics

Related Question