#R script for "Regression Models for Binary Dependent Variables"
#				 http://www.indiana.edu/~statmath/stat/all/cdvm/index.html
#
#Created on 09/17/2009
#Last Modified on September 25, 2009
#
#Written by Hun Myoung Park, Ph.D. 
#	Statistical Software Analyst and Programmer
#	UITS Center for Statistical and Mathematical Computing, Indiana University
#
#University Information Technology Services, Indiana University
#Center for Statistical and Mathematical Computing
#http://www.indiana.edu/~statmath
#
#**************************************************************************************/


df<-read.table('http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.csv', sep=',', header=T)

attach(df)
summary(df)


# Binary Logit Model

blm<-glm(trust~educate+income+age+male+www, data=df, family=binomial(link="logit"))
summary(blm)

blm$deviance/-2
AIC(blm)

LRtest<-blm$null.deviance-blm$deviance
LRtest
LRp<-dchisq(LRtest, blm$df.null - blm$df.residual)
1-bpm$deviance/bpm$null.deviance # McFadden's pseudo R square


meanX<-c(1, mean(educate), mean(income), mean(age), mean(male), mean(www))
sdX<-c(1, sd(educate), sd(income), sd(age), sd(male), sd(www))
bHat<-coef(blm) 		# vector of parameter estimates
K<-length(bHat) 		# the number of parameters

fcOdds<-cbind(bHat, exp(bHat), exp(bHat*sdX), meanX, sdX)
fcOdds<-fcOdds[2:K,]
colnames(fcOdds)<-c("b", "exp(b)", "exp(b*sd)", "Mean of X", "SD of X")
fcOdds

referX<-c(1, 16, mean(income), mean(age), 0, 1) # set reference points
xb<-bHat %*% referX # element by element product
prob<-exp(xb)/(1+exp(xb)) # compute a pridicted probability
prob

margEffect<-cbind(bHat, prob*(1-prob)*bHat, prob*(1-prob)*bHat*sdX, meanX, sdX)
margEffect<-margEffect[2:K,]
colnames(margEffect)<-c("b", "MargEffect", "MargEffect(SD)", "Mean of X", "SD of X")
margEffect






# Binary Probit Model

bpm<-glm(trust~educate+income+age+male+www, data=df, family=binomial(link="probit"))
summary(bpm)

bpm$deviance/-2
AIC(bpm)

LRtest<-bpm$null.deviance-blm$deviance
LRtest

dchisq(LRtest, bpm$df.null - bpm$df.residual)
1-bpm$deviance/bpm$null.deviance # McFadden's pseudo R square

meanX<-c(1, mean(educate), mean(income), mean(age), mean(male), mean(www))
sdX<-c(1, sd(educate), sd(income), sd(age), sd(male), sd(www))
bHat<-coef(bpm) 		# vector of parameter estimates
K<-length(bHat) 		# the number of parameters


referX<-c(1, 16, mean(income), mean(age), 0, 1)
xb<-bHat %*% referX # element by element product
prob<-pnorm(xb)
prob

margEffect<-cbind(bHat, dnorm(xb)*bHat, dnorm(xb)*bHat*sdX, meanX, sdX)

for (i in c(5, 6)) {
    referX0<-matrix(referX)
    referX1<-matrix(referX)
    referX0[i,1]<-0
    referX1[i,1]<-1

    xb0<-bHat %*% referX0
    xb1<-bHat %*% referX1

    dChange<-pnorm(xb1)-pnorm(xb0)
    margEffect[i,2]<-dChange # replace the marginal effect with the discrete change
}

margEffect<-margEffect[2:K,]
colnames(margEffect)<-c("b", "MargEffect", "MargEffect(SD)", "Mean of X", "SD of X")
margEffect

# End of the file