# R FILE NUMBER 3 FOR "THE POLITICAL ENTROPY OF VOTE CHOICE", WINTER 2004 # THIS FILE BUILDS THE STATISTICAL MODEL FOR THE ENTROPY ESTIMATION # ALL ANALYSIS ASSUMES THAT THE FILE entropy.data.s AND entropy.prior.s # HAS BEEN RUN TO CONDITION THE DATA AND SET THE ENTROPY TERMS # HERE'S A REMINDER OF THE VARIABLE NAMES # [1] "Apr.Clt" "Dem.Supp.Clt" "Rep.Supp.Clt" "Inc.Spt.Clt.alw" # [5] "Inc.Spt.Clt.nvr" "Inc.Party" "Supp.CB" "Dem.Vote.CB" # [9] "Rep.Vote.CB" "Gvt.Help" "Dem.Gvt.Help" "Rep.Gvt.Help" #[13] "Gvt.Spend" "Dem.Gvt.Spend" "Rep.Gvt.Spend" "Fed.Health" #[17] "Dem.Fed.Health" "Rep.Fed.Health" "Remember.Names" "Vote" #[21] "Party.Vote" "Register.Vote" "Party.ID" "NV.Preference" #[25] "Ideology" "Dem.Ideology" "Dem.Certainty" "Rep.Ideology" #[29] "Rep.Certainty" "Age" "Education" "Gender" #[33] "Race" "Pol.Interest" "Read.P" "Media.Exp" print("################## SETUP NEW SUBSET STRUCTURE: e.model ################################") e.model <- matrix(0,nrow=nrow(e.data2),ncol=18) dimnames(e.model) <- list(NULL,c("Dem.Supp.Clt", "Rep.Supp.Clt", "Dem.Vote.CB", "Rep.Vote.CB", "Dem.Gvt.Help", "Rep.Gvt.Help", "Dem.Gvt.Spend", "Rep.Gvt.Spend", "Dem.Fed.Health", "Rep.Fed.Health", "Party.ID", "Pol.Interest", "Education", "Gender","Race", "Media.Exp", "Age", "Party.Vote")) # DEMOCRAT SQUARED DISTANCE, SUPPORT FOR CLINTON + ENTROPY e.model[,"Dem.Supp.Clt"] <- (e.data2[,"Apr.Clt"] - e.data2[,"Dem.Supp.Clt"])^2 + entropy.Dem.Supp.Clt # REPUBLICAN SQUARED DISTANCE, SUPPORT FOR CLINTON + ENTROPY e.model[,"Rep.Supp.Clt"] <- (e.data2[,"Apr.Clt"] - e.data2[,"Rep.Supp.Clt"])^2 + entropy.Rep.Supp.Clt # DEMOCRAT SQUARED DISTANCE, CRIME BILL + ENTROPY e.model[,"Dem.Vote.CB"] <- (e.data2[,"Supp.CB"] - e.data2[,"Dem.Vote.CB"])^2 + entropy.Dem.Supp.CB # REPUBLICAN SQUARED DISTANCE, CRIME BILL + ENTROPY e.model[,"Rep.Vote.CB"] <- (e.data2[,"Supp.CB"] - e.data2[,"Rep.Vote.CB"])^2 + entropy.Rep.Supp.CB # DEMOCRAT SQUARED DISTANCE, JOB CREATION + ENTROPY e.model[,"Dem.Gvt.Help"] <- (e.data2[,"Gvt.Help"] - e.data2[,"Dem.Gvt.Help"])^2 + entropy.Dem.Gvt.Help # REPUBLICAN SQUARED DISTANCE, JOB CREATION + ENTROPY e.model[,"Rep.Gvt.Help"] <- (e.data2[,"Gvt.Help"] - e.data2[,"Rep.Gvt.Help"])^2 + entropy.Rep.Gvt.Help # DEMOCRAT SQUARED DISTANCE, GOVERNMENT SPENDING + ENTROPY e.model[,"Dem.Gvt.Spend"] <- (e.data2[,"Gvt.Spend"] - e.data2[,"Dem.Gvt.Spend"])^2 + entropy.Dem.Gvt.Spend # REPUBLICAN SQUARED DISTANCE, GOVERNMENT SPENDING + ENTROPY e.model[,"Rep.Gvt.Spend"] <- (e.data2[,"Gvt.Spend"] - e.data2[,"Rep.Gvt.Spend"])^2 + entropy.Rep.Gvt.Spend # DEMOCRAT SQUARED DISTANCE, GOVERNMENT SPENDING + ENTROPY e.model[,"Dem.Fed.Health"] <- (e.data2[,"Fed.Health"]-e.data2[,"Dem.Fed.Health"])^2 + entropy.Dem.Fed.Health # REPUBLICAN SQUARED DISTANCE, GOVERNMENT SPENDING + ENTROPY e.model[,"Rep.Fed.Health"] <- (e.data2[,"Fed.Health"]-e.data2[,"Rep.Fed.Health"])^2 + entropy.Rep.Fed.Health e.model <- cbind(e.model,matrix(0,nrow=nrow(e.model),2)) nc <- ncol(e.model) dimnames(e.model)[[2]][(nc-1):(nc)] <- c("Dem.NP.Ideology","Rep.NP.Ideology") # DEMOCRAT SQUARED DISTANCE, NO PRIOR IDEOLOGY e.model[,"Dem.NP.Ideology"] <- (e.data2[,"NP.Ideology"]-e.data2[,"NP.Dem.Ideology"])^2 + entropy.Dem.NP.Ideology # REPUBLICAN SQUARED DISTANCE, NO PRIOR IDEOLOGY e.model[,"Rep.NP.Ideology"] <- (e.data2[,"NP.Ideology"]-e.data2[,"NP.Rep.Ideology"])^2 + entropy.Rep.NP.Ideology # PARTY ID e.model[,"Party.ID"] <- e.data2[,"Party.ID"] # POLITICAL INTEREST e.model[,"Pol.Interest"] <- e.data2[,"Pol.Interest"] # EDUCATION IN YEARS, MEAN CENTERING ADDED TO THIS VARIABLE e.model[,"Education"] <- e.data2[,"Education"] - mean(e.data2[,"Education"]) # GENDER, MEAN CENTERING ADDED TO THIS VARIABLE e.model[,"Gender"] <- e.data2[,"Gender"] - mean(e.data2[,"Gender"]) # RACE, MEAN CENTERING ADDED TO THIS VARIABLE e.model[,"Race"] <- e.data2[,"Race"] - mean(e.data2[,"Race"]) # MEDIA EXPOSURE e.model[,"Media.Exp"] <- e.data2[,"Media.Exp"] # AGE, MEAN CENTERING ADDED TO THIS VARIABLE e.model[,"Age"] <- e.data2[,"Age"] - mean(e.data2[,"Age"]) # HOUSE VOTE BY PARTY: THE DEPENDENT VARIABLE e.model[,"Party.Vote"] <- e.data2[,"Party.Vote"] - 1 entropy.X <- cbind( rep(1,nrow(e.model)), e.model[,1:17] ) entropy.Y <- e.model[,"Party.Vote"] ######################## CODE THAT FOLLOWS IS FOR THE NON-HARVEY MODEL FIRST ########################## ################################################################################### # SPLUS:e.probit.out<-nlminb(start=e.like.start[1:14],objective=entropy.log.like, # # X=entropy.X,Y=entropy.Y) # ################################################################################### entropy.glm.probit.baseline <- glm(entropy.Y ~ entropy.X[,-1], family=binomial(probit)) summary(entropy.glm.probit.baseline) 18 - entropy.glm.probit.baseline$aic/2 # -641.3081 # THE LOG-LIKELIHOOD FUNCTION entropy.log.like <- function(beta,X,Y) { - sum( log(1-pnorm(X%*%beta))*(1-Y) ) - sum( log(pnorm(X%*%beta))*Y ) } # THE SECOND DERIVATIVE FUNCTION TO PRODUCE THE HESSIAN dd.log.like <- function(beta,X,Y) { lambda.0 <- (1-Y)*(-1*dnorm(X%*%beta)/(1-pnorm(X%*%beta))) lambda.1 <- Y*dnorm(X%*%beta)/(pnorm(X%*%beta)) ( -sum( lambda.0*(lambda.0 + X%*%beta)*(1-Y) )*t(X)%*%X - sum( lambda.1*(lambda.1 + X%*%beta)*(Y) )*t(X)%*%X ) } e.max.hist <- matrix(c(entropy.glm.probit$coefficients[1:12], entropy.glm.probit$coefficients[13:18]), nrow=1) for (i in 1:10) { e.max <- optim(fn=entropy.log.like,par=e.max.hist[i,],X=entropy.X[,1:18],Y=entropy.Y,hessian=TRUE) #e.max <- optim(fn=entropy.log.like,par=e.max.hist[i,], # X=cbind(entropy.X[,1:11],entropy.X[,13:14]),Y=entropy.Y,hessian=TRUE) print(e.max.hist[i,]) e.max.hist <- rbind(e.max.hist,e.max$par) print(paste("finished with iteration:",i,"loglike:",e.max$value,"convergence:",e.max$convergence)) } # [1] -1.1023385929 -0.0049456009 0.0180914701 0.0055367172 -0.0140183806 # [6] 0.0223805008 -0.0134281622 0.0781199069 -0.0799090317 0.0070205949 # [11] -0.0080678415 0.3690521905 -0.1621437145 0.0295289781 -0.3060869445 # [16] -0.3823269112 -0.0193017627 -0.0005972723 # [1] "finished with iteration: 10 loglike: 641.780627894882 convergence: 1" # CALCULATE THE STANDARD ERRORS FROM THE ASYMPTOTIC COVARIANCE # MATRIX (BY THE BE NEGATIVE INVERSE OF THE HESSIAN). entropy.hessian <- dd.log.like(e.max$par,entropy.X, entropy.Y) asym.cov.mat <- -solve(entropy.hessian)*nrow(entropy.X) sqrt(diag(asym.cov.mat)) # RUN A FULL GLM MODEL TO GET STARTING POINTS entropy.glm.probit.full <- glm(entropy.Y ~ entropy.X[,-1],family=binomial(probit)) summary(entropy.glm.probit.full) # Deviance Residuals: # Min 1Q Median 3Q Max # -2.4582 -0.7435 -0.3278 0.7214 2.7792 # # Coefficients: # Estimate Std. Error z value Pr(>|z|) # (Intercept) -1.3788779 0.2879028 -4.789 1.67e-06 *** # entropy.X[, -1]Dem.Supp.Clt -0.0046486 0.0069173 -0.672 0.501566 # entropy.X[, -1]Rep.Supp.Clt 0.0185810 0.0091050 2.041 0.041276 * # entropy.X[, -1]Dem.Vote.CB 0.0065277 0.0072550 0.900 0.368251 # entropy.X[, -1]Rep.Vote.CB -0.0119882 0.0071544 -1.676 0.093812 . # entropy.X[, -1]Dem.Gvt.Help 0.0231766 0.0090206 2.569 0.010191 * # entropy.X[, -1]Rep.Gvt.Help -0.0128598 0.0106864 -1.203 0.228833 # entropy.X[, -1]Dem.Gvt.Spend 0.0737514 0.0207614 3.552 0.000382 *** # entropy.X[, -1]Rep.Gvt.Spend -0.0732918 0.0206997 -3.541 0.000399 *** # entropy.X[, -1]Dem.Fed.Health 0.0075751 0.0071399 1.061 0.288713 # entropy.X[, -1]Rep.Fed.Health -0.0073288 0.0080043 -0.916 0.359872 # entropy.X[, -1]Party.ID 0.3733166 0.0228510 16.337 < 2e-16 *** # entropy.X[, -1]Pol.Interest -0.1212590 0.0608179 -1.994 0.046173 * # entropy.X[, -1]Education 0.0401078 0.0266079 1.507 0.131717 # entropy.X[, -1]Gender -0.3006725 0.0829970 -3.623 0.000292 *** # entropy.X[, -1]Race -0.3982267 0.1167557 -3.411 0.000648 *** # entropy.X[, -1]Media.Exp -0.0156525 0.0174187 -0.899 0.368865 # entropy.X[, -1]Age 0.0008031 0.0024248 0.331 0.740488 # --- # Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 # # (Dispersion parameter for binomial family taken to be 1) # # Null deviance: 1839.7 on 1334 degrees of freedom # Residual deviance: 1282.6 on 1317 degrees of freedom # AIC: 1318.6 aic2ll <- function(aic,p) { p - 0.5*aic } aic2ll(entropy.glm.probit.full$aic,entropy.glm.probit.full$rank) # -641.308 ####################### NOW MAXIMIZE THE HARVEY LIKELIHOOD FUNCTION ################################### ################################################################################### # SPLUS: h.e.2.probit.out <- nlminb(start=e.like.start,objective=harvey.log.like, # # X=entropy.X,Y=entropy.Y) # ################################################################################### harvey.log.like <- function(beta, X, Y) { g <- cbind(X[,1:11],X[,14:15],X[,17:18])%*%c(beta[1:11],beta[14:15],beta[17:18])/ exp(cbind(X[,12:13],X[,16])%*%c(beta[12],beta[13],beta[16])) - sum(log(1 - pnorm(g)) * (1 - Y)) - sum(log(pnorm(g)) * Y) } e.max.hist <- matrix(entropy.glm.probit.full$coefficients,nrow=1) for (i in 1:100) { e.max <- optim(fn=harvey.log.like,par=e.max.hist[i,],X=entropy.X,Y=entropy.Y,hessian=TRUE) print(e.max.hist[i,]) e.max.hist <- rbind(e.max.hist,e.max$par) print(paste("finished with iteration:",i,"loglike:",e.max$value,"convergence:",e.max$convergence)) } ################# CALCULATE THE STANDARD ERRORS FROM THE ASYMPTOTIC COVARIANCE ###################### ################### VC MATRIX (BY THE BE NEGATIVE INVERSE OF THE HESSIAN). ########################## harvey.entropy.hessian <- dd.log.like(e.max.hist[101,],entropy.X, entropy.Y) h.asym.cov.mat <- -solve(harvey.entropy.hessian)*nrow(entropy.X) ############################## BUILD THE OUTPUT TABLE ############################################## h.entropy.mat <- cbind( e.max.hist[101,],diag(chol(h.asym.cov.mat)), e.max.hist[101,]/diag(chol(h.asym.cov.mat)) ) dimnames(h.entropy.mat)[[2]] <- c("Coefficient","Std.Error","t-stat") h.entropy.mat h.entropy.mat <- cbind( e.max.hist[101,],diag(chol(h.asym.cov.mat)), e.max.hist[101,]-1.96*diag(chol(h.asym.cov.mat)), e.max.hist[101,]+1.96*diag(chol(h.asym.cov.mat)) ) dimnames(h.entropy.mat)[[2]] <- c("Coefficient","Std.Error","95% CI Lower","95% CI Upper") round(h.entropy.mat,4) # CALCULATE THE LRT AND ITS CHISQUARE VALUE 1-pchisq(2*(724-641),3) # PERCENT CORRECTLY CLASSIFIED coefs <- e.max.hist[101,] g <- cbind(entropy.X[,1:11],entropy.X[,14:15],entropy.X[,17:18])%*%c(coefs[1:11],coefs[14:15],coefs[17:18])/ exp(cbind(entropy.X[,12:13],entropy.X[,16])%*%c(coefs[12],coefs[13],coefs[16])) predicts <- pnorm(g) table(cut(predicts,b=c(0,0.37,1)),entropy.Y) # DEVIANCE, inlier problem d.sum <- 0 for (i in 1:length(entropy.Y)) { print(i) d.sum <- d.sum + (entropy.Y[i] * log(entropy.Y[i]/predicts[i])) #+ #((1 - entropy.Y[i]) * log((1-entropy.Y[i])/(1-predicts[i]))) } ############################## CHECK BY INCLUDING GAMMA VARS IN BETA VECTOR ######################## harvey.log.like <- function(beta, X, Y) { g <- cbind(X[,1:18])%*%c(beta[1:18])/ exp(cbind(X[,12:13],X[,16])%*%c(beta[12],beta[13],beta[16])) - sum(log(1 - pnorm(g)) * (1 - Y)) - sum(log(pnorm(g)) * Y) } e.max.hist <- matrix(entropy.glm.probit.full$coefficients,nrow=1) for (i in 1:500) { e.max <- optim(fn=harvey.log.like,par=e.max.hist[i,],X=entropy.X,Y=entropy.Y,hessian=TRUE) print(e.max.hist[i,]) e.max.hist <- rbind(e.max.hist,e.max$par) print(paste("finished with iteration:",i,"loglike:",e.max$value,"convergence:",e.max$convergence)) } harvey.entropy.hessian <- dd.log.like(e.max.hist[501,],entropy.X, entropy.Y) h.asym.cov.mat <- -solve(harvey.entropy.hessian)*nrow(entropy.X) h.entropy.mat <- cbind( e.max.hist[501,],diag(chol(h.asym.cov.mat)), e.max.hist[501,]/diag(chol(h.asym.cov.mat)) ) dimnames(h.entropy.mat)[[2]] <- c("Coefficient","Std.Error","t-stat") h.entropy.mat