---
title: "Introduction to R: Final Project"
output: word_document
---

## Demographic and Health Survey
The Demographic and Health Surveys are nationally representative population-based serial cross-sectional surveys sponsored by the United States Agency for International Development (USAID) and administered by ICF International.  Currently data is collected from 44 countries/regions in sub-Saharan Africa as well as countries in North Africa, West Asia, Europe, South and Southeast Asia, Latin America and the Caribbean.  Households in these regions are selected using a probability sample from census frames or, in cases where no census frame exists, from a complete list of villages or communities.  The researchers employed a two-stage cluster sampling procedure where a cluster of households was randomly selected from the sampling frame.  Subsequent to the household interviews, every eligible woman and man were interviewed.  The domestic violence module of the DHS was only administered to females in a subset of the countries that participated in the DHS.  Eligibility for the domestic violence module included being able to administer the survey in the absence of other people and currently/previously been married and/or living with a man that was her intimate partner.

## Sample
  The present study includes countries in sub-Saharan Africa that participated in the sixth DHS survey (DHS-6) between 2010-2014 and were selected to receive the domestic violence module, including information on alcohol use.  Hence, the study contains DHS-6 data from Burkina Faso, Cameroon, Comoros, the Democratic Republic of the Congo, Gabon, Cote d’Ivoire, Malawi, Mali, Mozambique, Nigeria, Sierra Leone, Tanzania, Uganda and Zimbabwe.  This resulted in a final sample of 103,453 women aged 15-49 years.
  
## Measures
	Information on basic demographics was taken from the general women’s questionnaire and includes data on current marital status, residency type, wealth and a literacy assessment.  Data on partner’s alcohol use and intimate partner violence (IPV) perpetrated by the husband/partner against the participant was collected in the domestic violence module.  For the items included in this study, all questions were asked of the female participant in regards to the behaviors of her current or most recent husband/partner.  Items assessing the participant’s current or most recent partner’s controlling behaviors, emotional violence, physical (severe and less severe) violence and sexual violence were collected.   

## Research Aims
  For the purpose of this project, I plan to:
    1. Characterize the sample based on social and demographic variables
    2. Examine the prevalence of alcohol use and various forms of IPV
    3. Evaluate the relationships between socio-demographic variables, alcohol use and IPV


```{r, cache=FALSE}
dhs = read.csv(file = "/Users/Desktop/DHS_CSV/DHSMERGE.csv", header=TRUE, as.is=TRUE, na.strings="")

```

To begin, I removed variables that are irrelevant to my research question.

```{r, cache=FALSE}
dhsvars = names(dhs) %in% c("D113", "country", "D104", "D106", "D107", "D108", "V149", "V501", "V140", "V155", "V190", "V447A", "D103A", "D103B", "D103C", "D105A", "D105B", "D105C", "D105D", "D105E", "D105G", "D105H", "D105I", "D102", "D101A", "D101B", "D101C", "D101D", "D101E", "D101F", "D114")
dhs = dhs[dhsvars]
```

I then added in variable names for the remaining `r ncol(dhs)` variables in my dataset.

```{r, cache=FALSE}
ColNames = c("Residence", "Education", "Literacy", "Wealth", "Age", "Marital_Status",  "Control_Jealous", "Control_Unfaithful", "Control_Friends", "Control_Family", "Control_Whereabouts", "Control_Trust","Control_N", "EV_Humiliate", "EV_Threaten", "EV_Insult", "Emotional_Violence", "PV_Push", "PV_Slap", "PV_Twist", "PV_Kick", "PV_StrangleBurn", "PV_Weapon", "SV_Sex", "SV_OthSex",  "Less_Severe_Violence", "Severe_Violence", "Sexual_Violence", "Alcohol", "Drunk_Freq", "Country" )
length(ColNames)
colnames(dhs) = ColNames
```

The majority of variables in this dataset are categorical and required some recoding.  In the following code chunk I also recoding missing data as NA:

```{r, cache=FALSE}
dhs[(dhs==7) | (dhs==8) | (dhs==9)] = NA
dhs$Urban = ifelse(dhs$Residence==2, 0, 1)
dhs$SomeEduc = ifelse(dhs$Education>=1, 1, 0)
dhs$Married = ifelse(dhs$Marital_Status==1, 2, ifelse(dhs$Marital_Status==2, 1, 0))
dhs$Literacy[(dhs$Literacy==3) | (dhs$Literacy==4)] = NA
dhs$Literacy = ifelse(dhs$Literacy==0, 0, ifelse(dhs$Literacy==1, 1, 2))
dhs$Country = ifelse(dhs$Country==0, "Burkina_Faso", ifelse(dhs$Country==1, "Cameroon", ifelse(dhs$Country==2, "Comoros", ifelse(dhs$Country==3, "DRC", ifelse(dhs$Country==4, "Gabon", ifelse(dhs$Country==5, "Ivory_Coast", ifelse(dhs$Country==6, "Malawi", ifelse(dhs$Country==7, "Mali", ifelse(dhs$Country==8, "Mozambique", ifelse(dhs$Country==9, "Nigeria", ifelse(dhs$Country==10, "Sierra_Leone", ifelse(dhs$Country==11, "Tanzania", ifelse(dhs$Country==12, "Uganda", "Zimbabwe")))))))))))))
```

##Aim 1: Characterize the sample based on social and demographic variables
To characterize the sample, I described the average age (in years), % illiterate, % residing in a rural locale, % without any formal education, and % divorced/widowed/separated stratified by country.

```{r, cache=FALSE}
Index = split(1:nrow(dhs), dhs$Country)
lapply(Index, head)
#Age
T1_Age_Mean = tapply(dhs$Age, dhs$Country, mean)
T1_Age_SD = tapply(dhs$Age, dhs$Country, sd)
table1 = matrix(paste(signif(T1_Age_Mean, 4), " (SD=", signif(T1_Age_SD,2),")",sep=""), nrow =1, ncol=11)
colnames(table1) = list("Burkina Faso", "Cameroon", "Comoros", "DRC", "Gabon", "Ivory Coast", "Malawi", "Sierra Leone", "Tanzania", "Uganda", "Zimbabwe")
rownames(table1) = list("Age (Yrs)")

#Education
T1_Educ = sapply(Index, function(x) table (dhs$SomeEduc[x]))
T1_Educ_P = (signif(prop.table(T1_Educ,2),3))*100
table1 = rbind(table1, T1_Educ_P[1,])
rownames(table1)[nrow(table1)] = "% Uneducated"

#Residence
T1_Urban = sapply(Index, function(x) table (dhs$Urban[x]))
T1_Urban_P = (signif(prop.table(T1_Urban,2),3))*100
table1 = rbind(table1, T1_Urban_P[1,])
rownames(table1)[nrow(table1)] = "% Rural"

#Marital Status
T1_Married = sapply(Index, function(x) table (dhs$Married[x]))
T1_Married_P = (signif(prop.table(T1_Married,2),3))*100
table1 = rbind(table1, T1_Married_P[1,])
rownames(table1)[nrow(table1)] = "% Divorced/Widowed/Separated"

#Literacy
T1_Literacy = sapply(Index, function(x) table (dhs$Literacy[x]))
T1_Literacy_P = (signif(prop.table(T1_Literacy,2),3))*100
table1 = rbind(table1, T1_Literacy_P[1,])
rownames(table1)[nrow(table1)] = "% Illiterate"

#Adding p-values
agereg = lm(formula = Age~Country, data=dhs) #p<.001
chisq.test(dhs$SomeEduc, dhs$Country, correct=TRUE, p=rep(1/length(dhs$Literacy), length(dhs$Literacy)), rescale.p=FALLSE, simulate.p.value=FALSE, B=2000)
chisq.test(dhs$Urban, dhs$Country, correct=TRUE, p=rep(1/length(dhs$Literacy), length(dhs$Literacy)), rescale.p=FALLSE, simulate.p.value=FALSE, B=2000)
chisq.test(dhs$Married, dhs$Country, correct=TRUE, p=rep(1/length(dhs$Literacy), length(dhs$Literacy)), rescale.p=FALLSE, simulate.p.value=FALSE, B=2000)
chisq.test(dhs$Literacy, dhs$Country, correct=TRUE, p=rep(1/length(dhs$Literacy), length(dhs$Literacy)), rescale.p=FALLSE, simulate.p.value=FALSE, B=2000)
pvalues = c("<.001", "<.001", "<.001", "<.001", "<.001")
dim(table1)
table1 = cbind(table1, pvalues)
colnames(table1)[ncol(table1)] = "p-value"
```

##Table 1: Social and demographic characteristics by country
```{r, cache=FALSE}
table1
```

To visualize this data, I created a box plot for age by country and bar plots with proportions for the remaining socio-demographic variables by country.  Several of the variables had more than one level.  For the purpose of the bar plots, I allowed for more than 2 levels to be displayed.  Also, several of the countries did not contain some of these questions in their survey which is why you will notice that some of the variables are missing data on entire countries.  For the purpose of this project, I'll assume that there are no systematic differences between countries or within-country correlation that could invalidate the inferential analyses.


```{r, cache=FALSE}
#Graphs
#Age
boxplot(dhs$Age ~ dhs$Country, las=3, main="Distribution of Age by Country", ylab="Age(Yrs)", cex.lab=0.6) #Not all countries asked about age

#Education
tableeduc = table(dhs$SomeEduc, dhs$Country)
barplot(prop.table(tableeduc, 2), las=3, main="% Educated by Country", ylab="Proportion", col=c("grey", "white"), legend=(c("No Education", "Some Education")), cex.names=0.6)

#Residence
tableurban = table(dhs$Urban, dhs$Country)
barplot(prop.table(tableurban, 2), las=3, main="Place of Residence by Country", ylab="Proportion", col=c("grey", "white"), legend=(c("Rural", "Urban")), cex.names=0.6)

#Marital Status
tablemarried = table(dhs$Married, dhs$Country)
barplot(prop.table(tablemarried, 2), las=3, main="Marital Status by Country", ylab="Proportion", col=c("darkgrey", "grey", "white"), legend=(c("Widowed/Divorced/Separated", "Co-Habitating", "Married")), cex.names=0.6)

#Literacy
tableliteracy = table(dhs$Literacy, dhs$Country)
barplot(prop.table(tableliteracy, 2), las=3, main="Literacy by Country", ylab="Proportion", col=c("darkgrey", "grey", "white"), legend=(c("Illiterate", "Somewhat Literate", "Literate")), cex.names=0.6)
```

Results from table 1 and the graphs above reveal that the average age of participants was between 28.2-29.5 years.  The majority of participants had received education and were literate, with the exception of respondents in Burkina Faso, Ivory Coast, and Sierra Leone.  The majority of participants lived in rural areas except in Gabon.  The majority of participants were married or co-habitating with their husband/partner.  There is clearly some variation in several of these sociodemographic indicators by country.  

##Aim 2: Examine the prevalence of alcohol use and various forms of IPV
To explore the prevalence of intimate partner violence (IPV) and partner's alcohol use, I calculated the prevalence of each type of violence and partner's alcohol use.
```{r, cache=FALSE}
#Calculating Prevalence
EV = table(dhs$Emotional_Violence)
EV_Prev = EV[2]/(EV[1]+EV[2])
LSV = table(dhs$Less_Severe_Violence)
LSV_Prev = LSV[2]/(LSV[1]+LSV[2])
PV = table(dhs$Severe_Violence)
PV_Prev = PV[2]/(PV[1]+PV[2])
SV = table(dhs$Sexual_Violence)
SV_Prev = SV[2]/(SV[1]+SV[2])
Alc = table(dhs$Alcohol)
Alc_Prev = Alc[2]/(Alc[1]+Alc[2])

#Creating Table 2 by IPV Type
table2 = matrix(paste(EV), nrow=1, ncol=2)
rownames(table2) = list("Emotional Violence")
table2 = rbind(table2, LSV)
rownames(table2)[nrow(table2)] = "Physical Violence (Less Severe)"
table2 = rbind(table2, PV)
rownames(table2)[nrow(table2)] = "Physical Violence (Severe)"
table2 = rbind(table2, SV)
rownames(table2)[nrow(table2)] = "Sexual Violence"
table2 = rbind(table2, Alc)
rownames(table2)[nrow(table2)] = "Partner's Alcohol Use"
Prevalence = c(EV_Prev, LSV_Prev, PV_Prev, SV_Prev, Alc_Prev)
table2 = cbind(table2, (Prevalence*100))
colnames(table2) = list("No", "Yes", "Prevalence (%)")
```

##Table 2: Prevalence of IPV and Partner's Alcohol Use by Country
```{r, cache=FALSE}
table2
```

As can be seen from Table 2 above, the prevalence of partner's alcohol use as reported by the female respondent was 32.4%.  The prevalence of the types of intimate partner violence are as follows:

Emotional Violence: 24.96%
Physical Violence (Less Severe): 25.44%
Physical Violence (Severe): 8.85%
Sexual Violence: 9.87%

```{r, cache=FALSE}
#Plot of IPV Prevalence by Type
prevmatrix = matrix(c(24.96, 25.44, 8.85, 9.87), 1, 4, byrow=TRUE)
colnames(prevmatrix) = c("EV", "PV", "SPV", "SV")
barplot(prevmatrix, main="Prevalence of IPV by Type", ylab = "Prevalence", ylim=c(0, 100))
```
Footnote: EV=Emotional Violence, PV = Physical Violence, SPV = Severe Physical Violence, SV=Sexual Violence

##Aim 3:Evaluate the relationships between socio-demographic variables, alcohol use and IPV

To investigate Aim 3, I examined the univariate relationship between partner's alcohol use and IPV using chi-squared analyses.  Subsequently, I examined this relationship controlling for potential confounders of this relationship including current marital status, literacy, education level, place of residence (urban vs. rural) and age.  While not all countries administered questions for all of the aforementioned potential confounders, for the purpose of this project we are going to assume that missing data is not systematic or informative.

```{r, cache=FALSE}
#Chi-Squared Analyses

#Emotional Violence
EV_ChiSq = chisq.test(dhs$Emotional_Violence, dhs$Alcohol)
names(EV_ChiSq)
table3 = matrix(paste(round(EV_ChiSq$statistic, digits=3)), nrow=1, ncol=1)

#Less Severe Physical Violence
LSV_ChiSq = chisq.test(dhs$Less_Severe_Violence, dhs$Alcohol)
table3 = rbind(table3, round(LSV_ChiSq$statistic, digits=3))

#Severe Physical Violence
PV_ChiSq = chisq.test(dhs$Severe_Violence, dhs$Alcohol)
table3 = rbind(table3, round(PV_ChiSq$statistic, digits=3))

#Sexual Violence
SV_ChiSq = chisq.test(dhs$Sexual_Violence, dhs$Alcohol)
table3 = rbind(table3, round(SV_ChiSq$statistic, digits=3))

pval_ChiSq = signif(c(EV_ChiSq$p.value, LSV_ChiSq$p.value, PV_ChiSq$p.value, SV_ChiSq$p.value), digits=3)
table3 = cbind(table3, c("<.001", "<.001", "<.001", "<.001"))
colnames(table3) = list("X-Squared", "p-value")
rownames(table3) = list("Emotional Violence", "Physical Violence(Less Severe)", "Physical Violence(Severe)", "Sexual Violence")
```

##Table 3: Univariate Relationship Between Partner's Alcohol Use and IPV
```{r, cache=FALSE}
table3
```

##Logistic Regression Models Evaluating the Relationship Between Partner's Alcohol Use and IPV Controlling for Education, Literacy, Marital Status, Residence, Age and Wealth

The first model evaluated the relationship between partner's alcohol use and emotional violence

```{r, cache=FALSE}
#Emotional Violence
EVlogistic = glm(Emotional_Violence ~ factor(Alcohol) + factor(SomeEduc) + factor(Literacy) + factor(Married) + factor(Urban) + Wealth + Age, data=dhs, family=binomial())
summary(EVlogistic)
OR_EV = exp(coef(EVlogistic))
OR_EV
```
Results from this model suggest that partner's alcohol use was associated with increased odds of Emotional Violence (OR=2.58, p<0.001).

The second model evaluated the relationship between partner's alcohol use and less severe physical violence.
```{r, cache=FALSE}
#Less Severe Physical Violence
LSVlogistic = glm(Less_Severe_Violence ~ factor(Alcohol) + factor(SomeEduc) + factor(Literacy) + factor(Married) + factor(Urban) + Wealth + Age, data=dhs, family=binomial())
summary(LSVlogistic)
OR_LSV = exp(coef(LSVlogistic))
OR_LSV
```
Results from this model suggest that partner's alcohol use was associated with increased odds of physical violence (less severe type; OR=23.18, p<0.001).

The third model evaluated the relationship between partner's alcohol use and severe physical violence.
```{r, cache=FALSE}
#Severe Physical Violence
PVlogistic = glm(Severe_Violence ~ factor(Alcohol) + factor(SomeEduc) + factor(Literacy) + factor(Married) + factor(Urban) + Wealth + Age, data=dhs, family=binomial())
summary(PVlogistic)
OR_PV = exp(coef(PVlogistic))
OR_PV
```
Results from this model suggest that partner's alcohol use was associated with increased odds of severe physical violence (OR=3.34, p<0.001).

The fourth model evaluated the relationship between partner's alcohol use and sexual violence.
```{r, cache=FALSE}
#Sexual Violence
SVlogistic = glm(Sexual_Violence ~ factor(Alcohol) + factor(SomeEduc) + factor(Literacy) + factor(Married) + factor(Urban) + Wealth + Age, data=dhs, family=binomial())
summary(SVlogistic)
OR_SV = exp(coef(SVlogistic))
OR_SV
```
Results from this model suggest that partner's alcohol use was associated with increased odds of sexual violence (OR=2.34, p<0.001).

In summary, alcohol use by ones partner appeared to increase the odds of emotional, physical and sexual violence.  