3
$\begingroup$

I'm trying to run an ordinal logistic regression model with 5 IV's and a DV with 4 levels. I'm using the function polr from the MASS package in R. My data consists of 46 observations and all the IV's are continuous.

After fitting the model, I run the summary command, but I get the following warning: In sqrt(diag(vc)) : NaNs produced

After checking the variance-covariance matrix I find that the calculated variance for 3 of the IV's are negative (and very close to zero). Hence, R is only able to calculate SE's for 2 IV's. One of these IV's is significant with a p-value of roughly 0,008.

After reading some other threads, a possible explanation for the negative variance seems to be that the true variance of the beta is zero, or very close to zero. Is this a plausible explanation in a situation like this?

Nevertheless, I still have two question that I'm wondering about:

  1. What would be a good approach for dealing with NaNs in a model like this (or in any type of regression model for that matter)? When I was a student in statistics I remember my teacher (it was a course on analysis of variance) telling us that one way of dealing with a negative variance would be to set it equal to zero, is that a reasonable approach here? That would however make it impossible to calculate a p-value, so I guess you would then assume the p-value to be basically zero as well?
  2. Considering that my model produces NaNs, does that mean that the calculated p-values for my other two IV's are wrong or can I trust them? My guess is that the NaNs doesn't affect the two calculated SE's, but I am not very sure about this.

Any help with this would be very appreciated!

Model:

polr(DV ~ IV1 + IV2 + IV3 + IV4 + IV5, data = df, Hess=TRUE)

Data:

structure(list(DV = structure(c(1L, 2L, 3L, 4L, 4L, 3L, 3L, 2L, 
2L, 2L, 3L, 4L, 4L, 4L, 3L, 3L, 2L, 4L, 2L, 3L, 1L, 4L, 3L, 3L, 
3L, 1L, 4L, 2L, 3L, 4L, 3L, 3L, 2L, 3L, 3L, 2L, 4L, 2L, 4L, 3L, 
4L, 2L, 3L, 3L, 2L, 3L), .Label = c("1", "2", "3", "4"), class = "factor"), 
    IV1 = c(31.4471225882832, 26.7995926461918, 6.15597391098257, 
    32.2641357244643, 5.22243966268843, 50.905131621076, 15.9862670353658, 
    20.3080359821474, 8.65788372462158, 18.505603942379, 6.29152570809156, 
    24.9245052889569, 29.7359693719515, 21.859516176042, 22.8111791071324, 
    23.2948895757526, 7.48299489410315, 25.7653834495302, 11.4419627627679, 
    11.6797724921892, 18.9905104994504, 12.8525067743421, 23.7294286791083, 
    17.7572582793217, 6.51496922762868, 53.5047647614377, 19.7405021265905, 
    27.7157343854444, 20.2388180530257, 9.03173525976901, 15.3013141742333, 
    12.7120495663999, 15.7178643351835, 10.0053600142934, 16.6321423711387, 
    11.8383765639605, 11.4583741140051, 22.9047246959457, 7.72280996279261, 
    20.3566081932521, 7.17681029058905, 8.27263288950658, 27.9236215467482, 
    14.8812599466755, 11.9252118854612, 13.2169871327123), IV2 = c(5.60931899641577, 
    4.99118165784832, 5.40925266903915, 5.11764705882353, 4.89913544668588, 
    5, 5.07709251101322, 4.84575835475578, 3.73287671232877, 
    5.30413625304136, 5.52755905511811, 4.97975708502024, 5.02013422818792, 
    4.81481481481481, 4.69059405940594, 5.07853403141361, 4.31137724550898, 
    5.30913978494624, 4.73282442748092, 4.7065868263473, 5.18518518518519, 
    4.44909344490934, 4.55344070278185, 5.33557046979866, 5.37695590327169, 
    4.96981891348088, 4.60559796437659, 5.21172638436482, 5.10909090909091, 
    4.96318114874816, 5.19753086419753, 4.83164983164983, 4.62227912932138, 
    4.54188481675393, 4.79136690647482, 4.99156829679595, 4.67706013363029, 
    4.2027027027027, 4.81313703284258, 4.7979797979798, 4.29109159347553, 
    4.69798657718121, 4.48275862068965, 4.81012658227848, 4.86254295532646, 
    4.85776805251641), IV3 = c(18.4117647058824, 16.6470588235294, 
    17.8823529411765, 17.4, 18.8888888888889, 16.6875, 18.44, 
    18.85, 16.7692307692308, 19.8181818181818, 18.4736842105263, 
    17.5714285714286, 17, 17.6428571428571, 18.047619047619, 
    18.4761904761905, 16, 19.75, 18.0833333333333, 14.0357142857143, 
    18.1176470588235, 17.7222222222222, 16.3684210526316, 17.6666666666667, 
    18, 17.6428571428571, 16.4545454545455, 17.7777777777778, 
    16.5294117647059, 17.7368421052632, 18.304347826087, 15.9444444444444, 
    18.05, 17.35, 16.65, 17.4117647058824, 17.5, 14.8095238095238, 
    15.7407407407407, 16.7647058823529, 19, 17.5, 16.25, 12.6666666666667, 
    18.8666666666667, 16.4444444444444), IV4 = c(36.1498372555604, 
    29.4083258687817, 30.8237670147563, 30.6498598318535, 33.4609173395659, 
    33.8926923968973, 38.2607303267824, 35.6276276774419, 28.0350077556936, 
    29.8014773776547, 34.7520103344474, 34.2427136374588, 28.4133063999497, 
    30.8454122703714, 25.3434749682913, 31.7741243673625, 26.4398342610694, 
    30.0799370487075, 25.6232533845839, 22.8932982631499, 27.5351660757584, 
    40.9484363905134, 29.9603223354011, 41.4730169772058, 46.6645863988013, 
    25.463413385222, 38.0071453548843, 37.8674109498922, 25.1711130821037, 
    31.105850465851, 28.6739482669466, 33.4120721890402, 35.7158076753924, 
    29.1595387321196, 35.8018355105914, 33.826314306969, 31.3280107842738, 
    28.6294396753682, 26.8677745442452, 33.1564492854815, 23.3944257683981, 
    32.3845970978919, 21.6765100540051, 32.8688639919967, 29.5649229611823, 
    29.9880309522896), IV5 = c(4338.66666666667, 5806, 2450, 
    2624, 2259, 2831.5, 2879.66666666667, 4106.5, 4867, 8250, 
    3824.8, 3431, 1651.2, 2902.25, 3118, 1857.66666666667, 8377.2, 
    2022.33333333333, 4010, 5754.6, 5641, 4042.25, 5745.4, 5008.5, 
    6145.5, 3233.5, 2007.6, 2475.75, 1884.25, 6938.25, 3475, 
    6575.83333333333, 5564, 5474.8, 5665.75, 6787.66666666667, 
    3560.14285714286, 5707.4, 6354.25, 6548.75, 6883.66666666667, 
    7444.5, 6051.33333333333, 2750.5, 4788, 7105.2)), row.names = c(NA, 
-46L), class = c("tbl_df", "tbl", "data.frame"))
$\endgroup$
2
  • 1
    $\begingroup$ Apparently we are talking about 6 columns of 48 observations each. So you could easily copy the dput of the original data here into the forum so people could reproduce the problem which is usually a great help in finding solutions. $\endgroup$
    – Bernhard
    Dec 29, 2020 at 11:08
  • $\begingroup$ The question is updated now $\endgroup$
    – Jonas8
    Dec 29, 2020 at 12:43

1 Answer 1

1
$\begingroup$

Thank you for adding the data. I tried rstanarm::stan_polr only to find the coefficients of three of the five predictors to be 0.0. Then I thought of a scale effect and apparently the large values in IV5 are the problem:

> summary(df)
 DV          IV1              IV2             IV3             IV4             IV5      
 1: 3   Min.   : 5.222   Min.   :3.733   Min.   :12.67   Min.   :21.68   Min.   :1651  
 2:12   1st Qu.:11.446   1st Qu.:4.692   1st Qu.:16.65   1st Qu.:28.47   1st Qu.:2885  
 3:19   Median :16.309   Median :4.860   Median :17.61   Median :30.83   Median :4563  
 4:12   Mean   :18.254   Mean   :4.878   Mean   :17.35   Mean   :31.46   Mean   :4592  
        3rd Qu.:23.197   3rd Qu.:5.101   3rd Qu.:18.11   3rd Qu.:34.16   3rd Qu.:5990  
        Max.   :53.505   Max.   :5.609   Max.   :19.82   Max.   :46.66   Max.   :8377  

If you scale the values of IV5 the functions runs with no problems:

> summary(polr(DV ~ IV1 + IV2 + IV3 + IV4 + scale(IV5), data = df, Hess=TRUE))
Call:
polr(formula = DV ~ IV1 + IV2 + IV3 + IV4 + scale(IV5), data = df, 
    Hess = TRUE)

Coefficients:
              Value Std. Error t value
IV1        -0.06206    0.03176 -1.9539
IV2        -0.85094    1.02221 -0.8325
IV3        -0.08293    0.22607 -0.3668
IV4         0.02918    0.05853  0.4986
scale(IV5) -0.94037    0.32998 -2.8498

Intercepts:
    Value   Std. Error t value
1|2 -8.8255  4.9492    -1.7832
2|3 -6.6140  4.8206    -1.3720
3|4 -4.5695  4.7956    -0.9529

Residual Deviance: 104.1743 
AIC: 120.1743 

Simply dividing IV5 values by 1000 works as well:

> summary(polr(DV ~ IV1 + IV2 + IV3 + IV4 + I(IV5/1000), data = df, Hess=TRUE))
Call:
polr(formula = DV ~ IV1 + IV2 + IV3 + IV4 + I(IV5/1000), data = df, 
    Hess = TRUE)

Coefficients:
               Value Std. Error t value
IV1         -0.06205    0.03176 -1.9536
IV2         -0.84851    1.02206 -0.8302
IV3         -0.08289    0.22606 -0.3667
IV4          0.02918    0.05853  0.4985
I(IV5/1000) -0.49940    0.17530 -2.8488

Intercepts:
    Value    Std. Error t value 
1|2 -11.1059   5.2975    -2.0964
2|3  -8.8946   5.1598    -1.7238
3|4  -6.8500   5.1160    -1.3389

Residual Deviance: 104.1744 
AIC: 120.1744 
$\endgroup$
3
  • $\begingroup$ Thank you! Great solution. Just one thing I'm wondering about. Changing the scale of IV5 changed the p-value for IV1 from 0.008 to slightly above 0.05 so it's no longer significant. I guess a model without NaNs is more reliable than one that produces NaNs, but do you have any idea what might cause such a big shift in the p-value? $\endgroup$
    – Jonas8
    Jan 1, 2021 at 12:50
  • 1
    $\begingroup$ No guaranties -- but cor(df[,2:6]) shows non-zero correlation between IV1 and IV5. So some of the variance of DV that can be explained via IV1 can also alternatively be explained via IV5, once IV5 is in the game. $\endgroup$
    – Bernhard
    Jan 1, 2021 at 17:38
  • $\begingroup$ Makes sense! Thanks a lot! $\endgroup$
    – Jonas8
    Jan 1, 2021 at 20:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.