postscript("cumlogit.ps") par(mfrow=c(2,2)) ### prop odds model such that category probabilities are ### (1/4,1/4,1/4,1/4) when x=0, and (1/20, ?, ?, ?) when x=1. ### so theta is quartiles of standard logistic distribution theta <- log(c(1/4,2/4,3/4)) - log(c(3/4,2/4,1/4)) ### and beta chosen to yield Pr(Y=1|X=1)=1/20 beta <- log(1/19)-theta[1] print(beta) x <- (0:200)/200 plot(-99,-99, xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Prob(Y=j)") ### j=1 points(x, 1/(1+exp(-(theta[1]+x*beta))), type="l") ### j=2,3 for (j in 2:3) { points(x, 1/(1+exp(-(theta[j]+x*beta))) - 1/(1+exp(-(theta[j-1]+x*beta))) ,type="l") } ### j=4 points(x, 1 - 1/(1+exp(-(theta[3]+x*beta))), type="l") ### prop hazards model such that category probs are ### (1/4,1/4,1/4,1/4) when x=0, and (1/20, ?, ?, ?) when x=1. ### so theta is quartiles of standard extreme value distribution theta <- log(-log(1-c(1/4,2/4,3/4))) ### and beta chosen to yield Pr(Y=1|X=1)=1/20 beta <- log(-log(1-.05)) - theta[1] print(beta) plot(-99,-99, xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Prob(Y=j)") ### j=1 points(x, 1-exp(-exp(theta[1]+x*beta)), type="l") ### j=2,3 for (j in 2:3) { points(x, exp(-exp(theta[j-1]+x*beta)) - exp(-exp(theta[j] +x*beta)), type="l") } ### j=4 points(x, exp(-exp(theta[3]+x*beta)), type="l") #### this time choose beta so that Pr(Y=4|X=1)=0.5 theta <- log(c(1/4,2/4,3/4)) - log(c(3/4,2/4,1/4)) beta <- log(1)-theta[3] print(beta) x <- (0:200)/200 plot(-99,-99, xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Prob(Y=j)") ### j=1 points(x, 1/(1+exp(-(theta[1]+x*beta))), type="l") ### j=2,3 for (j in 2:3) { points(x, 1/(1+exp(-(theta[j]+x*beta))) - 1/(1+exp(-(theta[j-1]+x*beta))) ,type="l") } ### j=4 points(x, 1 - 1/(1+exp(-(theta[3]+x*beta))), type="l") theta <- log(-log(1-c(1/4,2/4,3/4))) beta <- log(-log(1-.5)) - theta[3] print(beta) plot(-99,-99, xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Prob(Y=j)") ### j=1 points(x, 1-exp(-exp(theta[1]+x*beta)), type="l") ### j=2,3 for (j in 2:3) { points(x, exp(-exp(theta[j-1]+x*beta)) - exp(-exp(theta[j] +x*beta)), type="l") } ### j=4 points(x, exp(-exp(theta[3]+x*beta)), type="l") graphics.off()