2
$\begingroup$

I have conducted a logistic regression.

Model <- logistf(A ~ C, family ="binomial"(link = "logit"), data=Data)

The dependent variable is binary (1,0) The sole independent variable is ordinal with 3 levels. I've included the relevant parts of the summary below.


                   coef       se(coef)    Chisq           p
(Intercept)       -0.05       0.36        0.023       0.87
C.L               1.6         0.74        6.67        0.009
C.Q              -0.2         0.51        0.31        0.57

Please note here L and Q refer to linear and quadratic effects. This is a polynomial contrast. The independent variable is ordinal with 3 levels low<medium<high.

How can I interpret this summary? As I understand it the log odds of a "success" linearly increase by 1.6 with the independent variable "C". However what are the log odds of a success when my independent is at its lowest level, if the intercept is insignificant?

EDIT to include info from comment

> contr.poly(3)
                .L         .Q
[1,] -7.071068e-01  0.4082483
[2,] -7.850462e-17 -0.8164966
[3,]  7.071068e-01  0.4082483

EDIT 2 More info from comments

> head(Data)
                    C           A
1                Medium         0
2                Medium         0
3                  High         1
4                  High         1
5                Medium         1
6                Medium         0

#note C is 3 levels - Low, medium, high. Low simply isn't shown at the top of the dataframe.

Data$C <- factor(Data$C, levels = c("Low", "Medium", "High"), ordered = TRUE)
Data$A <- Data$A <- as.integer(Data$A)

Model <- logistf(A ~ C, family ="binomial"(link = "logit"), data=Data)

summary(Model)
                         coef  se(coef) lower 0.95 upper 0.95      Chisq           p
(Intercept)       -0.05228334 0.3693310 -0.8607010  0.5764470 0.02318317 0.878981724
C.L                1.60842992 0.7458429  0.3711041  3.2751171 6.67898157 0.009755635
C.Q               -0.26863239 0.5120064 -1.3463216  0.6337436 0.31863128 0.572431359

Likelihood ratio test=7.75 on 2 df, p=0.020, n=82
Wald test = 6.036922 on 2 df, p = 0.048

Covariance-Matrix:
           [,1]       [,2]       [,3]
[1,]  0.1364054 -0.2044094  0.1200786
[2,] -0.2044094  0.5562817 -0.2503493
[3,]  0.1200786 -0.2503493  0.2621505

EDIT 3 More info from comments


> Model$predict[1:5]
[1] 0.5416667 0.5416667 0.7261905 0.7261905 0.5416667
Improve this question
$\endgroup$
13
  • $\begingroup$ Apologies. Just realised that I misinterpreted what I read before, so I have removed that sentence. I was reading this support.sas.com/kb/23/136.html $\endgroup$
    – Harry
    Commented Aug 1, 2020 at 18:50
  • $\begingroup$ See also stats.stackexchange.com/questions/260209/… $\endgroup$
    – kjetil b halvorsen ♦
    Commented Aug 1, 2020 at 19:31
  • $\begingroup$ Thanks @kjetilbhalvorsen I'm afraid my actual understanding of the maths here is pretty poor. Would you be able to describe what I can actually conclude from my data here? $\endgroup$
    – Harry
    Commented Aug 1, 2020 at 20:33
  • 1
    $\begingroup$ It is telling you that the data is consistent with a linear effects. Going from low to medium has the same effect on the log-odds as going from medium to high, which are half of the effect you get from going from low to high. Cubic effects are not well supported. How the intercept fits into all that, I don't really know. $\endgroup$
    – dimitriy
    Commented Aug 1, 2020 at 22:09
  • 1
    $\begingroup$ This can be resolved if you show us what the first few rows of Data look like, show us the actual code that you typed, with the output that you want help in interpreting that results from the code you actually typed. We have wasted a fair bit of time here already trying to fill the gaps between what you say you are doing, what you think you are doing, and what you actually are doing. $\endgroup$
    – dimitriy
    Commented Aug 1, 2020 at 22:47

2 Answers 2

Reset to default
3
$\begingroup$

You defined C as an ordered variable, which it is, so logistf automatically converts C to an orthogonal polynomial representation, which makes things harder to interpret. Why does R do this?

Orthogonal polynomial encoding is a form of trend analysis in that it is looking for the linear, quadratic and cubic trends in the effect of a categorical variable. Here you only have three values of C, so we can only do linear and cubic. This type of coding system should be used only with an ordinal variable in which the levels are equally spaced (which may or may not be the case in your data, the coding is evenly spaced, but the underlying variables may not be). In R it is not necessary to compute these values since this contrast can be obtained for any categorical variable by using the contr.poly function that I asked you to use. This the default contrast used for ordered factor variables. This makes it easier to find a functional form that will allow you to model your categorical variables as if they were more continuous. This can help impose more structure on the effects, compared to something more flexible, where every category can have its own effect.

The fact that C.L coefficient is positive and significant, tells you that the data is consistent with a linear trend of C on the log odds: each time C goes up by 1, log odds go up by the same amount. The fact that quadratic is negative, means you have a deceleration in that trend (slower growth), but it is not significant.

Here is a reproducible example with some comments interspersed (reproducible code will be below):

    > library(logistf)
    > library(boot)
    > 
    > # Calculate your OP model predictions for C=1,2,3
    > inv.logit(-0.05228334 + 1.60842992*( -7.071068e-01) + 
               (-0.26863239)*0.4082483)
    [1] 0.2142856
    > inv.logit(-0.05228334 + 1.60842992*(-7.850462e-17) + 
                (-0.26863239)*-0.8164966)
    [1] 0.5416667
    > inv.logit(-0.05228334 + 1.60842992*(7.071068e-01) + 
               (-0.26863239)*0.4082483)
    [1] 0.7261905

As you can see, this matches 2 of the 3 predictions that you reported. Now for the toy example:

    > # Toy Example
    > raw <- 'A C   cL  cQ
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 3   .7071068    .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 3   .7071068    .4082483
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 3   .7071068    .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 3   .7071068    .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 3   .7071068    .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 3   .7071068    .4082483
    + 0 1   -.7071068   .4082483
    + 0 3   .7071068    .4082483
    + 0 3   .7071068    .4082483
    + 0 3   .7071068    .4082483
    + 0 3   .7071068    .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 1   -.7071068   .4082483
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 2   -7.85e-17   -.8164966
    + 0 3   .7071068    .4082483
    + 1 1   -.7071068   .4082483
    + 1 2   -7.85e-17   -.8164966
    + 1 3   .7071068    .4082483
    + 1 2   -7.85e-17   -.8164966
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 2   -7.85e-17   -.8164966
    + 1 2   -7.85e-17   -.8164966
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 1   -.7071068   .4082483
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 1   -.7071068   .4082483
    + 1 3   .7071068    .4082483
    + 1 1   -.7071068   .4082483
    + 1 2   -7.85e-17   -.8164966
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 3   .7071068    .4082483
    + 1 1   -.7071068   .4082483
    + '
    
    > Data <- read.table(text = raw, header = TRUE)
    > Data$C <- factor(Data$C, levels = c(1, 2, 3), ordered = TRUE)
    > 
    > # Show the orthogonal polynomial representation of C
    > unique(Data[,c("C","cL","cQ")])
      C            cL         cQ
    1 2 -7.850000e-17 -0.8164966
    2 1 -7.071068e-01  0.4082483
    7 3  7.071068e-01  0.4082483
    > contr.poly(3)
                    .L         .Q
    [1,] -7.071068e-01  0.4082483
    [2,] -7.850462e-17 -0.8164966
    [3,]  7.071068e-01  0.4082483

As you can see, the values returned by contr.poly(3) match what is in my data. We will fit the model, but not a Firth logit, just a normal logit (but that should not really matter; also not clear to me what you are doing at home):

    > # Plain Logit Model with Orthogonal Polynomial
    > op_model<-logistf(A~C, family=binomial(link="logit"), data=Data, 
                       firth=FALSE, pl=TRUE)
    > summary(op_model)
    logistf(formula = A ~ C, data = Data, pl = TRUE, firth = FALSE, 
        family = binomial(link = "logit"))
    
    Model fitted by Standard ML
    Confidence intervals and p-values by Profile Likelihood Profile Likelihood Profile Likelihood 
    
                      coef  se(coef) lower 0.95 upper 0.95     Chisq            p
    (Intercept) -0.9098647 0.2733188 -1.4739880 -0.3932878 12.372954 0.0004355978
    C.L          1.1091790 0.4580922  0.2432744  2.0621977  6.387326 0.0114938060
    C.Q          0.5206835 0.4882321 -0.4030296  1.5437748  1.191198 0.2750880661
    
    Likelihood ratio test=7.786674 on 2 df, p=0.02037724, n=74
    Wald test = 7.538426 on 2 df, p = 0.02307021
    
    Covariance-Matrix:
                [,1]        [,2]        [,3]
    [1,]  0.07470317 -0.01678486 -0.01164408
    [2,] -0.01678486  0.20984848 -0.02055717
    [3,] -0.01164408 -0.02055717  0.23837055

You can even see what R is actually doing to C with this bit of code:

    > head(model.matrix(logistf(Data$A~Data$C, family=binomial(link="logit"), firth=FALSE, pl=TRUE)))
      (Intercept)      Data$C.L   Data$C.Q
    1           1 -7.850462e-17 -0.8164966
    2           1 -7.071068e-01  0.4082483
    3           1 -7.850462e-17 -0.8164966
    4           1 -7.850462e-17 -0.8164966
    5           1 -7.071068e-01  0.4082483
    6           1 -7.071068e-01  0.4082483

Now we will do the same thing using the cL and cQ variables I entered myself. We should see the same model summary as above:

    > # Plain Logit Model with Orthogonal Polynomial by Hand
    > op_model_by_hand<-logistf(A~cL+cQ, family=binomial(link="logit"), 
                               data=Data, firth=FALSE, pl=TRUE)
    > summary(op_model_by_hand)
    logistf(formula = A ~ cL + cQ, data = Data, pl = TRUE, firth = FALSE, 
        family = binomial(link = "logit"))
    
    Model fitted by Standard ML
    Confidence intervals and p-values by Profile Likelihood Profile Likelihood Profile Likelihood 
    
                      coef  se(coef) lower 0.95 upper 0.95     Chisq            p
    (Intercept) -0.9098647 0.2733188 -1.4739880 -0.3932878 12.372954 0.0004355978
    cL           1.1091789 0.4580922  0.2432744  2.0621977  6.387326 0.0114938060
    cQ           0.5206835 0.4882321 -0.4030296  1.5437748  1.191198 0.2750880661
    
    Likelihood ratio test=7.786674 on 2 df, p=0.02037724, n=74
    Wald test = 7.538426 on 2 df, p = 0.02307021
    
    Covariance-Matrix:
                [,1]        [,2]        [,3]
    [1,]  0.07470317 -0.01678486 -0.01164408
    [2,] -0.01678486  0.20984847 -0.02055717
    [3,] -0.01164408 -0.02055717  0.23837054

This way to represent C is somewhat unintuitive. The specific values of cL and cQ themselves have no meaning to us because they were computed by R to make all the contrasts linearly independent of one another. This kind of encoding is just one of at least 9 possible ways to encode categorical variables. For OP, it probably does not make too much sense to think hard about the intercept insignificance. And it certainly does not mean that you should take it out. You should never decide which variables to omit from the model based on statistical significance only.

This is how you might do the one-hot encoding that lends itself to easier interpretation:

    > # Remove the constraint that C is ordered to get one-hot encoded version
    > Data$C <- factor(Data$C, levels = c(1, 2,3), ordered = FALSE)
    > oh_model<-logistf(A~C, family=binomial(link="logit"), data=Data, 
                       firth=FALSE, pl=TRUE)
    > summary(oh_model)
    logistf(formula = A ~ C, data = Data, pl = TRUE, firth = FALSE, 
        family = binomial(link = "logit"))
    
    Model fitted by Standard ML
    Confidence intervals and p-values by Profile Likelihood Profile Likelihood Profile Likelihood 
    
                      coef  se(coef) lower 0.95 upper 0.95       Chisq            p
    (Intercept) -1.4816045 0.4954337 -2.5755111 -0.5899917 11.55500406 0.0006756716
    C2           0.1466035 0.7057522 -1.2656065  1.5604230  0.04313599 0.8354694171
    C3           1.5686159 0.6478402  0.3440419  2.9163880  6.38732576 0.0114938060
    
    Likelihood ratio test=7.786674 on 2 df, p=0.02037724, n=74
    Wald test = 7.538426 on 2 df, p = 0.02307021
    
    Covariance-Matrix:
               [,1]       [,2]       [,3]
    [1,]  0.2454545 -0.2454545 -0.2454545
    [2,] -0.2454545  0.4980861  0.2454545
    [3,] -0.2454545  0.2454545  0.4196970

Personally, this is the model I would use myself with your data. You can use it to calculate the baseline odds from the intercept, and the exponentiated coefficients will give you the multiplicative effect on those baseline odds.

I don't think there is any way to interpret the OP model in the same way. I also don't think it makes sense to treat C as an ordered variable here and try to find trends when it only has three values.

Code:

    library(logistf)
    library(boot)
    
    # Calculate your OP model predictions for C=1,2,3
    inv.logit(-0.05228334 + 1.60842992*( -7.071068e-01) + (-0.26863239)*0.4082483)
    inv.logit(-0.05228334 + 1.60842992*(-7.850462e-17) + (-0.26863239)*-0.8164966)
    inv.logit(-0.05228334 + 1.60842992*(7.071068e-01) + (-0.26863239)*0.4082483)
    
    # Toy Example
    raw <- 'A   C   cL  cQ
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   3   .7071068    .4082483
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   3   .7071068    .4082483
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   3   .7071068    .4082483
    0   2   -7.85e-17   -.8164966
    0   3   .7071068    .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   3   .7071068    .4082483
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   3   .7071068    .4082483
    0   1   -.7071068   .4082483
    0   3   .7071068    .4082483
    0   3   .7071068    .4082483
    0   3   .7071068    .4082483
    0   3   .7071068    .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   1   -.7071068   .4082483
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   2   -7.85e-17   -.8164966
    0   3   .7071068    .4082483
    1   1   -.7071068   .4082483
    1   2   -7.85e-17   -.8164966
    1   3   .7071068    .4082483
    1   2   -7.85e-17   -.8164966
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   2   -7.85e-17   -.8164966
    1   2   -7.85e-17   -.8164966
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   1   -.7071068   .4082483
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   1   -.7071068   .4082483
    1   3   .7071068    .4082483
    1   1   -.7071068   .4082483
    1   2   -7.85e-17   -.8164966
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   3   .7071068    .4082483
    1   1   -.7071068   .4082483
    '
    
    Data <- read.table(text = raw, header = TRUE)
    Data$C <- factor(Data$C, levels = c(1, 2, 3), ordered = TRUE)
    
    # Show the orthogonal polynomial representation of C
    unique(Data[,c("C","cL","cQ")])
    contr.poly(3)
    
    # Plain Logit Model with Orthogonal Polynomial
    op_model<-logistf(A~C, family=binomial(link="logit"), data=Data, firth=FALSE, pl=TRUE)
    summary(op_model)
    head(model.matrix(logistf(Data$A~Data$C, family=binomial(link="logit"), firth=FALSE, pl=TRUE)))
    
    # Plain Logit Model with Orthogonal Polynomial by Hand
    op_model_by_hand<-logistf(A~cL+cQ, family=binomial(link="logit"), data=Data, firth=FALSE, pl=TRUE)
    summary(op_model_by_hand)
    
    # Remove the constraint that C is ordered to get one-hot encoded version
    Data$C <- factor(Data$C, levels = c(1, 2,3), ordered = FALSE)
    oh_model<-logistf(A~C, family=binomial(link="logit"), data=Data, firth=FALSE, pl=TRUE)
    summary(oh_model)
Improve this answer
$\endgroup$
2
$\begingroup$

The predictor is categorical (in fact ordinal) with three levels. Since you are representing it with a linear and quadratic term (I suppose), it matters which numerical values you have used to represent it. Three parameters gives a saturated model, and the intercept is an intrinsic part of that description, so it does not really make much sense to look at tests of the coefficients individually. If you used a dummy variable representation, you would get also three parameters, with identical fitted values (but different coefficients.)

The most practical way to represent your results would be by plotting the estimated probabilities as a function of the ordinal variable, say with numerical values 1,2,3 (if those are what you used.) Since you didn't tell us I cannot plot it for you.

If you want to look at individual hypothesis tests, start with the highest polynomial order, your C.Q. That seems to be insignificant, so you can consider a reduced linear model. Since C.L have a low p-value, you should stop the reduction there. You should never consider to remove the intercept!

Improve this answer
$\endgroup$
4
  • $\begingroup$ The numerical values used are 1, 2, 3. Model$predict[1:3] yields 0.5416667 0.5416667 0.7261905. Is this what you need? Sorry if it's not, I'm a little out of my depth. $\endgroup$
    – Harry
    Commented Aug 1, 2020 at 21:40
  • $\begingroup$ Can you please add that information to the Q as an edit, togetger with relevant code, so we can be sure what you are doing? $\endgroup$
    – kjetil b halvorsen ♦
    Commented Aug 1, 2020 at 23:18
  • $\begingroup$ I have added the section on Model$predict. I think I have added everything included in this model in the code now. $\endgroup$
    – Harry
    Commented Aug 1, 2020 at 23:29
  • $\begingroup$ That is good, but what you jave wrote in the Q is different from what in the comment .... ? can you correct, and explain? Maybe share (a link to) the data, so we can show some examples in answers, that way it should be clearer. $\endgroup$
    – kjetil b halvorsen ♦
    Commented Aug 2, 2020 at 3:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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