4 added graphical estimate of underlying distribution from proportional odds model edited Feb 19 '16 at 12:38 Frank Harrell 58.3k44 gold badges118118 silver badges252252 bronze badges 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='') 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='') 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 edited Feb 16 '16 at 12:52 Frank Harrell 58.3k44 gold badges118118 silver badges252252 bronze badges 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 edited Feb 16 '16 at 12:45 Frank Harrell 58.3k44 gold badges118118 silver badges252252 bronze badges 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 answered Feb 15 '16 at 19:48 Frank Harrell 58.3k44 gold badges118118 silver badges252252 bronze badges