#Chris Cochrane, June 27, 2012
library(ggplot2)
library(reshape2)
library(plyr)
library(scales)
library(grid)
library(RColorBrewer)


memory.limit(10000)

#loading Census Data 2001
AllImmigs <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/immigsonly.csv")
Muslims <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Muslims.csv")
Buddhist <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Buddhist.csv")
Catholics <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Catholics.csv")
HinduSikh <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/HinduSikh.csv")
Jews <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Jews.csv")
None <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/None.csv")
Other <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Other.csv")
Protestants <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Census/Protestants.csv")

#Variable for immigrant status 

Full$ImmigrantText[Full$Immigrant==1] <- "Immigrant"
Full$ImmigrantText[Full$Immigrant==0] <- "Native"

Full$ReligionText <- as.character(Full$Religion)

Full$Position[Full$ReligionText=="1.Muslim"] <- .53
Full$Position[Full$ReligionText=="2.Buddhist"] <- .51
Full$Position[Full$ReligionText=="3.Hindu/Sikh"] <- .40
Full$Position[Full$ReligionText=="4.Jewish"] <- -.2
Full$Position[Full$ReligionText=="5.Other"] <- -.22
Full$Position[Full$ReligionText=="6.None"] <- -.38
Full$Position[Full$ReligionText=="7.Protestant"] <- -.47
Full$Position[Full$ReligionText=="8.Catholic"] <- -.56


 #R seems incapable of rendering the graph that I want with the large data set
 
Religion <- c("1.Muslim", "2.Buddhist", "3.Hindu/Sikh", "4.Jewish", "5.Other", "6.None",
"7.Protestant", "8.Catholic")
Religion <- rep(Religion, each=2)

Religion2 <- c("Muslim", "Buddhist", "Hindu/Sikh", "Jewish", "Other", "None",
"Protestant", "Catholic")
Religion2 <- rep(Religion2, each=2)

ImmigrantText <- c("Native", "Immigrant")
ImmigrantText <- rep(Immigrant, 8)

Immigrant <- c(0,1)
Immigrant <- rep(Immigrant,8)

ImmigvNative <- c(.234, .766, .246, .754, .298, .702, .679, .320, 
.714, .286, .790, .210, .817, .183, .847, .152)

Position <- c(.43, .38, .38, -.28, -.28, -.28, -.38, -.38)
Position <- rep(Position, each=2)  
Position <- as.character(Position)
Position <- as.numeric(Position)

Full.Gimped <- data.frame(cbind(Religion, Religion2, Immigrant, ImmigrantText, ImmigvNative, Position))
Full.Gimped$Position <- as.character(Full.Gimped$Position)
Full.Gimped$Position <- as.numeric(Full.Gimped$Position)
Full.Gimped$ImmigvNative <- as.character(Full.Gimped$ImmigvNative)
Full.Gimped$ImmigvNative <- as.numeric(Full.Gimped$ImmigvNative)
Full.Gimped$Religion2 <- as.character(Full.Gimped$Religion2)

rm(Religion, Religion2, ImmigrantText, Immigrant, ImmigvNative, Position)


Immig.Graph <- ggplot(Full.Gimped, aes(Religion)) + geom_bar(subset = .(Immigrant==1), aes(y=ImmigvNative, fill=ImmigrantText)) +
    geom_bar(subset = .(Immigrant==0), aes(y = -ImmigvNative, fill=ImmigrantText)) +
    geom_text(aes(x=Religion, y=Position, label=as.character(Religion2)), angle=-90, size=3) +
  scale_y_continuous(breaks=seq(-.8, 8, .2), labels = percent_format()) +
  ylab(" ") +
  xlab("Religious Groups") +
  opts(title="Immigration Status") +
  opts(plot.title=theme_text(size =15)) +
  opts(legend.title=theme_blank()) +
  opts(axis.text.y = theme_text(size=15)) +
  opts(axis.text.x = theme_blank()) +
  geom_hline(yintercept=0) +
  opts(legend.key.size = unit(3, "mm")) +
  opts(legend.position = c(.85,.9)) +
  opts(axis.title.x = theme_text(size=15))
 
  
  


Muslim.Graph <- ggplot(Muslims, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Muslim Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.85,.70))
    
             
Buddhist.Graph <- ggplot(Buddhist, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Buddhist Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.85,.7)) +
     opts(legend.title = theme_blank())

HinduSikh.Graph <- ggplot(HinduSikh, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Hindu/Sikh Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.85,.75)) +
     opts(legend.title = theme_blank())
    
Jews.Graph <- ggplot(Jews, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Jewish Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.text = theme_text(size=10)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.15,.75)) +
    opts(legend.title = theme_blank())
    
None.Graph <- ggplot(None, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="No Religion Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.85,.75)) +
     opts(legend.title = theme_blank())
    
    
Other.Graph <- ggplot(Other, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Other Religion Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.85,.75)) +
     opts(legend.title = theme_blank())
    
    
Protestant.Graph <- ggplot(Protestants, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Protestant Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.15,.75)) +
     opts(legend.title = theme_blank())
    
    
Catholic.Graph <- ggplot(Catholics, aes(YearsinCountry)) + 
    geom_bar(aes(y = ..count../sum(..count..), fill=Birthplace)) +
    scale_x_continuous(limits=c(0, 36), breaks=seq(0,36,5) ) +
    ylab("") +
    xlab("Years in Country") +
    scale_y_continuous(limits=c(0,.12), labels=percent_format()) +
    opts(title="Catholic Immigrants") +
    opts(plot.title=theme_text(size =15)) + 
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(legend.text = theme_text(size=10)) +
    opts(legend.title = theme_text(size=10)) +
    opts(legend.key.size = unit(3, "mm")) +
    opts(legend.position = c(.15,.75)) +
    opts(legend.title = theme_blank())
    
#Multiplot function from http://wiki.stdout.org/rcookbook/Graphs/Multiple%20graphs%20on%20one%20page%20(ggplot2)/

multiplot <- function(..., plotlist=NULL, cols) {
    require(grid)

    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)

    numPlots = length(plots)

    # Make the panel
    plotCols = cols                          # Number of columns of plots
    plotRows = ceiling(numPlots/plotCols) # Number of rows needed, calculated from # of cols

    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
    vplayout <- function(x, y)
        viewport(layout.pos.row = x, layout.pos.col = y)

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
        curRow = ceiling(i/plotCols)
        curCol = (i-1) %% plotCols + 1
        print(plots[[i]], vp = vplayout(curRow, curCol ))
    }

}
    
All.Graphs <- multiplot(Immig.Graph, Muslim.Graph, Buddhist.Graph, HinduSikh.Graph, 
Jews.Graph, Other.Graph, None.Graph, Protestant.Graph, Catholic.Graph, cols=3)




CountryData <- read.csv(file = "C:/Dropbox/Projects/Muslims/Data/WVSHDI/CountryLevelData.csv")
CountryData$CountryName <- as.character(CountryData$Country)

  


#p <- ggplot(CountryData, aes(x=HDI, y=Homosexuality, size=Muslims_per)) +
#      geom_smooth(legend=FALSE) +
#      geom_point(aes(colour=Region), alpha=.9) +
#      labs(colour= "Region") +
#      geom_point(aes(size=Muslims_per), shape=1) +
#      labs(size="Muslims (%)") +
#      scale_area(range=c(2,20)) +
#      xlab("Human Development Index") +
#      ylab("Mean Opinion re: Homosexuality") +
#      opts(legend.title = theme_text(size=15)) +
#      opts(legend.text = theme_text(size=10)) +
#      opts(axis.title.y = theme_text(size=15, angle=90)) +
#      opts(axis.title.x = theme_text(size=15))


      p <- ggplot(CountryData, aes(x=HDI, y=Homosexuality)) +
      geom_smooth(legend=FALSE) +
      geom_point(aes(shape=Region, fill=Region, size=Muslims_per), alpha=.8, legend=TRUE) +
      scale_shape_manual(values=c(22, 25, 21, 24, 12, 23)) +
      labs(colour= "Region") +
      labs(size="Muslims") +
      scale_area(range=c(2,20), label=percent) +
      xlab("\nHuman Development Index") +
      ylab("Mean Opinion about Homosexuality\n") +
      opts(legend.title = theme_text(size=15)) +
      opts(legend.text = theme_text(size=15)) +
      opts(axis.title.y = theme_text(size=15, angle=90)) +
      opts(axis.title.x = theme_text(size=15)) +
      opts(axis.text.x = theme_text(size=15)) +
      opts(axis.text.y = theme_text(size=15)) +
       scale_x_continuous(breaks=c(.3, .4, .5, .6, .7, .8, .9)) 


#Coding in Stata, and then import the merged file


Ipsos <- read.csv(file="C:/Dropbox/Projects/Muslims/Data/Ipsos/Merged0611.csv")

rm(Religions, ReligionsChar, GayMarriage, Support)

Religions <- c("4.Catholic (n=24e3)", "2.Protestant (n=27e3)", "1.Muslim (n=363)", "6.Jewish (n=1065)", "3.Hindu/Sikh (n=232)", "5.Other (n=3457)", "7.No Religion (n=14e3)")
Religions <- rep(Religions, each=3)
ReligionsChar <- c("Catholic", "Protestant", "Muslim", "Jewish", "Hindu/Sikh", "Other", "No Religion")
ReligionsChar <- rep(ReligionsChar, each=3)
GayMarriage <- c(1,2,3)
GayMarriage <- rep(GayMarriage, 7)

#Support is proportion supporting + 1/2 indiff, proportion opposing + 1.2 indiff
Support <- c(.716, .169, .285, 
            .625, .173, .374, 
            .438, .121, .562, 
            .84, .105, .16, 
            .701, .106, .3,
            .785, .09, .215,
            .873, .08, .128)

Ipsos.Gimped <- data.frame(cbind(Religions, ReligionsChar, GayMarriage, Support))
Ipsos.Gimped$Support <- as.character(Ipsos.Gimped$Support)
Ipsos.Gimped$Support <- as.numeric(Ipsos.Gimped$Support)
Ipsos.Gimped$ReligionsChar <- as.character(Ipsos.Gimped$ReligionsChar)
Ipsos.Gimped$GayMarriage <- as.character(Ipsos.Gimped$GayMarriage)
Ipsos.Gimped$GayMarriage <- as.numeric(Ipsos.Gimped$GayMarriage)

attach(Ipsos.Gimped)

Ipsos.Gimped$GayMarriageChar[GayMarriage==1] <- "1.Marriage"
Ipsos.Gimped$GayMarriageChar[GayMarriage==2] <- "2.Civil Unions"
Ipsos.Gimped$GayMarriageChar[GayMarriage==3] <- "3.No Recognition"

detach(Ipsos.Gimped)


 GayMarriage.Graph <- ggplot(Ipsos.Gimped, aes(Religions)) + geom_bar(subset = .(GayMarriage==1), aes(y=Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==3), aes(y = -Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==2), aes(y = Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==2), aes(y = -Support, fill=GayMarriageChar), color="black") +  
    scale_y_continuous(breaks=seq(-.9, .9, .2), labels = percent_format()) +
  ylab(" ") +
  xlab(" ") +
  opts(legend.title=theme_blank()) +
  opts(legend.text=theme_text(size=15)) +
  opts(axis.text.y = theme_text(size=15)) +
  opts(axis.text.x = theme_text(size=15)) +
  geom_hline(yintercept=0, legend=FALSE) +
  opts(axis.title.x = theme_text(size=15)) +
  scale_fill_brewer(palette= "YlOrRd") +
    coord_cartesian(ylim= c(-.9, .9)) +
    opts(legend.position="bottom") +
      coord_flip() 
  
GayMarriage.Graph




 GayMarriage.Graph <- ggplot(Ipsos.Gimped, aes(Religions)) + geom_bar(subset = .(GayMarriage==1), aes(y=-Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==3), aes(y = Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==2), aes(y = -Support, fill=GayMarriageChar), color="black") +
    geom_bar(subset = .(GayMarriage==2), aes(y = Support, fill=GayMarriageChar), color="black") +  
    scale_y_continuous(breaks=seq(-.9, .9, .2), labels = percent_format()) +
  ylab(" ") +
  xlab(" ") +
  opts(legend.title=theme_blank()) +
  opts(legend.text=theme_text(size=15)) +
  opts(axis.text.y = theme_text(size=15)) +
  opts(axis.text.x = theme_text(size=15)) +
  geom_hline(yintercept=0, legend=FALSE) +
  opts(axis.title.x = theme_text(size=15)) +
  scale_fill_brewer(palette= "YlOrRd") +
    coord_cartesian(ylim= c(-.9, .9)) +
    opts(legend.position="bottom") +
    coord_flip() 
  
GayMarriage.Graph







#Generating data for Opposition to Gay Marriage by Church Attendance, with size
#variable for proportion of each religious group within each category of church
#attendance

#Removing all objects from workspace
rm(list=ls())

#Religion Variable, each of six religion labels repeating four times
Religion <- c("2.Catholic", "5.Protestant", "6.Muslim", "1.Jewish", "4.Hindu/Sikh", "3.Other")
Religion <- rep(Religion, each=4)

#Level of Church Attendance, each of four levels repeating for 6 religions

ChurchAtt.Level <- c("1.Less Often", "2.Monthly", "3.Weekly", "4.More Often")
ChurchAtt.Level <- rep(ChurchAtt.Level, 6)

#Proportion of each religion within each category of Church Att
Density <- c(.7645, .0859, .1246, .0250,
                    .6639, .0938, .1673, .0750,
                    .5455, .1253, .1376, .1919,
                    .8014, .0957, .0713, .0316,
                    .6450, .1641, .1374, .0534,
                    .8117, .0763, .0694, .0426)

#Proportion of each level of Churchatt within each Religion opposing all recognition of gay marriage
GayMarr.Oppose <- c(.0884, .1383, .2177, .3618,
                    .1122, .1838, .3972, .5551,
                    .3420, .50,   .5577, .5775,
                    .0377, .0577, .1316, .3030,
                    .1533, .1538, .2759, .5385,
                    .0777, .1907, .3908, .4490)

#A gimped data set for the graph
Ipsos.Gimped <- data.frame(cbind(Religion, ChurchAtt.Level, Density, GayMarr.Oppose))
Ipsos.Gimped$Religion <- as.character(Ipsos.Gimped$Religion)
Ipsos.Gimped$Density <- as.character(Ipsos.Gimped$Density)
Ipsos.Gimped$Density <- as.numeric(Ipsos.Gimped$Density)
Ipsos.Gimped$GayMarr.Oppose <- as.character(Ipsos.Gimped$GayMarr.Oppose)
Ipsos.Gimped$GayMarr.Oppose <- as.numeric(Ipsos.Gimped$GayMarr.Oppose)
label(Ipsos.Gimped$ChruchAtt.Prop) <- "Density"

p <- ggplot(Ipsos.Gimped, aes(x=ChurchAtt.Level, y=GayMarr.Oppose, size=Density), facets=~.Religion) +
      scale_y_continuous(labels=percent) +
      xlab("\nFrequency of Religious Attendance") +
      ylab("% Opposed to All Forms of Legal Recognition\n") +
      geom_point(aes(colour=Religion), show_guide=FALSE) +
      geom_point(shape=1, name="Density") +
      opts(legend.title=theme_blank()) +
      opts(legend.text=theme_text(size=15)) +
      opts(axis.text.y = theme_text(size=15)) +
      opts(axis.text.x = theme_text(size=15)) +
      opts(axis.title.x = theme_text(size=15)) +
      opts(axis.title.y = theme_text(size=15, angle=90)) +
      scale_area(range = c(2,10), labels=percent) +
      facet_grid(. ~ Religion) +
      opts(axis.text.x = theme_text(size=15, angle=90)) +
      opts(strip.text.x = theme_text(size=15)) +
      geom_smooth(aes(group=Religion, colour=Religion), show_guide=FALSE)



      

      

      #Generating data for Opposition to Gay Marriage by Education, with size
#variable for proportion of each religious group within each category of education

#Removing all objects from workspace
rm(list=ls())

#Religion Variable, each of six religion labels repeating forteen times
Religion <- c("2.Catholic", "5.Protestant", "6.Muslim", "1.Jewish", "4.Hindu/Sikh", "3.Other")
Religion <- rep(Religion, each=14)

#Level of Education, each of seven levels repeating twice for 6 religions to make room for data about native born and foreign born

Educ.Level <- c("1.<HS", "2.HS", "3.Some College", "4.College", "5.Some University", "6.University", "7.Grad/Prof Degree")
Educ.Level <- rep(Educ.Level, 12)

Immigrant <- c("Native", "Immigrant")
Immigrant <- rep(Immigrant, each = 7) # repeating immigrant and nb for each educ category
Immigrant <- rep(Immigrant, 6) # repeating immigrant for each religion

#Proportion of each religion, native born and foreign born, within each category of Education
Density <- c(.0624, .1911, .1471, .2129, .1082, .2041, .0743,
              .0291, .1063, .1499, .2015, .0992, .2466, .1674,
             .0536, .1734, .1695, .2236, .1039, .1950, .0810,
              .0410, .1281, .1662, .2128, .1118, .2073, .1328,
             .0620, .1240, .0853, .1550, .2016, .2868, .0853,
              .0250, .0857, .0821, .1464, .0786, .3714, .2017,
             .0150, .0808, .0924, .1120, .1443, .3245, .2309,
              .0244, .0447, .0854, .1301, .1179, .3049, .2927,
             .0211, .1474, .0842, .2000, .1368, .3053, .1053,
              .0060, .0714, .1250, .1488, .0833, .2857, .2798,
             .0565, .1528,.1823, .2268, .1337, .1785, .0693,
              .0274, .0822, .1429, .1879, .1311, .2407, .1879)



#Proportion of each level of Churchatt within each Religion opposing all recognition of gay marriage
GayMarr.Oppose <- c(.1974, .1453, .1300, .1066, .0833, .0704, .0746,
                     .3878, .2356, .2324, .1782, .1376, .1624, .1278,
                    .3204, .2450, .2375, .2078, .1463, .1355, .0994,
                     .4186, .2676, .2597, .2760, .1803, .2021, .2235,
                    .8750, .3750, .3333, .4737, .2917, .2188, .1250,
                     1.000, .5000, .4091, .4474, .5882, .4045, .5893,
                    .0833, .0923, .0667, .0638, .0000, .0412, .0466,
                     .0000, .1111, .1579, .1250, .1724, .0833, .0278,
                    .0000, .0769, .0000, .0769, .0000, .1111, .2000,
                     .0000, .1818, .3333, .3810, .1667, .2105, .2093,
                    .2642, .1636, .1394, .1149, .0918, .0795, .0370,
                     .3077, .2368, .2000, .2444, .0455, .1441, .0978)


Hline.Intercept <- c(.1103, .1812,
                     .1952, .2446,
                     .3448, .4858,
                     .0460, .0879,
                     .1136, .2431,
                     .1192, .1602)

Hline.Intercept <- rep(Hline.Intercept, each=7)

#A gimped data set for the graph
Ipsos.Gimped <- data.frame(cbind(Religion, Immigrant, Educ.Level, Density, GayMarr.Oppose, Hline.Intercept))
Ipsos.Gimped$Religion <- as.character(Ipsos.Gimped$Religion)
Ipsos.Gimped$Density <- as.character(Ipsos.Gimped$Density)
Ipsos.Gimped$Density <- as.numeric(Ipsos.Gimped$Density)
Ipsos.Gimped$GayMarr.Oppose <- as.character(Ipsos.Gimped$GayMarr.Oppose)
Ipsos.Gimped$GayMarr.Oppose <- as.numeric(Ipsos.Gimped$GayMarr.Oppose)
Ipsos.Gimped$Immigrant <- as.character(Ipsos.Gimped$Immigrant)
Ipsos.Gimped$Hline.Intercept <- as.character(Ipsos.Gimped$Hline.Intercept)
Ipsos.Gimped$Hline.Intercept <- as.numeric(Ipsos.Gimped$Hline.Intercept)

      p <- ggplot(Ipsos.Gimped, aes(x=Educ.Level, y=GayMarr.Oppose, size=Density, group=Immigrant)) +
      scale_y_continuous(labels=percent) +
      coord_cartesian(ylim = c(-.05, 1.05)) +
      xlab("Level of Education") +
      ylab("% Opposed to All Forms of Legal Recognition\n") +
      geom_point(aes(colour=Religion), show_guide=FALSE) +
      geom_point(shape=1, name="Density") +
      opts(legend.title=theme_blank()) +
      opts(legend.text=theme_text(size=15)) +
      opts(axis.text.y = theme_text(size=15)) +
      opts(axis.text.x = theme_text(size=15, angle=90)) +
      opts(axis.title.x = theme_text(size=15)) +
      opts(axis.title.y = theme_text(size=15, angle=90)) +
      scale_area(range = c(2,10), labels=percent) +
      facet_grid(. ~ Religion) +
      opts(axis.text.x = theme_text(size=15, angle=90)) +
      opts(strip.text.x = theme_text(size=15)) +
      opts(strip.text.y = theme_text(size=15, angle=90)) +
      stat_smooth(aes(colour=Religion), se = FALSE, show_guide=FALSE, alpha=.4) +
      facet_grid(Immigrant ~ Religion) +
      geom_hline(aes(yintercept=Hline.Intercept, colour=Religion, alpha=.9), show_guide=FALSE)











##Figure re: level and effect of Religiosity Among Muslims

Estimate.Type <- c("1.Observed Level/Observed Effect", "2.Observed Level/Average Effect", "3.Average Level/Observed Effect")

Prediction.Opposed <- c(-.4410, -.4417, -.2337)
Lower.Bound.Opposed <- c(-.3716, -.3736, -.1530)
Upper.Bound.Opposed <- c(-.5105, -.5097, -.3144)
Prediction.Support <- c(.2034, .2046, .3937)
Lower.Bound.Support <- c(.1576, .1597, .2898)
Upper.Bound.Support  <- c(.2491, .2496, .5056)



Estimates <- data.frame(cbind(Estimate.Type, Prediction.Opposed, Lower.Bound.Opposed, Upper.Bound.Opposed, Prediction.Support, Lower.Bound.Support, Upper.Bound.Support), stringsAsFactors=FALSE)
Estimates$Prediction.Opposed <- as.numeric(Estimates$Prediction.Opposed)
Estimates$Lower.Bound.Opposed <- as.numeric(Estimates$Lower.Bound.Opposed)
Estimates$Upper.Bound.Opposed <- as.numeric(Estimates$Upper.Bound.Opposed)
Estimates$Prediction.Support <- as.numeric(Estimates$Prediction.Support)
Estimates$Lower.Bound.Support <- as.numeric(Estimates$Lower.Bound.Support)
Estimates$Upper.Bound.Support <- as.numeric(Estimates$Upper.Bound.Support)

Estimates$Testing <- c("1.Observed Level\nObserved Effect", "2.Observed Level\nAverage Effect", "3.Average Level\n Observed Effect")

p <- ggplot(Estimates, aes(Estimate.Type)) +
      geom_errorbar(aes(ymin=Lower.Bound.Opposed, ymax=Upper.Bound.Opposed), width=.05, color="Red") +
      geom_point(aes(y=Prediction.Opposed), size=5, color="White") +
      geom_point(aes(y=Prediction.Opposed), size=5, shape=1, color="Red") +
      scale_y_continuous(limits = c(-.7, .7), breaks=seq(-.7,.7, .1), labels=percent) +
      geom_hline(yintercept=-.1039, color="Red", linetype="dashed") +
      ylab("Estimates\nOppose All Forms                                                              Support Marriage") +
      xlab(\n"Scenarios") + 
      geom_errorbar(aes(ymin=Lower.Bound.Support, ymax=Upper.Bound.Support), width=.05, color="Blue") +
      geom_point(aes(y=Prediction.Support), size=5, color="White") +
      geom_point(aes(y=Prediction.Support), size=5, shape=1, color="Blue") + 
      geom_hline(yintercept= .6347, color="Blue", linetype="dashed") +
      geom_hline(yintercept= 0, color="Black") +
      opts(axis.title.y = theme_text(size=15, angle=90)) +
      opts(axis.title.x = theme_text(size=15))
       
     
ggplot(Estimates, aes(Estimate.Type)) +
      geom_errorbar(aes(ymin=Lower.Bound.Opposed, ymax=Upper.Bound.Opposed), width=.1, color="Black") +
      geom_bar(aes(y=Prediction.Opposed), width=0.4, fill="Red", color="Black", alpha=.3) +
      scale_y_continuous(limits = c(-.7, .7), breaks=seq(-.7,.7, .1), labels=percent) +
      geom_hline(yintercept=-.1039, color="Red", linetype="dashed") +
      ylab("Estimates\nOppose All Forms                                                     Support Marriage") +
      xlab("\nScenarios") + 
      geom_errorbar(aes(ymin=Lower.Bound.Support, ymax=Upper.Bound.Support), width=.1, color="Black") +
        geom_bar(aes(y=Prediction.Support), width=0.4, fill="Blue", color="Black", alpha=.3) + 
      geom_hline(yintercept= .6347, color="Blue", linetype="dashed") +
      geom_hline(yintercept= 0, color="Black") +
      opts(axis.title.y = theme_text(size=15, angle=90)) +
      opts(axis.title.x = theme_text(size=15)) +
      opts(legend.text=theme_text(size=15)) +
      opts(axis.text.y = theme_text(size=15)) +
      opts(axis.text.x = theme_text(size=15)) +
      opts(axis.title.x = theme_text(size=15)) +
      opts(axis.title.y = theme_text(size=15, angle=90)) 




      ##Figure re: Effect of Length of Exposure to the Country on Opinions about Immigrants

LengthofResidence <- c(0, 1,2, 3, 4, 5, 6)
LengthofResidence.Estimate <- c(.2793, .3641, .4358, .4879, .5170, .5237, .5077)
LengthofResidence.UpperBound <- c(.3228, .3966, .4624, .5110, .5374, .5497, .5556)
LengthofResidence.LowerBound <- c(.2359, .3317, .4093, .4643, .4967, .4977, .4598)
MuslimsPer <- c(.3288, .2603, .1575, .1986, .0411, .0068, .0068) # Muslim immigs within each category
ImmigsPer <- c(.0935, .1168, .1119, .1917, .2128, .2124, .0619) # All immigs within each category



Residence <- data.frame(cbind(LengthofResidence, LengthofResidence.Estimate, LengthofResidence.UpperBound, LengthofResidence.LowerBound, MuslimsPer, ImmigsPer))
Residence$LengthofResidence.Estimate <- as.character(Residence$LengthofResidence.Estimate )
Residence$LengthofResidence.Estimate <- as.numeric(Residence$LengthofResidence.Estimate )
Residence$LengthofResidence.UpperBound <- as.character(Residence$LengthofResidence.UpperBound)
Residence$LengthofResidence.UpperBound <- as.numeric(Residence$LengthofResidence.UpperBound )
Residence$LengthofResidence.LowerBound <- as.character(Residence$LengthofResidence.LowerBound )
Residence$LengthofResidence.LowerBound <- as.numeric(Residence$LengthofResidence.LowerBound )
Residence$MuslimsPer <- as.character(Residence$MuslimsPer )
Residence$MuslimsPer <- as.numeric(Residence$MuslimsPer)
Residence$ImmigsPer <- as.character(Residence$ImmigsPer )
Residence$ImmigsPer <- as.numeric(Residence$ImmigsPer)
Residence$LengthofResidence <- as.character(LengthofResidence)
Residence$LengthofResidence <- as.numeric(LengthofResidence)


p <- ggplot(Residence, aes(LengthofResidence)) +
    geom_hline(yintercept=.61, color="Black", linetype="dotted") +
    geom_hline(yintercept=.33, color="Blue", linetype="dashed") +
    geom_hline(yintercept=.48, color="Red") +
    geom_errorbar(aes(ymin=LengthofResidence.LowerBound, ymax=LengthofResidence.UpperBound), width=.05, color="Red") +
    geom_line(aes(y=LengthofResidence.Estimate), color="Blue") +
    geom_point(aes(y=LengthofResidence.Estimate, size=MuslimsPer), color="White") +
    geom_point(aes(y=LengthofResidence.Estimate, size=MuslimsPer), shape=1, color="Blue") +
    geom_point(aes(y=LengthofResidence.Estimate, size=ImmigsPer), shape=0, color="Red") +
     scale_size_manual(name="Density") +
    scale_y_continuous(limits=c(.2,.65), labels=percent_format()) +
    xlab("Decades of Residence in Canada") +
    ylab("% Supporting Same-Sex Marriage\n") + 
    opts(axis.title.y = theme_text(size=15, angle=90)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    scale_area(range=c(2,20), name="Density", labels=percent_format()) +
    opts(legend.title = theme_text(size=15)) +
    opts(legend.text = theme_text(size=15))
   


       

##Figure RE: Effects of Education 

Education.Label <- c("0.<HS", "1.HS", "2.Some College", "3.College Diploma", "4.Some University", "5.University", "6.Grad/Prof Degree")
Education.Level <- c(0, 1, 2, 3, 4, 5, 6)
Group <- c("Native-Born Muslims", "Foreign-Born Muslims")
Group <- rep(Group, each=7)
Education.Label <- rep(Education.Label, 2)
Education.Level <- rep(Education.Level, 2)

Mus.Support.Est <- c(.1125, .1591, .2202, .2966, .3863, .4844, .5837, .1926, .2134, .2358, .2597, .2851, .312, .3403)
Mus.Support.Est.Lo <- c(0, .04, .1042, .1853, .2617, .3203, .3731, .0287, .0734, .1184, .1666, .2067, .2251, .2208)
Mus.Support.Est.Hi <- c(.226, .2782, .3363, .4078, .5108, .6485, .7943, .3565, .3554, .3531, .3528, .3636, .399, .4597)
AllCdn.Support.Est <- c(.5391, .5702, .6007, .6304, .6592, .6869, .7133)
AllCdn.Support.Est <- rep(AllCdn.Support.Est, 2)
AllCdn.Support.Est.Lo <- c(.5267, .5607, .5937, .6297, .6531, .6793, .7038)
AllCdn.Support.Est.Lo <- rep(AllCdn.Support.Est.Lo, 2)
AllCdn.Support.Est.Hi <- c(.5516, .5797, .6077, .6392, .6653, .6945, .7227)
AllCdn.Support.Est.Hi <- rep(AllCdn.Support.Est.Hi, 2)

AllCdn.Oppose.Est <- c(-.1445, -.1296, -.1161, -.1038, -.0927, -.0826, -.0736)
AllCdn.Oppose.Est <- rep(AllCdn.Oppose.Est, 2)
AllCdn.Oppose.Est.Lo <- c(-.1512, -.1346, -.1199, -.107, -.0957, -.0859, -.0772)
AllCdn.Oppose.Est.Lo <- rep(AllCdn.Oppose.Est.Lo, 2)
AllCdn.Oppose.Est.Hi <- c(-.1378, -.1246, -.1123, -.1006, -.0896, -.0793, -.0699) 
AllCdn.Oppose.Est.Hi <- rep(AllCdn.Oppose.Est.Hi, 2)
Mus.Oppose.Est <- c(-.6091, -.5108, -.4166, -.3191, -.2389, -.1738, -.1235, -.453, -.4214, -.3904, -.3603, -.3312, -.3034, -.2769)
Mus.Oppose.Est.Hi <- c(-.3385, -.2884, -.248, -.2032, -.1434, -.0794, -.0296, -.1918, -.2151, -.2355, -.2487, -.246, -.2178, -.1704)
Mus.Oppose.Est.Lo <- c(-.8797, -.7331, -.5751, -.4349, -.3345, -.2681, -.2173, -.7142, -.6277, -.5453, -.4719, -.4164, -.389, -.3835)



Education.Effects <- data.frame(cbind(Group, Education.Label, Education.Level, Mus.Support.Est, Mus.Support.Est.Lo, Mus.Support.Est.Hi, AllCdn.Support.Est, 
                      AllCdn.Support.Est.Hi, AllCdn.Support.Est.Lo, Mus.Oppose.Est, Mus.Oppose.Est.Hi, Mus.Oppose.Est.Lo, AllCdn.Oppose.Est,
                      AllCdn.Oppose.Est.Hi, AllCdn.Oppose.Est.Lo), stringsAsFactors=FALSE)


Education.Effects$Mus.Support.Est <- as.numeric(Education.Effects$Mus.Support.Est)
Education.Effects$Mus.Support.Est.Hi <- as.numeric(Education.Effects$Mus.Support.Est.Hi)
Education.Effects$Mus.Support.Est.Lo <- as.numeric(Education.Effects$Mus.Support.Est.Lo)
Education.Effects$AllCdn.Support.Est <- as.numeric(Education.Effects$AllCdn.Support.Est)
Education.Effects$AllCdn.Support.Est.Hi <- as.numeric(Education.Effects$AllCdn.Support.Est.Hi)
Education.Effects$AllCdn.Support.Est.Lo <- as.numeric(Education.Effects$AllCdn.Support.Est.Lo)
Education.Effects$Mus.Oppose.Est  <- as.numeric(Education.Effects$Mus.Oppose.Est)
Education.Effects$Mus.Oppose.Est.Hi  <- as.numeric(Education.Effects$Mus.Oppose.Est.Hi)
Education.Effects$Mus.Oppose.Est.Lo <- as.numeric(Education.Effects$Mus.Oppose.Est.Lo)
Education.Effects$AllCdn.Oppose.Est <- as.numeric(Education.Effects$AllCdn.Oppose.Est)
Education.Effects$AllCdn.Oppose.Est.Hi <- as.numeric(Education.Effects$AllCdn.Oppose.Est.Hi)
Education.Effects$AllCdn.Oppose.Est.Lo <- as.numeric(Education.Effects$AllCdn.Oppose.Est.Lo)
Education.Effects$Education.Level <- as.numeric(Education.Effects$Education.Level)


p <- ggplot(Education.Effects, aes(Education.Level)) +
    geom_ribbon(aes(ymin=Mus.Support.Est.Lo, ymax=Mus.Support.Est.Hi), color="Blue", fill="Blue", alpha=.2) +
    geom_line(aes(y=Mus.Oppose.Est)) +
    geom_line(aes(y=Mus.Oppose.Est)) +
    geom_ribbon(aes(ymin=AllCdn.Oppose.Est.Lo, ymax=AllCdn.Oppose.Est.Hi), color="Red", fill="Red", alpha=.5) +
    geom_line(aes(y=AllCdn.Oppose.Est), color="White") +
    geom_ribbon(aes(ymin=Mus.Oppose.Est.Lo, ymax=Mus.Oppose.Est.Hi), color="Blue", fill="Blue", alpha=.2) +
    geom_line(aes(y=Mus.Support.Est)) +
    geom_ribbon(aes(ymin=AllCdn.Support.Est.Lo, ymax=AllCdn.Support.Est.Hi), color="Red", fill="Red", alpha=.5) +
    geom_line(aes(y=AllCdn.Support.Est), color="White") +
    scale_y_continuous(limits=c(-.9,.9), breaks=c(-1, -.8, -.6, -.4, -.2, 0, .2, .4, .6, .8, 1), labels= percent_format()) +
    scale_x_continuous(labels = c("< HS", "HS", "< College", "College", "< Univ", "Univ", "> Univ")) +
    xlab("\nLevel of Education") +
    ylab("Estimates\nOppose All Forms                                          Support Marriage") +
    geom_hline(yintercept=0) +
    opts(axis.title.y = theme_text(size=15, angle=90)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    facet_grid(.~Group) +
    opts(strip.text.x = theme_text(size = 15)) +
    geom_text(aes(X, Y, label = Text), 
        data.frame(X=1, Y=.65, Text="All Canadians"))
      
       

    
    +
    geom_line(aes(y=FBMus.Oppose.Est), color="Green") +
    geom_line(aes(y=AllCdn.Oppose.Est), color="Red")
    
    
    geom_point(aes(y=LengthofResidence.Estimate, size=MuslimsPer), color="White") +
    geom_point(aes(y=LengthofResidence.Estimate, size=MuslimsPer), shape=1, color="Blue") +
    geom_point(aes(y=LengthofResidence.Estimate, size=ImmigsPer), shape=0, color="Red") +
     scale_size_manual(name="Density") +
    scale_y_continuous(limits=c(.2,.65), labels=percent_format()) +
    xlab("Decades of Residence in Canada") +
    ylab("% Supporting Same-Sex Marriage\n") + 
    opts(axis.title.y = theme_text(size=15, angle=90)) +
    opts(axis.title.x = theme_text(size=15)) +
    opts(axis.text.x = theme_text(size=15)) +
    opts(axis.text.y = theme_text(size=15)) +
    scale_area(range=c(2,20), name="Density", labels=percent_format()) +
    opts(legend.title = theme_text(size=15)) +
    opts(legend.text = theme_text(size=15))