# -------------------------------------------------------
# 			Micheal de Barra
# Project:	Bangladesh / disgustPaper.
# Purpose: 	Dataset and analysis for publication
# Created: 	July 23 2013
#--------------------------------------------------------

############################################### RUN FROM HERE 1ST


library("splines") 
library("survival") 
library("foreign")
library("Hmisc")
library("lsr")
library("stats") 
library("psych")
library("lattice")
library("nFactors") 
library("GPArotation") 
library("Matrix")
library("lme4") 
library("pwr")

###############################################
## LOADING DATA  ##############################
###############################################
setwd('/Users/itservices/Dropbox/Projects/BdeshDisgust/dataPublication') 
 	 
disgust <- spss.get("disgustData.por", use.value.labels=TRUE)
# orgional home: /Users/michealdebarra/Dropbox/Projects/FacePaper/Data/Matlab1FromEpiData.sav 
# names, dates and underscores removed. 

###############################################
## POWER CALC		 	   ####################
###############################################

pwr.r.test(r = .2, sig.level = .05, power = .9)	 

###############################################
## NEW VARIABLES  #############################
###############################################

## disgust vars   #############################
# 2 factor solution
disgust$fwNum <- (disgust$FW0 * 1) + (disgust$FW1 * 2) + (disgust$FW2 * 3) + (disgust$FW3 * 4)
disgust$fwNum <- factor(disgust$fwNum)


#### efa 3 below
disgust1 <- disgust[,c("HY8","C2","G3","C3","HY4","FA3","C4","HY7","FA2","G5","HY1","G2")]

disgust2 <- disgust[,c("G4", "DI7","FA1", "FA6","C1","FA4","DI4","DI3","DI6")]

disgust$DisgAv1 <- rowMeans(disgust1, na.rm=TRUE)
disgust$DisgAv2 <- rowMeans(disgust2, na.rm=TRUE)

alpha(disgust2)
alpha(disgust1)


# number of recent infections
disgust$InfectFreq <- disgust$ACUFLU1 + disgust$CHRGAST1 + disgust$ACUCGH1 + disgust$ACUVOM1 + disgust$ACUDAW1 + disgust$ACUEYE1 + disgust$ACUTTH1 + disgust$ACUSKI1 + disgust$ACUFEV1 + disgust$ACUEAR1 + disgust$ACUSIN1 + disgust$ACUDAM1 + disgust$CHRTB1


# to recode similar to stevenson. 
attach(disgust)
disgust$ACUFLU3[ACUFLU3 == 0] 					<- 4
disgust$ACUFLU3[ACUFLU3 > 0  & ACUFLU3 <= 7] 	<- 3
disgust$ACUFLU3[ACUFLU3 > 7  & ACUFLU3 <= 30] 	<- 2
disgust$ACUFLU3[ACUFLU3 > 30 & ACUFLU3 <= 365] 	<- 1
disgust$ACUFLU3[ACUFLU3 > 365]					<- 0
disgust$ACUFLU3[is.na(ACUFLU3)]					<- 0
disgust$ACUCGH3[ACUCGH3 == 0] 					<- 4
disgust$ACUCGH3[ACUCGH3 > 0  & ACUCGH3 <= 7] 	<- 3
disgust$ACUCGH3[ACUCGH3 > 7  & ACUCGH3 <= 30] 	<- 2
disgust$ACUCGH3[ACUCGH3 > 30 & ACUCGH3 <= 365] 	<- 1
disgust$ACUCGH3[ACUCGH3 > 365]					<- 0
disgust$ACUCGH3[is.na(ACUCGH3)]					<- 0
disgust$ACUVOM3[ACUVOM3 == 0] 					<- 4
disgust$ACUVOM3[ACUVOM3 > 0  & ACUVOM3 <= 7] 	<- 3
disgust$ACUVOM3[ACUVOM3 > 7  & ACUVOM3 <= 30] 	<- 2
disgust$ACUVOM3[ACUVOM3 > 30 & ACUVOM3 <= 365] 	<- 1
disgust$ACUVOM3[ACUVOM3 > 365]					<- 0
disgust$ACUVOM3[is.na(ACUVOM3)]					<- 0
disgust$ACUDAW3[ACUDAW3 == 0] 					<- 4
disgust$ACUDAW3[ACUDAW3 > 0  & ACUDAW3 <= 7] 	<- 3
disgust$ACUDAW3[ACUDAW3 > 7  & ACUDAW3 <= 30] 	<- 2
disgust$ACUDAW3[ACUDAW3 > 30 & ACUDAW3 <= 365] 	<- 1
disgust$ACUDAW3[ACUDAW3 > 365]					<- 0
disgust$ACUDAW3[is.na(ACUDAW3)]					<- 0
disgust$ACUEYE3[ACUEYE3 == 0] 					<- 4
disgust$ACUEYE3[ACUEYE3 > 0  & ACUEYE3 <= 7] 	<- 3
disgust$ACUEYE3[ACUEYE3 > 7  & ACUEYE3 <= 30] 	<- 2
disgust$ACUEYE3[ACUEYE3 > 30 & ACUEYE3 <= 365] 	<- 1
disgust$ACUEYE3[ACUEYE3 > 365]					<- 0
disgust$ACUEYE3[is.na(ACUEYE3)]					<- 0
disgust$ACUTTH3[ACUTTH3 == 0] 					<- 4
disgust$ACUTTH3[ACUTTH3 > 0  & ACUTTH3 <= 7] 	<- 3
disgust$ACUTTH3[ACUTTH3 > 7  & ACUTTH3 <= 30] 	<- 2
disgust$ACUTTH3[ACUTTH3 > 30 & ACUTTH3 <= 365] 	<- 1
disgust$ACUTTH3[ACUTTH3 > 365]					<- 0
disgust$ACUTTH3[is.na(ACUTTH3)]					<- 0
disgust$ACUSKI3[ACUSKI3 == 0] 					<- 4
disgust$ACUSKI3[ACUSKI3 > 0  & ACUSKI3 <= 7] 	<- 3
disgust$ACUSKI3[ACUSKI3 > 7  & ACUSKI3 <= 30] 	<- 2
disgust$ACUSKI3[ACUSKI3 > 30 & ACUSKI3 <= 365] 	<- 1
disgust$ACUSKI3[ACUSKI3 > 365]					<- 0
disgust$ACUSKI3[is.na(ACUSKI3)]					<- 0
disgust$ACUFEV3[ACUFEV3 == 0] 					<- 4
disgust$ACUFEV3[ACUFEV3 > 0  & ACUFEV3 <= 7] 	<- 3
disgust$ACUFEV3[ACUFEV3 > 7  & ACUFEV3 <= 30] 	<- 2
disgust$ACUFEV3[ACUFEV3 > 30 & ACUFEV3 <= 365] 	<- 1
disgust$ACUFEV3[ACUFEV3 > 365]					<- 0
disgust$ACUFEV3[is.na(ACUFEV3)]					<- 0
disgust$ACUEAR3[ACUEAR3 == 0] 					<- 4
disgust$ACUEAR3[ACUEAR3 > 0  & ACUEAR3 <= 7] 	<- 3
disgust$ACUEAR3[ACUEAR3 > 7  & ACUEAR3 <= 30] 	<- 2
disgust$ACUEAR3[ACUEAR3 > 30 & ACUEAR3 <= 365] 	<- 1
disgust$ACUEAR3[ACUEAR3 > 365]					<- 0
disgust$ACUEAR3[is.na(ACUEAR3)]					<- 0
disgust$ACUSIN3[ACUSIN3 == 0] 					<- 4
disgust$ACUSIN3[ACUSIN3 > 0  & ACUSIN3 <= 7] 	<- 3
disgust$ACUSIN3[ACUSIN3 > 7  & ACUSIN3 <= 30] 	<- 2
disgust$ACUSIN3[ACUSIN3 > 30 & ACUSIN3 <= 365] 	<- 1
disgust$ACUSIN3[ACUSIN3 > 365]					<- 0
disgust$ACUSIN3[is.na(ACUSIN3)]					<- 0
disgust$ACUDAM3[ACUDAM3 == 0] 					<- 4
disgust$ACUDAM3[ACUDAM3 > 0  & ACUDAM3 <= 7] 	<- 3
disgust$ACUDAM3[ACUDAM3 > 7  & ACUDAM3 <= 30] 	<- 2
disgust$ACUDAM3[ACUDAM3 > 30 & ACUDAM3 <= 365] 	<- 1
disgust$ACUDAM3[ACUDAM3 > 365]					<- 0
disgust$ACUDAM3[is.na(ACUDAM3)]					<- 0
detach(disgust)

disgust$ARIALLZE <- factor(disgust$ARIALLZE)


# most recent infection (- "CHRGAST3" , "CHRTB3" )
infect <- disgust[,c("ACUFLU3" , "ACUCGH3" , "ACUVOM3" , "ACUDAW3" , "ACUEYE3" , "ACUTTH3" , "ACUSKI3" , "ACUFEV3" , "ACUEAR3" , "ACUSIN3" , "ACUDAM3" )]

disgust$recentSum <- apply(infect,1,sum)
hist(disgust$recentSum)


############################################### RUN TO HERE 1ST
###############################################
## DESCRIPTIVE STATISTICS  ####################
###############################################

summary(disgust$DEMSEX)
summary(disgust$AGEATINT)
sd(disgust$AGEATINT)
table(disgust$DEMMAR)
summary(as.factor(disgust$DEMREL))


sum(disgust$ACUFLU1); 
sum(disgust$CHRGAST1); 
sum(disgust$ACUCGH1); 
sum(disgust$ACUVOM1); 
sum(disgust$ACUDAW1); 
sum(disgust$ACUEYE1); 
sum(disgust$ACUTTH1); 
sum(disgust$ACUSKI1); 
sum(disgust$ACUFEV1); 
sum(disgust$ACUEAR1); 
sum(disgust$ACUSIN1); 
sum(disgust$ACUDAM1); 
sum(disgust$CHRTB1)


mean(disgust$DIAALLZE)
SD(disgust$DIAALLZE)

mean(as.numeric(disgust$ARIALLZE))
SD(disgust$ARIALLZE)

attach(disgust)
cor.test(InfectFreq,as.numeric(ARIALLZE), method = "s")	 
cor.test(InfectFreq,DIAALLZE, method = "s")	 
cor.test(as.numeric(ARIALLZE),DIAALLZE, method = "s")
detach(disgust)
###############################################
## FACTOR ANALYSIS 1 - first run   ############
###############################################

# http://www.statmethods.net/advstats/factor.html


attach(disgust)
items <- c(
"C1", #Touching an animal
"C2", #Touching the inside of a toilet
"C3", #Eating from dirty plate
"C4", #Accidentally using other persons toothbrush
"C5", #Eating dropped sweet
"C6", #A stranger touching your things
"C7", #Old/dirty money
"DI1", # blind person
"DI2", #Single legged woman
"DI3", #Very obese man
"DI4", # unkempt beggar
"DI5", #Person coughing with mucus
"DI6", #Child with diarrhoea 
"DI7", #Deformed body
"FA1", #Perished / decomposed fish
"FA2", #Dead animal 
"FA3", #Sour milk
"FA4", #Snake meat 
"FA5", #Eating last nights food
"FA6", #Animal feces in yard
"FA7", #Over-ripe fruit
"G2", #Infected eye
"G3", #Skin with scabies
"G4", #Hand without a finger
"G5", #Small acne
"HY1", #Person who never washes himself
"HY2", #Human feces in yard
"HY3", #Urination beside household 
"HY4", #Spit on the road
"HY5", #Leaving dirty toilet 
"HY6", #Untidy house-hold/family
"HY7", #Eating something with left hand
"HY8" #Picking your nose
# "O1" #Hot boiled rice
#O3	 	# Pure water
#O2		# Dry road
)
detach(disgust)
disgItems <- disgust[,items]
disgItemsNA <- na.omit(disgItems)
# reorder columns
disgItemsNA[,order(names(disgItemsNA))]

names(disgItemsNA) <- c("C1 Touching an animal" , "C2 Touching the inside of a toilet" , "C3 Eating from dirty plate" , "C4 Accidentally using other persons toothbrush" , "C5 Eating dropped sweet" , "C6 A stranger touching your things" , "C7 Old/dirty money" , "DI1  blind person" , "DI2 Single legged woman" , "DI3 Very obese man" , "DI4  unkempt beggar" , "DI5 Person coughing with mucus" , "DI6 Child with diarrhoea " , "DI7 Deformed body" , "FA1 Perished / decomposed fish" , "FA2 Dead animal " , "FA3 Sour milk" , "FA4 Snake meat " , "FA5 Eating last nights food" , "FA6 Animal feces in yard" , "FA7 Over-ripe fruit" , "G2 Infected eye" , "G3 Skin with scabies" , "G4 Hand without a finger" , "G5 Small acne" , "HY1 Person who never washes himself" , "HY2 Human feces in yard" , "HY3 Urination beside household " , "HY4 Spit on the road" , "HY5 Leaving dirty toilet " , "HY6 Untidy house-hold/family" , "HY7 Eating something with left hand" , "HY8 Picking your nose")



## PARALLEL ANALYSIS				   ############
ev <- eigen(cor(disgItemsNA))
ap <- parallel(subject=nrow(disgItemsNA),var=ncol(disgItemsNA),rep=100,cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)
ap
# two factors. 
ev

## EFA								   ############

fa(disgItemsNA, nfactors=2,  rotate="oblimin")

# communality < 2: C7 DI1 DI2 DI5 FA4 FA7 HY2 HY3 HY5 HY6



########################################################################
## FACTOR ANALYSIS 2 - omitting low communality items  #################
########################################################################

# http://www.statmethods.net/advstats/factor.html


attach(disgust)
items <- c(
"C1", #Touching an animal
"C2", #Touching the inside of a toilet
"C3", #Eating from dirty plate
"C4", #Accidentally using other persons toothbrush
"C5", #Eating dropped sweet
"C6", #A stranger touching your things
#"C7", #Old/dirty money
#"DI1", # blind person
#"DI2", #Single legged woman
"DI3", #Very obese man
"DI4", # unkempt beggar
#"DI5", #Person coughing with mucus
"DI6", #Child with diarrhoea 
"DI7", #Deformed body
"FA1", #Perished / decomposed fish
"FA2", #Dead animal 
"FA3", #Sour milk
#"FA4", #Snake meat 
"FA5", #Eating last nights food
"FA6", #Animal feces in yard
#"FA7", #Over-ripe fruit
"G2", #Infected eye
"G3", #Skin with scabies
"G4", #Hand without a finger
"G5", #Small acne
"HY1", #Person who never washes himself
#"HY2", #Human feces in yard
#"HY3", #Urination beside household 
"HY4", #Spit on the road
#"HY5", #Leaving dirty toilet 
#"HY6", #Untidy house-hold/family
"HY7", #Eating something with left hand
"HY8" #Picking your nose
# "O1" #Hot boiled rice
#O3	 	# Pure water
#O2		# Dry road
)
detach(disgust)
disgItems <- disgust[,items]
disgItemsNA <- na.omit(disgItems)
# reorder columns
disgItemsNA[,order(names(disgItemsNA))]

names(disgItemsNA) <- c("C1 Touching an animal" , "C2 Touching the inside of a toilet" , "C3 Eating from dirty plate" , "C4 Accidentally using other persons toothbrush" , "C5 Eating dropped sweet" , "C6 A stranger touching your things"  , "DI3 Very obese man" , "DI4  unkempt beggar" , "DI6 Child with diarrhoea " , "DI7 Deformed body" , "FA1 Perished / decomposed fish" , "FA2 Dead animal " , "FA3 Sour milk" , "FA5 Eating last nights food" , "FA6 Animal feces in yard" , "G2 Infected eye" , "G3 Skin with scabies" , "G4 Hand without a finger" , "G5 Small acne" , "HY1 Person who never washes himself" , "HY4 Spit on the road" , "HY7 Eating something with left hand" , "HY8 Picking your nose")



## PARALLEL ANALYSIS				   ############
ev <- eigen(cor(disgItemsNA))
ap <- parallel(subject=nrow(disgItemsNA),var=ncol(disgItemsNA),rep=100,cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)
ap
# two factors. 
ev

## EFA								   ############

fa(disgItemsNA, nfactors=2,  rotate="oblimin")





########################################################################
## FACTOR ANALYSIS 3 - omitting low loading items  #################
#######################################################################

# http://www.statmethods.net/advstats/factor.html


attach(disgust)
items <- c(
"C1", #Touching an animal
"C2", #Touching the inside of a toilet
"C3", #Eating from dirty plate
"C4", #Accidentally using other persons toothbrush
#"C5", #Eating dropped sweet
#"C6", #A stranger touching your things
#"C7", #Old/dirty money
#"DI1", # blind person
#"DI2", #Single legged woman
"DI3", #Very obese man
"DI4", # unkempt beggar
#"DI5", #Person coughing with mucus
"DI6", #Child with diarrhoea 
"DI7", #Deformed body
"FA1", #Perished / decomposed fish
"FA2", #Dead animal 
"FA3", #Sour milk
#"FA4", #Snake meat 
"FA5", #Eating last nights food
"FA6", #Animal feces in yard
#"FA7", #Over-ripe fruit
"G2", #Infected eye
"G3", #Skin with scabies
"G4", #Hand without a finger
"G5", #Small acne
"HY1", #Person who never washes himself
#"HY2", #Human feces in yard
#"HY3", #Urination beside household 
"HY4", #Spit on the road
#"HY5", #Leaving dirty toilet 
#"HY6", #Untidy house-hold/family
"HY7", #Eating something with left hand
"HY8" #Picking your nose
# "O1" #Hot boiled rice
#O3	 	# Pure water
#O2		# Dry road
)
detach(disgust)
disgItems <- disgust[,items]
disgItemsNA <- na.omit(disgItems)
# reorder columns
disgItemsNA[,order(names(disgItemsNA))]

names(disgItemsNA) <- c("C1 Touching an animal" , "C2 Touching the inside of a toilet" , "C3 Eating from dirty plate" , "C4 Accidentally using other persons toothbrush" ,  "DI3 Very obese man" , "DI4  unkempt beggar" , "DI6 Child with diarrhoea " , "DI7 Deformed body" , "FA1 Perished / decomposed fish" , "FA2 Dead animal " , "FA3 Sour milk" , "FA5 Eating last nights food" , "FA6 Animal feces in yard" , "G2 Infected eye" , "G3 Skin with scabies" , "G4 Hand without a finger" , "G5 Small acne" , "HY1 Person who never washes himself" , "HY4 Spit on the road" , "HY7 Eating something with left hand" , "HY8 Picking your nose")



## PARALLEL ANALYSIS				   ############
ev <- eigen(cor(disgItemsNA))
ap <- parallel(subject=nrow(disgItemsNA),var=ncol(disgItemsNA),rep=100,cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)
ap
# two factors. 
ev

## EFA								   ############

fa(disgItemsNA, nfactors=2,  rotate="oblimin")


## CALCULATE ALPHA				   ############
head(disgust)



###############################################
## CORRELATIONS Disgust            ############
###############################################

# two disgust factors
r <- matrix(c(disgust$DisgAv1,disgust$DisgAv2),ncol=2)
rcorr(r,type="pearson")

# disgust and sex
t.test(disgust$DisgAv1~disgust$DEMSEX, var.equal = TRUE)
cohensD(disgust$DisgAv1~disgust$DEMSEX)

t.test(disgust$DisgAv2~disgust$DEMSEX, var.equal = TRUE)
cohensD(disgust$DisgAv2~disgust$DEMSEX)

# FWG ANOVA
fit <- aov(DisgAv1 ~ fwNum, data=disgust)
# layout(matrix(c(1,2,3,4),2,2))
# plot(fit)
summary(fit)

fit <- aov(DisgAv2 ~ fwNum, data=disgust)
layout(matrix(c(1,2,3,4),2,2))
plot(fit)
summary(fit)

###############################################
## DISGUST AND CURRENT DISEASE     ############
###############################################


# http://www.r-bloggers.com/linear-mixed-models-in-r/
# http://www.unt.edu/rss/class/Jon/Benchmarks/LinearMixedModels_JDS_Dec2010.pdf

# Flu Gastric Pain Cough Vommiting Diarrhea Eye Infection Tooth Ache Skin Infection Fever Ear Ache Sinus Infection Dysentery Genitourinary Infection Tuberculosis


m1.lmm = lmer(InfectFreq ~ (1|fwNum), data=disgust)
summary(m1.lmm)
lower <- coef(summary(m1.lmm))[,1] + qnorm(.025)*coef(summary(m1.lmm))[,2]
upper <- coef(summary(m1.lmm))[,1] + qnorm(.975)*coef(summary(m1.lmm))[,2]
cbind(coef(summary(m1.lmm)), lower, upper)


m2.lmm = lmer(InfectFreq ~ DisgAv1 + DisgAv2 + (1|fwNum), data=disgust)
summary(m2.lmm)
lower <- coef(summary(m2.lmm))[,1] + qnorm(.025)*coef(summary(m2.lmm))[,2]
upper <- coef(summary(m2.lmm))[,1] + qnorm(.975)*coef(summary(m2.lmm))[,2]
cbind(coef(summary(m2.lmm)), lower, upper)


m3.lmm = lmer(InfectFreq ~ DisgAv1 + DisgAv2 + DEMSEX + AGEATINT + (1|fwNum), data=disgust)
summary(m3.lmm)
lower <- coef(summary(m3.lmm))[,1] + qnorm(.025)*coef(summary(m3.lmm))[,2]
upper <- coef(summary(m3.lmm))[,1] + qnorm(.975)*coef(summary(m3.lmm))[,2]
cbind(coef(summary(m3.lmm)), lower, upper)



## DISGUST AND DISEASE RECENCY   ############

m1.lmm = lmer(recentSum ~ (1|fwNum), data=disgust)
summary(m1.lmm)
lower <- coef(summary(m1.lmm))[,1] + qnorm(.025)*coef(summary(m1.lmm))[,2]
upper <- coef(summary(m1.lmm))[,1] + qnorm(.975)*coef(summary(m1.lmm))[,2]
cbind(coef(summary(m1.lmm)), lower, upper)


m2.lmm = lmer(recentSum ~ DisgAv1 + DisgAv2 + (1|fwNum), data=disgust)
summary(m2.lmm)
lower <- coef(summary(m2.lmm))[,1] + qnorm(.025)*coef(summary(m2.lmm))[,2]
upper <- coef(summary(m2.lmm))[,1] + qnorm(.975)*coef(summary(m2.lmm))[,2]
cbind(coef(summary(m2.lmm)), lower, upper)

m3.lmm = lmer(recentSum ~ DisgAv1 + DisgAv2 + DEMSEX + AGEATINT +  (1|fwNum), data=disgust)
summary(m3.lmm)
lower <- coef(summary(m3.lmm))[,1] + qnorm(.025)*coef(summary(m3.lmm))[,2]
upper <- coef(summary(m3.lmm))[,1] + qnorm(.975)*coef(summary(m3.lmm))[,2]
cbind(coef(summary(m3.lmm)), lower, upper)


###############################################
## DISGUST AND CHILDHOOD DISEASE   ############
###############################################


fit <- lm((DisgAv1) ~ (DIAALLZE) + as.numeric(ARIALLZE)  + fwNum + AGEATINT + DEMSEX, data=disgust) 			
summary(fit)
confint(fit, level=0.95)

fit <- lm((DisgAv2) ~ (DIAALLZE) + as.numeric(ARIALLZE) + fwNum + AGEATINT + DEMSEX, data=disgust) 			
summary(fit)
confint(fit, level=0.95)