Probit with indicator variable

dummy probit.r

###### Ordered Categorical DV: #####
#This uses "importance of aid to blacks" as a DV and race of R as the indicator IV.#
     library(MASS)
     white <- c(1,1,0,1,1,1,1,0,1,1,1,0,1,1,1,1,1,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,
1,1,0,1,1,1,1,1,1,1,1,1,0,0,0,1,0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,1,0,1,1,1,1,1,0,1,1,1,0,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,0,1,0,1,1,1,1,0,1,0,0,1,0,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,
1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,1,1,0,1,1,1,1,1,0,1,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,0,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,1,
0,1,1,0,1,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,1,1,0,1,1,0,1,1,1,1,1,0,1,1,1,0,1,0,1,0,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,1,
1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,0,1,1,0,
1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,0,1,0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,1,0,0,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,0,1,1,1,1,1,1,1,1,1,0,1,
0,1,1,1,1,1,0,0,0,1,0,1,1,1,1,0,0,1,0,1,1,0,1,1,1,0,1,1,1,1,0,1,0,0,0,1,1,0,1,1,1,1,1,1,0,1,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,
1,1,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,
0,1,1,1,1,1,0,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,0,1,0,1,1,1,0,0,0,1,0,1,1,1,1,1,0,1,1,0,1,0,1,1,0,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,
1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1)

aid<-c(4,1,3,2,4,3,2,4,3,2,1,4,2,NA,2,4,4,4,1,4,2,2,4,2,3,2,4,1,3,2,4,2,4,0,1,2,1,1,3,3,3,1,2,4,2,3,4,3,4,3,4,2,1,
2,2,2,3,2,3,3,1,2,3,2,2,1,2,1,2,3,4,3,3,3,2,2,2,1,3,2,3,4,4,2,4,2,2,3,2,3,NA,2,3,3,2,4,4,4,2,1,2,2,3,0,3,4       ,
4,2,1,1,2,2,4,3,2,4,1,3,0,2,2,4,2,2,4,2,2,2,2,3,2,1,2,3,4,3,0,2,3,2,4,3,2,3,2,1,4,4,2,2,2,1,4,4,2,2,2,4,2        ,
2,3,1,2,4,1,3,3,2,3,1,3,2,4,2,2,2,2,3,2,4,3,2,4,2,3,4,4,3,2,2,3,2,3,2,1,1,3,3,1,4,4,4,4,1,2,2,2,1,2,2,2,2        ,
3,3,2,4,4,3,2,3,2,2,2,2,3,3,1,4,1,2,4,2,2,4,1,2,2,1,2,2,2,2,3,1,3,2,3,4,2,2,3,1,3,2,2,3,2,2,2,3,1,4,2,2,4        ,
3,3,4,4,4,2,4,0,2,1,1,0,3,1,4,2,2,2,3,1,3,2,0,1,2,3,4,3,1,4,3,1,1,2,3,3,1,2,2,3,2,1,2,3,1,1,3,3,1,3,2,4,4        ,
2,2,2,2,1,4,1,0,0,0,2,3,3,2,3,3,4,3,2,3,3,2,2,3,4,1,2,3,2,4,4,2,2,4,3,1,3,3,2,4,4,3,3,2,4,3,1,2,2,2,1,4,3        ,
3,2,2,2,2,2,2,3,4,3,4,1,4,0,3,4,2,2,3,4,2,2,4,2,4,2,1,3,0,3,2,2,2,2,2,3,4,2,3,1,3,1,4,2,1,1,1,3,2,3,3,1,2        ,
2,1,2,2,4,2,2,1,1,1,3,4,4,4,4,2,3,3,3,3,4,3,3,2,2,3,4,2,4,3,1,4,2,NA,3,0,2,4,3,2,3,1,4,3,2,4,1,2,1,2,1,3,3       ,
3,3,4,1,3,0,3,1,1,2,1,1,2,4,3,0,2,1,2,2,1,2,2,2,2,3,4,4,4,2,2,3,3,3,1,4,4,1,4,2,4,4,2,2,4,2,3,2,3,2,1,2,2        ,
2,4,2,3,4,2,2,2,1,4,3,4,2,1,2,3,1,2,2,3,4,2,2,1,1,3,2,3,2,4,3,2,2,2,3,2,3,3,3,4,4,0,2,3,3,NA,2,2,3,4,NA,3,3      ,
4,2,3,1,3,4,2,2,1,4,3,NA,1,3,3,2,4,1,2,3,4,2,4,3,1,4,2,3,4,2,3,4,2,1,2,4,4,3,1,3,4,2,3,3,2,4,3,2,2,2,0,2,3       ,
2,1,4,2,4,4,2,2,2,3,3,1,4,2,3,2,1,3,1,2,3,1,4,4,2,4,2,3,2,4,4,1,4,3,4,2,2,4,2,2,2,4,1,1,4,2,1,3,4,0,3,2,4        ,
4,2,4,3,2,4,4,2,2,0,1,1,2,3,2,1,2,2,1,4,NA,4,3,2,2,2,2,NA,2,4,3,3,2,2,3,4,2,2,2,3,2,1,2,3,3,3,2,3,3,3,3,3,1      ,
3,2,1,3,2,3,2,3,3,2,3,4,2,3,1,2,3,3,4,3,3,3,1,1,2,4,1,2,1,3,4,4,4,2,2,3,2,4,3,1,3,3,3,1,1,2,3,2,3,2,4,4,0        ,
2,4,2,3,4,3,1,4,4,0,2,4,0,2,4,4,2,2,4,2,2,2,4,2,NA,4,2,1,2,4,2,2,1,1,4,3,3,1,3,4,3,2,3,1,2,3,0,4,3,4,2,4,3       ,
2,2,1,1,2,2,3,1,3,2,1,2,2,3,2,2,3,1,4,2,3,3,3,2,4,3,1,4,0,1,2,2,2,1,2,2,3,2,4,2,1,2,1,2,2,2,1,1,2,2,2,4,3        ,
0,2,4,NA,2,2,0,1,2,1,1,4,2,4,3,2,2,2,4,4,2,2,3,2,1,1,4,1,3,1,1,4,2,3,1,2,2,3,4,1,1,3,1,3,NA,4,3,1,2,2,2,2,3      ,
3,4,3,1,2,3,1,3,3,2,0,3,1,3,0,4,0,3,4,4,3,3,3,3,4,1,2,2,2,3,4,1,4,3,1,4,2,0,0,1,2,2,4,1,4,2,2,3,4,3,3,3,2        ,
2,4,2,3,4,2,3,1,1,1,2,2,3,1,4,1,3,2,2,3,2,2,2,1,4,4,1,2,0,3,2,4,2,3,1,1,3,2,3,1,3,0,3,1,2,2,2,2,4)




     table(aid, white)

  byRace <- table(aid, white)
margin.table(byRace, 1) # A frequencies (summed over B)
margin.table(byRace, 2) # B frequencies (summed over A)
singleDV <-table(aidtoblacks)
prop.table(singleDV)


  #  0     1     2     3     4
#   0.032 0.160 0.357 0.253 0.199
#So, we can see that in our data set, knowing nothing about
#a particular respondent, we'd guess they'd fall into category
#0 with probability 0.032, and category 4 with probability 0.199
#Knowing this, let's look at the breakdown by race:

  ftable(byRace)
prop.table(byRace) # cell percentages

#We see here the respective probabilities for each racial
#group are given using this command:
prop.table(byRace, 2) # column percentages
#So, knowing only a respodnent's race, we'd estimate that
#whites fall into category 0 3% of the time, blacks 0% of
#the time, and category 4 with probability 0.15 and 0.42,
#respectively for whites and blacks. So, we can estimate
#a simple bi-variate logit model to check this.
#because our DV has five levels, we'd expect to see 4
#cutpoints estimated (Taus).

o.probit <- polr(as.ordered(aid)~ W, method=c("probit"))
summary(o.probit)


names(o.probit)
length(o.probit$zeta)
o.probit$zeta  #Our Cutpoints
o.probit$coef #Our coef. for being white.
o.probit$zeta[1]



pnorm(0)   #Just to remind you what pnorm is doing.
P.blacks <- c()
P.blacks[1] <- pnorm(o.probit$zeta[1])
P.blacks
P.blacks[2] <-  pnorm(o.probit$zeta[2]) - pnorm(o.probit$zeta[1])
P.blacks
P.blacks[3] <-  pnorm(o.probit$zeta[3]) - pnorm(o.probit$zeta[2])
P.blacks
P.blacks[4] <-  pnorm(o.probit$zeta[4]) - pnorm(o.probit$zeta[3])
P.blacks
P.blacks[5] <-  1 -  pnorm(o.probit$zeta[4])
P.blacks
sum(P.blacks)
print( t(P.blacks))
#Compare this to our original results:
 prop.table(byRace, 2) # column percentages
 o.probit$zeta
 o.probit
white.zetas <- o.probit$zeta - o.probit$coef
#Subtract because Model is estimating -xB
white.zetas

P.whites <- c()
P.whites[1] <- pnorm(white.zetas[1])
P.whites
P.whites[2] <-  pnorm(white.zetas[2]) - pnorm(white.zetas[1])
P.whites
P.whites[3] <-  pnorm(white.zetas[3]) - pnorm(white.zetas[2])
P.whites
P.whites[4] <-  pnorm(white.zetas[4]) - pnorm(white.zetas[3])
P.whites
P.whites[5] <-  1 -  pnorm(white.zetas[4])
P.whites
sum(P.whites)

our.estimates <- cbind(P.blacks, P.whites)
our.estimates

#compare to original probabilities:
 prop.table(byRace, 2) # column percentages