## ## Drew Linzer ## dlinzer@emory.edu ## August 5, 2008 ## ## BlaydesLinzer-WP08.R ## ## Statistical analyses in paper by Lisa Blaydes & Drew Linzer: ## "The Political Economy of Women's Support for Fundamentalist Islam" ## World Politics. 2008, 60(4). ## library(poLCA) # for latent class analysis library(Design) # needed to use ols() and robcov() for Efron robust standard errors library(car) # for recode function library(gdata) # for renaming vars in a data frame library(vcd) # for spine plots, other "visualizations of categorical data" library(lattice) # for wireframe 3d plot # Opinion items # ///////////// # WVS label Description # ====== =============== ================================================================= # A006 religimpt How important is religion in your life (1 very impt...4 not impt) # C001 jobscarcemen When jobs are scarce men should have more right to a job than women (1 agree, 2 disagree, 3 neither) # D019 womneedschild Woman has to have children (1 not necessary, 2=needs children [after recode]) # D059 menbetterpolit Men make better political leaders (1 agree strongly...4 disagree strongly) # D060 univboy University is more important for a boy (1 agree strongly...4 disagree strongly) # D067 womveil Importance of the traits in a woman: Wearing a veil in public places (1 very impt...5 not impt at all) # D070 womrelig Important for a woman to be religious (1 agree strongly...5 disagree strongly) # D076 polygamy More than one wife ok (1 agree strongly...5 disagree strongly) # D077 wifeobey Wife must obey (1 agree strongly...5 disagree strongly) # F063 godimpt God important in your life (1 not at all...10 very) [nearly unanimous] # F064 relcomft Religion gives comfort (1 no comfort from religion, 2 yes comfort from religion [after recode]) [nearly unanimous] # F066 prayperday Pray per day (1 every day...7 never) # F102 polgod Politician must believe in God (1 agree strongly...5 disagree strongly) # F103 religvote Religious leaders influence vote (1 strongly should not...5 strongly should) # F105 religgovt Religious leaders influence govt (1 strongly should not...5 strongly should) # F111 sharia Implement sharia law (1 agree strongly...5 disagree strongly) # # Demographic characteristics # /////////////////////////// # WVS label Description # ====== =============== ================================================================= # F025 religion Religious denomination # S003 country Country # X001 sex Sex -- 0 Male, 1 Female [after recode] # X003 age Age # X003R age Age recoded into six categories # X007 married Marital status # X011 children Have you had any children # X025 educ Highest educational level attained # X028 employed Are you employed now # X045 class Social class (subjective) # X047 income Income scale; 10 category low to high # S017 weight Weight ## Download data from http://www.worldvaluessurvey.org/ ## Notes on data. ## [1] This analysis uses the "v4" WVS data file, "WVSEVS_sb_v4.SAV". ## [2] Because the data file is so large, we created a subset of the original file ## containing only the 18 countries under investigation, and only the questions ## needed for analysis. We then converted this file to R format as wvs_muslim_18.put. setwd("C:/.../") # replace /.../ with working directory wvs.file <- dget("wvs_muslim_18.put") wvs <- wvs.file ## Recode variables as needed wvs$D019 <- recode(wvs$D019,"0=1;1=2") wvs[wvs$S003==12,]$F025 <- 49 # recode Algeria where religion q. wasn't asked wvs[wvs$S003==586 & wvs$F025==996,]$F025 <- 49 # recode Pakistan which asked religion q. differently wvs$F025 <- recode(wvs$F025,"c(2,49,70,75)=1;c(53,996,NA)=3;else=2") # 1=Muslim, 2=non-Muslim, 3=no religion wvs$F064 <- recode(wvs$F064,"0=1;1=2") wvs$X001 <- recode(wvs$X001,"1=0;2=1") # 0=Male, 1=Female wvs$X007 <- recode(wvs$X007,"c(1,3,4,5)=1;c(2,6)=0") # 1=Ever married, 0=Never married wvs$X045 <- recode(wvs$X045,"1=5;2=4;4=2;5=1") # Social class, low to high wvs$S017 <- NULL wvs <- rename.vars(wvs,c("S003", "F025", "X001","X003","X007", "X025","X028", "X045", "X047" ), c("country","religion","sex", "age", "married","educ","employed","class","income")) cnames <- c("Algeria","Azerbaijan","Bangladesh","BosniaHerz","India","Indonesia","Iran","Jordan","Morocco", "Nigeria","Pakistan","Singapore","Turkey","Uganda","Macedonia","Egypt","Tanzania","Montenegro") cnum <- as.numeric(names(table(wvs$country))) # c(12,31,50,70,356,360,364,400,504,566,586,702,792,800,807,818,834,912) ## create dichotomous umemployment variable (housewife or unemployed) wvs$unempl <- ((wvs$employed==5) | (wvs$employed==7)) + 0 ## remove respondents with no sex reported wvs <- wvs[!is.na(wvs$sex),] ## ID COUNTRY ## === ========== ## 12 Algeria ## 31 Azerbaijan ## 50 Bangladesh ## 70 Bosnia and Herzegovina ## 356 India ## 360 Indonesia ## 364 Iran ## 400 Jordan ## 504 Morocco ## 566 Nigeria ## 586 Pakistan ## 702 Singapore ## 792 Turkey ## 800 Uganda ## 807 Macedonia ## 818 Egypt ## 834 Tanzania ## 912 Montenegro ## ## Appendix A. Countries in study, and total number of individuals interviewed ## in each country by the WVS, tabulated by sex and Muslim religion. ## tabA <- cbind(Total =table(wvs$country), MuslimTot=table(wvs$country[wvs$religion==1]), MuslimM =table(wvs$country[wvs$religion==1 & wvs$sex==0]), MuslimF =table(wvs$country[wvs$religion==1 & wvs$sex==1]), PctMuslim=100*round(table(wvs$country[wvs$religion==1])/table(wvs$country),3), PctSamp =100*round(table(wvs$country[wvs$religion==1])/nrow(wvs[wvs$religion==1,]),3), PctWorldPop=100*c(0.0209,0.0048,0.0842,0.0015,0.1135,0.1386,0.0433,0.0036,0.0212, 0.0418,0.1044,0.0005,0.0448,0.0027,0.0004,0.0473,0.0107,0.0001)) tabA <- as.data.frame(tabA) tabA <- tabA[order(cnames),] tabA <- rbind(tabA,colSums(tabA)) tabA$PctMuslim[nrow(tabA)] <- 100*round(tabA$MuslimTot[nrow(tabA)]/tabA$Total[nrow(tabA)],3) rownames(tabA) <- c(sort(cnames),"Total") cat("Table (Appendix A): Countries in study, sample sizes. \n") print(tabA) ## ## Remove non-Muslims from data set. ## ## Note that 2541 respondents report no religion -- either "other" (96), ## "no religion" (1052), or "NA" (1393). Below, we check the possibility ## that secular Muslims are failing to identify themselves as Muslim since ## they are not necessarily observant of the religion. ## wvs.norelig <- wvs[wvs$religion==3,] wvs <- wvs[wvs$religion==1,] ## ## Most women are both employed and married. ## wvs.f <- wvs[wvs$sex==1,] tab.empl.marr <- table(wvs.f$unempl,wvs.f$married) rownames(tab.empl.marr) <- c("empl","hw/un") colnames(tab.empl.marr) <- c("single","married") tab.empl.marr <- round(tab.empl.marr/sum(tab.empl.marr),4)*100 cat("Women's employment and marital status. \n") print(tab.empl.marr) # table(wvs.f$unempl,wvs.f$married,wvs.f$country) ## ## First, decide how many latent classes to fit. ## ## The goal is to have enough classes to be able to identify and isolate the ## "fundamentalist" grouping: people who hold a traditionalist view on ## gender roles, and are also religious and believe that religion should be ## play a role in political affairs. ## ## To do this, gradually increment the number of latent classes until adding ## an additional class does not further narrow the classification of respondents ## with fundamentalist belief systems. ## f <- cbind(A006,C001,D019,D059,D060,D067,D070,D076,D077,F063,F064,F066,F102,F103,F105,F111)~1 lcfit2 <- poLCA(f,wvs,na.rm=F,nclass=2) table(lcfit2$posterior[,1]>0.75 | lcfit2$posterior[,1]<0.25)/nrow(wvs) table(lcfit2$posterior[,2]>0.75 | lcfit2$posterior[,2]<0.25)/nrow(wvs) # log-likelihood: -265671.8 # Estimated class population shares (percent w/ posterior prob. > 0.75) # secular 0.255 (88.5%) # fundamentalist 0.745 (88.5%) lcfit3 <- poLCA(f,wvs,na.rm=F,nclass=3) table(lcfit3$posterior[,1]>0.75 | lcfit3$posterior[,1]<0.25)/nrow(wvs) table(lcfit3$posterior[,2]>0.75 | lcfit3$posterior[,2]<0.25)/nrow(wvs) table(lcfit3$posterior[,3]>0.75 | lcfit3$posterior[,3]<0.25)/nrow(wvs) # log-likelihood: -261582.3 # Estimated class population shares (percent w/ posterior prob. > 0.75) # secular 0.142 (93.7%) # conservative 0.394 (76.0%) # fundamentalist 0.464 (79.3%) lcfit4 <- poLCA(f,wvs,na.rm=F,nclass=4,nrep=5) table(lcfit4$posterior[,1]>0.75 | lcfit4$posterior[,1]<0.25)/nrow(wvs) table(lcfit4$posterior[,2]>0.75 | lcfit4$posterior[,2]<0.25)/nrow(wvs) table(lcfit4$posterior[,3]>0.75 | lcfit4$posterior[,3]<0.25)/nrow(wvs) table(lcfit4$posterior[,4]>0.75 | lcfit4$posterior[,4]<0.25)/nrow(wvs) # log-likelihood: -259133.9 # Estimated class population shares (percent w/ posterior prob. > 0.75) # secular 0.120 (94.9%) # religious 0.246 (81.8%) # traditionalist 0.298 (78.8%) # fundamentalist 0.336 (82.1%) lcfit5 <- poLCA(f,wvs,na.rm=F,nclass=5,tol=1e-4,nrep=5) table(lcfit5$posterior[,1]>0.75 | lcfit5$posterior[,1]<0.25)/nrow(wvs) table(lcfit5$posterior[,2]>0.75 | lcfit5$posterior[,2]<0.25)/nrow(wvs) table(lcfit5$posterior[,3]>0.75 | lcfit5$posterior[,3]<0.25)/nrow(wvs) table(lcfit5$posterior[,4]>0.75 | lcfit5$posterior[,4]<0.25)/nrow(wvs) table(lcfit5$posterior[,5]>0.75 | lcfit5$posterior[,5]<0.25)/nrow(wvs) # log-likelihood: -257910.3 # Estimated class population shares (percent w/ posterior prob. > 0.75) # secular 0.048 (98.0%) # moderate 0.170 (86.9%) # conservative 0.253 (80.1%) # conservative 0.235 (81.9%) # fundamentalist 0.294 (83.5%) ## ## There is a small, identifiable secular core of about 5% in the five- ## class model. In the four-class model, the secular group comprises 12%. ## The question, though, is how to partition the non-secular respondents. ## For this purpose, the four-class model provides a more useful fit; the ## five-class model divides up the traditionalist-religious respondents ## from the four-class model into three more indistinguishable groups. Most ## importantly, the individuals identified as fundamentalist in the four- ## class model are nearly idendical to those identified as fundamentalist ## in the five-class model. Even when a sixth class is added, the size ## of the fundamentalist group is still about 30%, and comprises more ## or less the same survey respondents. ## # table(lcfit4$predclass,lcfit5$predclass) ## ## Next, add covariates to the four-class model to see what variables ## predict membership in each of the latent classes. ## f.cov <- cbind(A006,C001,D019,D059,D060,D067,D070,D076,D077,F063,F064,F066,F102,F103,F105,F111)~sex*educ+sex*I(educ^2)+sex*class+sex*unempl+sex*married probs.start <- dget("LCmodel-start.put") lcfit.cov <- poLCA(f.cov,wvs,na.rm=F,nclass=4,probs.start=probs.start) # llik -230596.8 coeff.m <- lcfit.cov$coeff[c(1,3:7),] coeff.f <- lcfit.cov$coeff[c(1,3:7),] + lcfit.cov$coeff[c(2,8:12),] coeff.se.m <- lcfit.cov$coeff.se[c(1,3:7),] coeff.se.f <- sqrt( (lcfit.cov$coeff.se[c(1,3:7),]^2) + (lcfit.cov$coeff.se[c(2,8:12),]^2) + (2*matrix(lcfit.cov$coeff.V[matrix(c(2,8,9,10,11,12, 14,20,21,22,23,24, 26,32,33,34,35,36, 1,3,4, 5, 6, 7, 13,15,16,17,18,19, 25,27,28,29,30,31), nrow=18,ncol=2,byrow=F)],nrow=6,ncol=3,byrow=F)) ) # Latent Classes # 1 2 3 4 # Coefficients (secular) Fund. Trad. Relig. # 1 Intercept # 2 sex # 3 educ # 4 I(educ^2) # 5 class # 6 unempl # 7 married # 8 sex:educ # 9 sex:I(educ^2) # 10 sex:class # 11 sex:unempl # 12 sex:married # as reported in paper: print(coeff.m) print(coeff.se.m) print(coeff.f) print(coeff.se.f) ## ## Figure 1. Illustrating the effects of education, for middle class women ## who are unemployed and married, the group most fundamentalist. ## ## 1. Secular dash-dot ## 2. Fundamentalist solid ## 3. Traditional dash ## 4. Religious dot plot.pv <- function(xc,b,xaxs) { palette(rep("black",4)) par(mar=c(5,4,1,1)+0.1) exb <- exp(xc %*% b) matplot(xaxs$scale,(cbind(1,exb)/(1+rowSums(exb))),ylim=c(0,0.6), type="l",cex.axis=1.5,cex.lab=1.5,xaxt="n",lwd=4,lty=c(4,1,2,3), xlab=xaxs$lab,ylab="Probability of belief system type") print(cbind(1,exb)/(1+rowSums(exb))) axis(1,c(1:8),labels=c("< Elem.","Elementary","< Tech.","Technical","< Uni prep","Uni prep","Some Uni","Uni grad")) } # 1 1 1 # 1 2 3 4 5 6 7 8 9 0 1 2 xF.unem.marr <- cbind(1,1, c(1:8),c(1:8)^2,3,1,1, c(1:8),c(1:8)^2,3,1,1) windows(9,7) plot.pv(xF.unem.marr,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) text(7,0.51,"Religious",cex=1.5) text(5,0.28,"Traditional",cex=1.5) text(3.35,0.05,"Secular-Liberal",cex=1.5) text(2.3,0.175,"Fundamentalist",cex=1.5) #savePlot(filename="educ-unem-mc-predval",type="ps") #savePlot(filename="educ-unem-mc-predval",type="png") ## ## The effects of education, class, marital status, and employment for men. ## #xM.unem.marr.low <- cbind(1,0, c(1:8),c(1:8)^2,1,1,1, 0,0,0,0,0) #xM.empl.marr.low <- cbind(1,0, c(1:8),c(1:8)^2,1,0,1, 0,0,0,0,0) #xM.unem.marr.hi <- cbind(1,0, c(1:8),c(1:8)^2,4,1,1, 0,0,0,0,0) #xM.empl.marr.hi <- cbind(1,0, c(1:8),c(1:8)^2,4,0,1, 0,0,0,0,0) #xM.unem.unmr.low <- cbind(1,0, c(1:8),c(1:8)^2,1,1,0, 0,0,0,0,0) #xM.empl.unmr.low <- cbind(1,0, c(1:8),c(1:8)^2,1,0,0, 0,0,0,0,0) #xM.unem.unmr.hi <- cbind(1,0, c(1:8),c(1:8)^2,4,1,0, 0,0,0,0,0) #xM.empl.unmr.hi <- cbind(1,0, c(1:8),c(1:8)^2,4,0,0, 0,0,0,0,0) # #plot.pv(xM.unem.marr.low,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.empl.marr.low,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.unem.marr.hi,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.empl.marr.hi,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.unem.unmr.low,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.empl.unmr.low,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.unem.unmr.hi,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) #plot.pv(xM.empl.unmr.hi,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) ## ## Table 2. Effects of non-employment, marriage on the tendency of ## Muslim women to hold fundamentalist belief systems. ## # unempl 1 1 1 # 1 2 3 4 5 6 7 8 9 0 1 2 xF.unmar.elem.lowsc <- cbind(1,1, 2,4, 1,c(0:1),0, 2,4, 1,c(0:1),0) xF.unmar.elem.hisc <- cbind(1,1, 2,4, 4,c(0:1),0, 2,4, 4,c(0:1),0) xF.unmar.secd.lowsc <- cbind(1,1, 6,36,1,c(0:1),0, 6,36,1,c(0:1),0) xF.unmar.secd.hisc <- cbind(1,1, 6,36,4,c(0:1),0, 6,36,4,c(0:1),0) xF.mar.elem.lowsc <- cbind(1,1, 2,4, 1,c(0:1),1, 2,4, 1,c(0:1),1) xF.mar.elem.hisc <- cbind(1,1, 2,4, 4,c(0:1),1, 2,4, 4,c(0:1),1) xF.mar.secd.lowsc <- cbind(1,1, 6,36,1,c(0:1),1, 6,36,1,c(0:1),1) xF.mar.secd.hisc <- cbind(1,1, 6,36,4,c(0:1),1, 6,36,4,c(0:1),1) # marrid 1 1 1 # 1 2 3 4 5 6 7 8 9 0 1 2 xF.empl.elem.lowsc <- cbind(1,1, 2,4, 1,0,c(0:1), 2,4, 1,0,c(0:1)) xF.empl.elem.hisc <- cbind(1,1, 2,4, 4,0,c(0:1), 2,4, 4,0,c(0:1)) xF.empl.secd.lowsc <- cbind(1,1, 6,36,1,0,c(0:1), 6,36,1,0,c(0:1)) xF.empl.secd.hisc <- cbind(1,1, 6,36,4,0,c(0:1), 6,36,4,0,c(0:1)) xF.unemp.elem.lowsc <- cbind(1,1, 2,4, 1,1,c(0:1), 2,4, 1,1,c(0:1)) xF.unemp.elem.hisc <- cbind(1,1, 2,4, 4,1,c(0:1), 2,4, 4,1,c(0:1)) xF.unemp.secd.lowsc <- cbind(1,1, 6,36,1,1,c(0:1), 6,36,1,1,c(0:1)) xF.unemp.secd.hisc <- cbind(1,1, 6,36,4,1,c(0:1), 6,36,4,1,c(0:1)) fund.fd <- function(xc,b) { # fundamentalism first difference exb <- exp(xc %*% b) p <- cbind(1,exb)/(1+rowSums(exb)) return(round(100*(c(p[2,2],p[1,2],p[2,2]-p[1,2])),2)) } ## effect of unemployment on fundamentalism fund.fd(xF.mar.elem.lowsc,lcfit.cov$coeff) # 4.99 fund.fd(xF.mar.secd.lowsc,lcfit.cov$coeff) # 4.52 fund.fd(xF.mar.elem.hisc,lcfit.cov$coeff) # 4.09 fund.fd(xF.unmar.elem.lowsc,lcfit.cov$coeff) # 3.72 fund.fd(xF.mar.secd.hisc,lcfit.cov$coeff) # 3.44 fund.fd(xF.unmar.secd.lowsc,lcfit.cov$coeff) # 2.98 fund.fd(xF.unmar.elem.hisc,lcfit.cov$coeff) # 2.91 fund.fd(xF.unmar.secd.hisc,lcfit.cov$coeff) # 2.06 ## effect of marriage on fundamentalism fund.fd(xF.unemp.elem.hisc,lcfit.cov$coeff) # 5.04 fund.fd(xF.unemp.elem.lowsc,lcfit.cov$coeff) # 4.05 fund.fd(xF.empl.elem.hisc,lcfit.cov$coeff) # 3.86 fund.fd(xF.unemp.secd.hisc,lcfit.cov$coeff) # 3.37 fund.fd(xF.empl.elem.lowsc,lcfit.cov$coeff) # 2.79 fund.fd(xF.unemp.secd.lowsc,lcfit.cov$coeff) # 2.77 fund.fd(xF.empl.secd.hisc,lcfit.cov$coeff) # 1.99 fund.fd(xF.empl.secd.lowsc,lcfit.cov$coeff) # 1.22 ## ## Check whether the model fits differently for people who report no religion. ## There is a class of 22.3% that appears similar to the "secular" class, ## and another class of 19.9% that is even more secular. ## lcfit.nr <- poLCA(f,wvs.norelig,na.rm=F,nclass=4,nrep=20) # llik -24208.22 #lcfit.nr <- poLCA(f,wvs.norelig,na.rm=F,nclass=4,probs.start=lcfit.nr$probs.start,graphs=T) ## ## Also check Iran to see if same effect of education variable. ## Omit questions F102,F103,F105,F111 from the model since they were not ## asked in Iran. The effect on fundamentalism is the same, except there ## is no uptick in support at higher education levels. ## wvs.iran <- wvs[(wvs$country==364) & (wvs$age <= 40),] # select Iranians younger than 40; revolution in 1979 f.iran.cov <- cbind(C001,D019,D059,D060,D067,D076,D077,A006,D070,F063,F064,F066)~sex*educ+sex*I(educ^2)+sex*class+sex*unempl+sex*married lcfit.iran <- poLCA(f.iran.cov,wvs.iran,na.rm=F,nclass=4,nrep=5) # llik -15015.03 # be aware of rearranging of latent class ordering. plot.pv(xF.unem.marr,lcfit.cov$coeff,list(scale=c(1:8),lab="Education level")) # all respondents plot.pv(xF.unem.marr,lcfit.iran$coeff,list(scale=c(1:8),lab="Education level")) # young Iranians only ## ## Figure 2: Percent in each class, by country. ## wvs.nm <- wvs[!is.na(wvs$sex) & !is.na(wvs$educ) & !is.na(wvs$class) & !is.na(wvs$unempl) & !is.na(wvs$married),] # nm="non-missing" predclass <- lcfit.cov$predclass csum <- round(table(wvs.nm$country,predclass)/rowSums(table(wvs.nm$country,predclass)),3) csum.m <- round(table(wvs.nm$country[wvs.nm$sex==0],predclass[wvs.nm$sex==0])/rowSums(table(wvs.nm$country[wvs.nm$sex==0],predclass[wvs.nm$sex==0])),3) csum.f <- round(table(wvs.nm$country[wvs.nm$sex==1],predclass[wvs.nm$sex==1])/rowSums(table(wvs.nm$country[wvs.nm$sex==1],predclass[wvs.nm$sex==1])),3) rownames(csum) <- rownames(csum.m) <- rownames(csum.f) <- cnames colnames(csum) <- colnames(csum.m) <- colnames(csum.f) <- c("secular","fund","trad","relig") # lcfit.cov$P # 0.133 secular (1) # 0.332 fundamentalist (2) # 0.260 traditional (3) # 0.275 religious (4) par(ask=T,mar=c(7,4.1,1,1)+0.1) # Secular-Liberal class (1) o <- order(csum[,1],decreasing=T) matplot(c(1:18),cbind(csum.m[o,1],csum.f[o,1]),ylim=c(0,1),pch=c("M","F"),cex=1.2, col="black",xaxt="n",xlab="",ylab="Probability of secular-liberal class membership") points(c(1:18),csum[o,1],pch=19,cex=1.2) for (i in 1:18) { if (csum.m[o,][i,1] > csum.f[o,][i,1]) { text(i,csum.m[o,][i,1]+0.04,round(csum.m[o,][i,1],2),cex=0.8) text(i,csum.f[o,][i,1]-0.04,round(csum.f[o,][i,1],2),cex=0.8) } else { text(i,csum.m[o,][i,1]-0.04,round(csum.m[o,][i,1],2),cex=0.8) text(i,csum.f[o,][i,1]+0.04,round(csum.f[o,][i,1],2),cex=0.8) } } axis(1,1:18,labels=rownames(csum)[o],las=3); axis(2,c(0:10)/10) segments(c(1:18),csum.m[o,1],c(1:18),csum.f[o,1]) #savePlot(filename="class-secular",type="ps") #savePlot(filename="class-secular",type="png") #Fundamentalist class (2) o <- order(csum[,2],decreasing=T) matplot(c(1:18),cbind(csum.m[o,2],csum.f[o,2]),ylim=c(0,1),pch=c("M","F"),cex=1.2, col="black",xaxt="n",xlab="",ylab="Probability of fundamentalist class membership") points(c(1:18),csum[o,2],pch=19,cex=1.2) text(c(1:18),csum.m[o,2]+0.04,round(csum.m[o,2],2),cex=0.8) text(c(1:18),csum.f[o,2]-0.04,round(csum.f[o,2],2),cex=0.8) axis(1,1:18,labels=rownames(csum)[o],las=3); axis(2,c(0:10)/10) segments(c(1:18),csum.m[o,2],c(1:18),csum.f[o,2]) #savePlot(filename="class-fund.ps",type="ps") #savePlot(filename="class-fund",type="png") # Traditional class (3) o <- order(csum[,3],decreasing=T) matplot(c(1:18),cbind(csum.m[o,3],csum.f[o,3]),ylim=c(0,1),pch=c("M","F"),cex=1.2, col="black",xaxt="n",xlab="",ylab="Probability of traditional class membership") points(c(1:18),csum[o,3],pch=19,cex=1.2) for (i in 1:18) { if (csum.m[o,][i,3] > csum.f[o,][i,3]) { text(i,csum.m[o,][i,3]+0.04,round(csum.m[o,][i,3],2),cex=0.8) text(i,csum.f[o,][i,3]-0.04,round(csum.f[o,][i,3],2),cex=0.8) } else { text(i,csum.m[o,][i,3]-0.04,round(csum.m[o,][i,3],2),cex=0.8) text(i,csum.f[o,][i,3]+0.04,round(csum.f[o,][i,3],2),cex=0.8) } } axis(1,1:18,labels=rownames(csum)[o],las=3); axis(2,c(0:10)/10) segments(c(1:18),csum.m[o,3],c(1:18),csum.f[o,3]) #savePlot(filename="class-trad.ps",type="ps") #savePlot(filename="class-trad",type="png") # Religious class (4) o <- order(csum[,4],decreasing=T) matplot(c(1:18),cbind(csum.m[o,4],csum.f[o,4]),ylim=c(0,1),pch=c("M","F"),cex=1.2, col="black",xaxt="n",xlab="",ylab="Probability of religious class membership") points(c(1:18),csum[o,4],pch=19,cex=1.2) text(c(1:18),csum.m[o,4]-0.04,round(csum.m[o,4],2),cex=0.8) text(c(1:18),csum.f[o,4]+0.04,round(csum.f[o,4],2),cex=0.8) axis(1,1:18,labels=rownames(csum)[o],las=3); axis(2,c(0:10)/10) segments(c(1:18),csum.m[o,4],c(1:18),csum.f[o,4]) #savePlot(filename="class-relig.ps",type="ps") #savePlot(filename="class-relig",type="png") ############################# ## Cross-National Analysis ## ############################# # Wage gap and GDP data # ##################### # source: UN Human Development Reports, Estimated earned income (PPP US$) 2000 # and World Bank WDI 2002, GDP per capita, PPP (constant 2000 international $) # Data are from years in which WVS was fielded in each study country, except for # when data are missing, in which case they are taken from the closest available year: # Algeria 2002; Azerbaijan 1996 (missing, use 1997); Bangladesh 2002; Bosnia-Herzegovina 2001 (missing, use 2003); # India 2001; Indonesia 2001; Iran 2000; Jordan 2001; Morocco 2001; Nigeria 2000; Pakistan 2001; Singapore 2002; # Turkey 2000; Uganda 2001; Macedonia 2001 (missing, use 2002); Egypt 2000; Tanzania 2001; Montenegro 2001 (not available). wagegap <- c(0.3052,0.5963,0.5651,0.4568,0.3758,0.5104,0.2777,0.3053,0.4003,0.4242,0.3219,0.4956,0.4602,0.6587,0.5546,0.3832,0.7082,NA) logGDPpc <- c(3.7445,3.2668,3.2124,3.7725,3.3967,3.4978,3.7463,3.5892,3.5593,2.9436,3.2838,4.3725,3.8142,3.1221,3.7606,3.5483,2.7280,NA) ## ## Table 3. Cross-national regression models. ## c.womfund <- csum.f[,2] c.womfund.lo <- log(c.womfund/(1-c.womfund)) xnat <- robcov(ols(c.womfund.lo~wagegap+logGDPpc,x=T,y=T), method="efron") print(xnat) cor(wagegap,logGDPpc,use="complete.obs") ## ## Figure 3. Fitted surface from cross-national regression. ## g <- expand.grid(x=seq(0,01,0.05),y=seq(2.5,4.5,0.2)) xb <- as.matrix(cbind(1,g)) %*% matrix(coefficients(xnat)) g$z <- exp(xb)/(1+exp(xb)) pts <- data.frame(x=wagegap,y=logGDPpc,z=c.womfund) pts <- pts[-nrow(pts),] xb.hat <- as.matrix(cbind(1,pts[,-3])) %*% matrix(coefficients(xnat)) pts$hat <- exp(xb.hat)/(1+exp(xb.hat)) panel.3dscatter.yhatup <- function (x, y, z, z.hat, rot.mat = diag(4), distance, groups = NULL, type = "p", xlim.scaled, ylim.scaled, zlim.scaled, zero.scaled, col, col.point = if (is.null(groups)) plot.symbol$col else superpose.symbol$col, col.line = if (is.null(groups)) plot.line$col else superpose.line$col, lty = if (is.null(groups)) plot.line$lty else superpose.line$lty, lwd = if (is.null(groups)) plot.line$lwd else superpose.line$lwd, cex = if (is.null(groups)) plot.symbol$cex else superpose.symbol$cex, pch = if (is.null(groups)) "+" else superpose.symbol$pch, cross, ..., subscripts = TRUE) { plot.symbol <- trellis.par.get("plot.symbol") plot.line <- trellis.par.get("plot.line") superpose.symbol <- trellis.par.get("superpose.symbol") superpose.line <- trellis.par.get("superpose.line") if (!missing(col)) { col.point <- col col.line <- col } n <- length(x) col.point <- rep(col.point, length = n) col.line <- rep(col.line, length = n) lty <- rep(lty, length = n) lwd <- rep(lwd, length = n) cex <- rep(cex, length = n) pch <- rep(pch, length = n) if ("h" %in% type) { m <- ltransform3dto3d(rbind(x, y, z), rot.mat, distance) zero.scaled <- if (zero.scaled < zlim.scaled[1]) zlim.scaled[1] else if (zero.scaled > zlim.scaled[2]) zlim.scaled[2] else zero.scaled other.end <- ltransform3dto3d(rbind(x, y, z.hat), rot.mat, distance) lsegments(m[1, ], m[2, ], other.end[1, ], other.end[2, ], col = col.line, lty = lty, lwd = lwd) } } trellis.par.set("axis.line",list(col="transparent")) wireframe(z~x*y,data=g,xlab=list(label="Male-female\n wage parity"),ylab=list(label="GDP per capita \n"),xaxt="n", zlab=list(label=" Percent \n of women \n fundamentalist\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n"), xlim=c(0,1),zlim=c(0,1),screen=list(z=-55,x=-70), scales=list(y=list(at=log10(c(500,1000,2500,5000,10000,25000)),labels=c(500,1000,2500,5000,10000,25000)),arrows=F,tck=0.8,col="black"), shade=T,shade.colors=function(irr,ref,height,w=0.75) grey(w*irr+(1-w)+(1-(1-ref)^0.2)), pts=pts, panel.3d.wireframe = function(x, y, z, xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts, ...) { panel.3dwire(x = x, y = y, z = z, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) xx <- xlim.scaled[1] + diff(xlim.scaled) * (pts$x - xlim[1]) / diff(xlim) yy <- ylim.scaled[1] + diff(ylim.scaled) * (pts$y - ylim[1]) / diff(ylim) zz <- zlim.scaled[1] + diff(zlim.scaled) * (pts$z - zlim[1]) / diff(zlim) above <- (pts$z-pts$hat)>0 zz.hat <- as.vector(zlim.scaled[1] + diff(zlim.scaled) * (pts$hat - zlim[1]) / diff(zlim)) panel.3dscatter(x = xx[above], y = yy[above], z = zz[above], xlim = xlim, ylim = ylim, zlim = zlim, col="black", pch=19, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) panel.3dscatter(x = xx[!above], y = yy[!above], z = zz[!above], xlim = xlim, ylim = ylim, zlim = zlim, col="black", pch=19, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) panel.3dscatter(x = xx[!above], y = yy[!above], z = zz[!above], xlim = xlim, ylim = ylim, zlim = zlim, type="h", lty="12", lwd=2, col="black", xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) panel.3dscatter(x = xx[above], y = yy[above], z = zz.hat[above], xlim = xlim, ylim = ylim, zlim = zlim, type="h", lty="12", lwd=2, col="black", xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) panel.3dscatter.yhatup(x = xx[above], y = yy[above], z = zz[above], z.hat = zz.hat[above], xlim = xlim, ylim = ylim, zlim = zlim, type="h", lty=1, lwd=2, col="black", xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) } ) #savePlot(filename="xnat3d.ps",type="ps") #savePlot(filename="xnat3d",type="png") ## end of file.