4 added graphical estimate of underlying distribution from proportional odds model
source | link

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex. New: a plot of the estimated underlying conditional distribution of y is added.

require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

# Show intercepts as a function of y to estimate the underlying
# conditional distribution.  Result: more uniform than Gaussian
alphas <- coef(f)[1 : num.intercepts(f)]
yunique <- f$yunique[-1]
par(mfrow=c(1,2))
plot(yunique, alphas)
# Compare to distribution of residuals
plot(ecdf(resid(ols(y ~ cohort * sex, data=d))), main='')

Intercepts vs. y and residuals from OLS

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, summarize(y, llist(cohort, sex), smean.cl.normal))
  cohort    sex         y      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex.

require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, summarize(y, llist(cohort, sex), smean.cl.normal))
  cohort    sex         y      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex. New: a plot of the estimated underlying conditional distribution of y is added.

require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

# Show intercepts as a function of y to estimate the underlying
# conditional distribution.  Result: more uniform than Gaussian
alphas <- coef(f)[1 : num.intercepts(f)]
yunique <- f$yunique[-1]
par(mfrow=c(1,2))
plot(yunique, alphas)
# Compare to distribution of residuals
plot(ecdf(resid(ols(y ~ cohort * sex, data=d))), main='')

Intercepts vs. y and residuals from OLS

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, summarize(y, llist(cohort, sex), smean.cl.normal))
  cohort    sex         y      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924
3 improved last part of example
source | link
require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, tapplysummarize(y, listllist(cohort, sex), meansmean.cl.normal))
  cohort    sex     male    femaley      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924
require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means:
with(d, tapply(y, list(cohort, sex), mean))
         male    female
one  2.055708  5.024132
two 12.711545 15.348114
require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, summarize(y, llist(cohort, sex), smean.cl.normal))
  cohort    sex         y      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924
2 added full example
source | link

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex.

require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means:
with(d, tapply(y, list(cohort, sex), mean))
         male    female
one  2.055708  5.024132
two 12.711545 15.348114

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex.

require(rms)
d1 <- data.frame(cohort='one', sex='male', y=c(.476,
.84,
1.419,
0.4295,
0.083,
2.9595,
4.20125,
1.6605,
3.493,
5.57225,
0.076,
3.4585))
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
4.591,
3.138,
2.699,
6.622,
6.8795,
5.5925,
1.6715,
4.92775,
6.68525,
4.25775,
8.677))
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
16.252,
15.30175,
8.66325,
15.6935,
16.214,
4.056,
8.316,
17.95725,
13.644,
15.76475))
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
22.22775,
18.00466667,
12.80925,
16.15425,
14.88133333,
12.0895,
16.5335,
17.68925,
15.00425,
12.149))
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

anova(f)
            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means:
with(d, tapply(y, list(cohort, sex), mean))
         male    female
one  2.055708  5.024132
two 12.711545 15.348114
1
source | link