Solved – R: Interpretation ordinal regression with clm

ordinal-datarregression

I have a question regarding the interpretation of an ordinal regression in R with the clm function.

I am new to this and used this formula:

RegADHD <- clm(ADHD$Score ~ ADHD$Group)

Score is a ordered factor from 0-3 (0= no symptoms, 3= maximal symptoms)

Group is a group variable with 3 categories (no exposure, less exposure, severe exposure).

With summary I get this:

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
ADHD$Group3   1.0422     0.4647   2.242   0.0249 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
0|1   1.3612     0.3720   3.659
1|2   2.5735     0.4465   5.764
2|3   3.3582     0.5367   6.258
(2 observations deleted due to missingness)

But I don´t really understand what is the interpretation of this coefficient?

More specifically, I struggle because I get only one value (for group 3).
These are the contrasts of group:

             2    3
    1        0   0
    2        1   0
    3        0   1

To get the odds I did:

round(exp(RegADHD$beta), 3)

This is 2.835

with(RegADHD, table(ADHD$Score, ADHD$Group))

leads to

   1   2   3
0  0   34  30
1  0   7  11
2  0   1   5
3  0   1   5

Does this mean that the probability of having a higher score (in ADHD) increases 2.835 times with being in the severe exposure group, compared with no exposure? But then I don´t know anything about less exposure… I don´t understand why I get only one coefficient!

I hope anyone can help me, I am thankful for every advice!

Best Answer

The odds ratio you have calculated is correct but your interpretation of it is not quite right. You did not tell R that Group was a factor so it has assumed it is a continuous variable with values 0, 1, 2. So the OR you have is for a unit change in Group, from 0 to 1, or from 1 to 2. I suspect that what you want to do is to change the class of Group and re-run. If you do and you still have problems then edit the new information into your question and see if someone can explain further.

With the new information we see that in fact there was nobody in the lowest Group (no exposure) so the coefficient is comparing severe with less, not with no. With that change your interpretation is correct.

Related Solutions

Solved – Lsmeans output for clmm models (R)

First, you should switch to the emmeans package, which is the continuation of lsmeans; the latter will be archived on CRAN soon.

Regarding the instruction to ignore certain results. I invite you to ask yourself what kind of statistical software would output results that should be ignored. The more likely possibility is that the writer of that advice does not understand the output so doesn't want to try to explain it, or disagrees with that method of analysis for an unstated reason (but why then show the methods at all?).

The default results of lsmeans() (or emmeans::emmeans()) are on a latent-variable scale; that latent model asserts that there is a continuous but unobservable response having a logistic distribution with a mean that depends on the predictors, and that there is also a set of cut points that define a set of intervals on the latent scale. In our data, the ordinal response corresponds to which interval each latent response falls in. This latent scale is not completely identifible, so any linear function of that scale is equally valid. The lsmeans() output picks one version of that, and its results are meaningful relative to that scale.

Moreover, the results with mode = "mean.class" are completely interpretable on an absolute scale. Regarding the ordinal values as 1, 2, ..., the mean.scale predictions comprise the mean of the ordinal variable, based on the predicted probabilities of its ordinal distribution.

For an example detailing the analysis of a clm model, see the relevant emmeans vignette: https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html#ordinal. Documentation for the different options available for ordinal models is provided in https://cran.r-project.org/web/packages/emmeans/vignettes/models.html#O. There is also an index of vignette topics where you may find yet more information.

Solved – Standard Error of simple linear regression coefficients

Let a simple linear regression model

$$ y_i = \beta_1 + \beta_2x_i + \epsilon_i $$ from $n$ observations, where $\epsilon_i$ are iid and of same variance $\sigma^2$.

OLS estimators of $\beta_1$ and $\beta_2$ are given by $$ \hat{\beta}_2 = \frac{\sum(x_i-\bar{x})y_i}{\sum(x_i - \bar{x}^2} $$ and $$ \hat{\beta}_1 = \bar{y} - \hat{\beta}_2 \bar{x} $$ where $\bar{x}$ denotes sample mean. From each parameter we only have one value (since we have one sample). We do not need to estimate $\sigma^2$ to compute both $\hat{\beta_1}$ and $\hat{\beta}_2$.
However, it can be estimated with \begin{align*} \hat{\sigma}^2 &= \frac{1}{n-2} \sum( y_i - \hat{y}_i )^2 \\ &= \frac{1}{n-2} \sum( y_i - \hat{\beta}_1 - \hat{\beta}_2 x_i )^2 \end{align*}

But even if we only have one value for each estimator, we have the following results : for OLS estimates $\hat{\beta}_1$ and $\hat{\beta}_2$.

$$ \text{Var}(\hat{\beta}_1) =\frac{1}{n} \frac{\sigma^2 \sum x_i^2}{\sum(x_i - \bar{x})^2} $$ and $$ \text{Var}(\hat{\beta}_2) = \frac{\sigma^2}{\sum(x_i - \bar{x})^2} $$

Moreover,

$$ \text{Cov}(\hat{\beta}_1,\hat{\beta}_2) = - \frac{\sigma^2 \bar{x}}{\sum(x_i - \bar{x})^2} $$ Proof of these results can be found in any textbook on linear regression.

Since we can estimate $\sigma^2$ from the original sample, we can estimate variances of estimators, even if we only have one value of them.

These variances estimates can also be computed from bootstrapping the sample, i.e by taking samples of the original sample and by computing estimators for each sub sample. For example for $\beta_1$, if you take $K$ sub sample then you get a sample $\hat{\beta}_1^{(1)},\dots,\hat{\beta}_1^{(K)}$ from which you can empirically estimate the variance of $\hat{\beta}_1$.

Here is a simple code (in R) showing how the two methods can be used

data<-data.frame(do.call(rbind,lapply(1:500,function(i){
id=i
x<-rexp(1,1)
y<- 1 + 3*x + rnorm(1,0,1)
return(c(id,y,x))
})))
names(data)<-c("id", "y","x")



summary(lm(y~x,data)) ## OLS std estimates from full sample


## by bootstrapping
Nboot<-500
list<-lapply(1:Nboot,function(i){
Ids <- sort(sample(data$id,replace=TRUE))
data.s = data[Ids,]
mod.s<-lm(y~x,data.s)
return(mod.s$coefficients)
})
## Bootstrap variances estimates
var(sapply(list,function(x) x[[1]]))**0.5  # beta1
var(sapply(list,function(x) x[[2]]))**0.5  # beta2

in the command

summary(lm(y~x,data)) ## OLS std estimates

the values in the "Std. Error" column should be close to the values computed in the last two lines

Best Answer

Related Solutions

Solved – Lsmeans output for clmm models (R)

Solved – Standard Error of simple linear regression coefficients

Related Question