############################################################################# ## ## PROGRAM TO IMPLEMENT RANDOM EFFECTS IDEAL POINT MODEL WITH COVARIATES ## Michael Bailey, Georgetown University Dept of Government ## March 1998; Updated September 2000 ## For complete list of references, see Bailey (2000). ## source("C:\\Mydocs\\Research\\REPAL\\BayesPA2000.q") ## ############################################################################## ##################################### ## ## Manually Set Parameters ## ##################################### H_200 ## Maximum number of iterations asig_0.25 ## Dispersion parameter for alpha distribution (!= variance) CONVERGE_0.00001 ## Convergence crit. (max difference in parameter estimates) GL_30 ## 10, 30 point Gauss-Legendre quadrature; ##################################### ## ## Read in and Organize Data ## ##################################### data_read.table("C:\\MyDocs\\Research\\REPAL\\Senate1993PAData.txt", header=T) avex_scale( data$Ex92 + data$Ex93 + data$Ex94 ) avic_scale( data$ImCh92 +data$ImCh93+data$ImCh94 ) Skill_scale((data$College +data$Grad)/ (data$NoHS+ data$HS+data$SomeColl+data$AssDeg+ data$College+data$Grad)) xcov_cbind(Skill, avex, scale(data$MNC9094PCT), avic, scale(data$union93), data$dem*scale(data$union93), scale(data$Labor9094PCT), scale(data$Perot92), data$dem, rep(1, N)) cov.names_c(" Skill", "ex", " MNC9094PCT", " imch", " Union", " Dem*Union", " Labor9094PCT", " Perot", "dem", " const") yobs_ cbind(data$Nafta, data$NSA, data$GATTFT, data$Gatt94, data$GattBW ) yobs_replace(yobs, yobs==9, NA) vote.names_ c("Nafta", "NSA", " GattFT", " Gatt94", " GattBW" ) N_length(data$dem) ## Number of Senators TT_length(yobs[1,]) ## Number of votes P_length(xcov[1,]) ## Number of covariates if(length(b)!= P) {b_rep(1, P)} ## Starting values for b if(length(a)!=TT) a_rep(1, TT) ## Starting values for a if(length(ak)!=TT) ak_rep(0, TT) ## Starting values for ak k_ak/a ## Starting values for k ximinus1_xi_0 ## For use in convergence check ####################################### ## ## Initialize variables and functions ## ####################################### ## Select quadrature points if(GL == 10){ ## Gauss-Legendre weights for 10 point (Judd, 260) w_c(0.066713443, .1494513491, .210863625, 0.2692667193, 0.2955242247, 0.2955242247, 0.2692667193, .210863625, .1494513491, 0.066713443 ) xiGL_c(-0.9739065285, -0.8650633667, -0.6794095682, -0.4333953941, - 0.1488743389, 0.1488743389, 0.4333953941, 0.6794095682, 0.8650633667, 0.9739065285) thetaq_-5+5*(xiGL+1) ## thetaq with max possible = 5 } # End GL==10 if(GL == 30){ ## Gauss-Legendre weights for 30 point (Stroud and Secrest ) xiGL_c(-0.996893484, -0.9836651232, -0.9600218649, -0.9262000474, -0.8825605357, -0.8295657623, -0.7677774321, -0.6978504947, -0.6205261829, -0.5366241481, -0.44703376695, -0.3527047255, -0.2546369261, -0.1538699136, -0.05147184255, 0.051471843, 0.153869914, 0.254636926, 0.352704726, 0.447033767, 0.536624148, 0.620526183, 0.697850495, 0.767777432, 0.829565762, 0.882560536, 0.926200047, 0.960021865, 0.983665123, 0.996893484) w_c(0.00796819, 0.01846646, 0.028784, 0.038799, 0.048402, 0.0574931, 0.06597422, 0.07375559, 0.08075589, 0.086899, 0.092122522, 0.096368, 0.09959342, 0.1017623, 0.1028526, 0.1028526, 0.1017623, 0.09959342, 0.096368, 0.092122522, 0.086899, 0.08075589, 0.07375559, 0.06597422, 0.0574931, 0.048402, 0.038799, 0.028784, 0.01846646, 0.00796819) thetaq_-5+5*(xiGL+1) ## thetaq with max possible = 5 } # End GL==30 Q_length(thetaq) ## Number of quadrature nodes ## Useful variables onesTT_rep(1, TT) onesN_rep(1, N) ## Initialize matrices Ptq_matrix(NA, TT, Q) PTtheta_array(NA, dim=c(N, TT, Q)) ## Prob(y_it|q, a, k) Liq_matrix(NA, N, Q) ## L_i(X_q) - prob Y_i|theta=theta_q ritq_array(0, dim=c(N, TT, Q)) ## yti * L*W/sum of L*W over q fitq_array(0, dim=c(N, TT, Q)) ## P(i is thetaq) for yobs!=NA ## Create function used in M-Step ## NOTE: Splus minimizes maxfun; so set to negative to get max maxfun_function(xi, q, r, f, w){ log(xi[1]) + (.5/asig^2)*(log(xi[1])-asig^2)^2 + sum(w*( r*log(1 + exp(xi[2] - xi[1]*q)) + (f-r)*log(1+exp(xi[1]*q - xi[2])) ) ) } ## Mode of a lognormal r.v. is exp(mu - sigma^2) ## --> for mode = 1 ==>, mu = sigma^2 ## Alpha is dist'ed LN(1, var) where mu = -.5, sigma = asig (above) rm(thetai) ## Delete thetai with names attached from previous estimation ####################################### ## ## Iterations ## ####################################### for(h in 1:H){ print (c("Iteration Number ", h)) ####################################### ## ## Expectation Step ## ####################################### for(q in 1:Q) { ## Start first q loop Ptq[,q]_(1+exp(-a*thetaq[q] +ak))^(-1) ## TTx1 matrix for each q PTtheta[,,q]_(outer(rep(1, N), Ptq[,q])^yobs)*( outer(rep(1, N), (1-Ptq[,q]))^(1-yobs)) ## Prob y_it given q and vt parameters PTtheta[,,q]_replace(PTtheta[,,q], PTtheta[,,q] =="NA", 1) ## Cumulative product, so ok to replace with 1's Liq[, q]_exp(log(PTtheta[,,q])%*%rep(1, TT)) ## prob of i's t votes given q } ## End first q loop ProbThetaq_(1/sqrt(2*pi))*exp((-.5)*(outer(rep(1, N), thetaq)- crossprod(t(xcov), b)%*%rep(1, Q))^2) ## Prob Thetaq given X'beta (numerator (n x Q) Ptheta_(Liq*ProbThetaq)/(((Liq*ProbThetaq)%*%w)%*%t(rep(1, Q))) for (t in 1:TT){ ## Start first t loop ritq[,t,]_outer(yobs[,t],rep(1, Q))*Ptheta ## prob(i is q)*y_it fitq[,t,]_Ptheta*((1+outer(yobs[,t],rep(1, Q)))/(1+outer(yobs[,t],rep(1, Q)))) ## jerry-rigged way to del. missing obs from fitq matrix (which is ## length(beta)(i is q)|y_it!=NA) } ## End first t loop ritq_replace(ritq, ritq=="NA", 0) ## Am summing, so ok to replace w/ fitq_replace(fitq, fitq=="NA", 0) ## zeroes rtq_apply(ritq, c(2, 3), sum) ## Sum y_it*E(type q) over i for each t ftq_apply(fitq, c(2, 3), sum) ## Sum of fiqt over i's for each t & q thetai_Ptheta%*%(thetaq*w) ## E(theta_i) - no w b/c Ptheta sums to 1 thetai_replace(thetai, Liq[,1]==1, NA) ## Liq[,1]=1 is for gthetai_thetai[!is.na(thetai)] gid_data$id[1:N][!is.na(thetai)] gxcov_matrix(xcov[!is.na(thetai)], ,length(b)) b_solve(crossprod(gxcov)) %*%crossprod(gxcov, gthetai) print(b) ####################################### ## ## Maximization Step ## ####################################### ## IDENTIFY WITH A1_1, K1_0 a[1]_1 ak[1]_0 for(t in 2:TT){ max_nlminb(start=c(a[t], ak[t]), objective = maxfun, lower=c(0, -Inf), r=rtq[t,], q=thetaq, f=ftq[t,], w=w) ak[t]_max$parameters[2] a[t]_ max$parameters[1] } ## End for t in 2:TT k_ak/a print(cbind(a, k)) ####################################### ## ## Convergence Loop ## ####################################### ximinus1_xi ## For convergence test xi_c(a, k, b) ## For convergence test print(c(" convergence criterion ", max(abs((xi-ximinus1))) )) if((max(abs((xi-ximinus1))))