if (FALSE) {
library(agridat)
data(foulley.calving)
dat <- foulley.calving
## Plot
d2 <- transform(dat,
                age=ordered(age, levels=c("0.0-2.0","2.0-2.5","2.5-3.0",
                                          "3.0-3.5","3.5-4.0",
                                          "4.0-4.5","4.5-5.0","5.0-8.0","8.0+")),
                score=ordered(score, levels=c('S1','S2','S3')))
libs(reshape2)
d2 <- acast(dat, sex+age~score, value.var='count')
d2 <- prop.table(d2, margin=1)
libs(lattice)
thm <- simpleTheme(col=c('skyblue','gray','pink'))
barchart(d2, par.settings=thm, main="foulley.calving",
         xlab="Frequency of calving difficulty", ylab="Calf gender and dam age",
         auto.key=list(columns=3, text=c("Easy","Assited","Difficult")))
## Ordinal multinomial model
libs(ordinal)
m2 <- clm(score ~ sex*age, data=dat, weights=count, link='probit')
summary(m2)
##   Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## sexM             0.500605   0.015178  32.982  < 2e-16 ***
## age2.0-2.5      -0.237643   0.013846 -17.163  < 2e-16 ***
## age2.5-3.0      -0.681648   0.018894 -36.077  < 2e-16 ***
## age3.0-3.5      -0.957138   0.018322 -52.241  < 2e-16 ***
## age3.5-4.0      -1.082520   0.024356 -44.446  < 2e-16 ***
## age4.0-4.5      -1.146834   0.022496 -50.981  < 2e-16 ***
## age4.5-5.0      -1.175312   0.028257 -41.594  < 2e-16 ***
## age5.0-8.0      -1.280587   0.016948 -75.559  < 2e-16 ***
## age8.0+         -1.323749   0.024079 -54.974  < 2e-16 ***
## sexM:age2.0-2.5  0.003035   0.019333   0.157  0.87527    
## sexM:age2.5-3.0 -0.076677   0.026106  -2.937  0.00331 ** 
## sexM:age3.0-3.5 -0.080657   0.024635  -3.274  0.00106 ** 
## sexM:age3.5-4.0 -0.135774   0.032927  -4.124 3.73e-05 ***
## sexM:age4.0-4.5 -0.124303   0.029819  -4.169 3.07e-05 ***
## sexM:age4.5-5.0 -0.198897   0.038309  -5.192 2.08e-07 ***
## sexM:age5.0-8.0 -0.135524   0.022804  -5.943 2.80e-09 ***
## sexM:age8.0+    -0.131033   0.031852  -4.114 3.89e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Threshold coefficients:
##       Estimate Std. Error z value
## S1|S2  0.82504    0.01083   76.15
## S2|S3  1.52017    0.01138  133.62
## Note 1.52017 - 0.82504 = 0.695 matches Foulley's '2-3' threshold estimate
predict(m2) # probability of each category
}
Run the code above in your browser using DataLab