#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