# Title  Sex and Religion  7 June 2014
# This R script for "Religion and Parental Cooperation: an empirical test of Slone’s sexual signaling model"
# The study will appear in a collection of essays edited by Jason Slone, (information to follow)
# For information please contact: Joseph Bulbulia:

# Bulbulia, J. A., Shaver, J., *Greaves, L., Sosis, R., & Sibley, C. G. (2015). Religion and parental 
# cooperation: an empirical test of Slone’s sexual signaling model. In D. Slone & J. Van Slyke (Eds.), 
# The Attraction of Religion: A Sexual Selectionist Account (pp. 29-62). Bloomsbury Press. 

#Joseph Bulbulia
#joseph.bulbulia@vuw.ac.nz
#Religious Studies
#Faculty of Humanities and Social Sciences
#Victoria University
#PO Box 600
#Wellington, NZ, 6140


========================================================

29 Zaara St Newcastle NSW 2300
EthL2.Euro.T1


#################
### PACKAGES ####
##################

#install.packages("glmmADMB",repos=c("http://glmmadmb.r-forge.r-project.org/repos", getOption("repos")),type="source")
#install.packages("coefplot2",repos="http://www.math.mcmaster.ca/bolker/R",type="source")
#library("devtools"); install_github("lme4",user="lme4")
#install.packages("lubridate")
library(car)
library(coefplot2)
library(corrgram)
library(data.table)
library(ggplot2)
library(gmodels)
library(gridExtra)
require(Hmisc)
require(knitr)
library(mice)
library(MCMCglmm)
library(MCMCpack)
library(plyr)
require(xtable)


###################################
#### Notes on citing R packages
##################################

# Select packages to cite:
#citPkgs <- names(sessionInfo()$otherPkgs)
# Write the bibtex file:
#write_bib(citPkgs, file="R-Pckgs")

###############################
# Information about the NZAVS
##############################

# TITLE:
# The New Zealand Attitudes and Values Study 
# Supporting Information for Research Collaboration 
# NZAVS 2012 Time 3.5 Data
# 
# This is the base syntax file for importing and analyzing
# anonymized data from the 2012 (Time 3.5) New Zealand Attitudes
# and Values Study in Mplus. This is designed to complement the 
# SPSS Base Dataset containing the same variables. For a full 
# list of the variables and item labels contained in this wave
# of the study please visit the NZAVS website or contact me directly.
# 
# The Time 3.5 NZAVS survey was designed to provide more in-depth
# measures about specific constructs in addition to those tracked 
# each year in the full NZAVS questionnaire. These measures were 
# assessed only at this mid-year time point, and participants from
# different demographic groups completed different version of this 
# questionnaire. People who affiliated with a religion completed 
# additional questions about their religious attitudes and beliefs.
# Maori completed a version of the questionnaire containing the 
# MMM-ICE. Pacific peoples completed a version of the questionnaire
# containing the PIWBS. These data should not be considered 
# representative and were administered using an online survey. 
# 
# Updated February 2012
# 
# http://www.psych.auckland.ac.nz/uoa/NZAVS
# Dr. Chris Sibley, c.sibley@auckland.ac.nz
# 




```{r} 


## Read in data from machine, recode missing values
system.time(
  longi4 <- as.data.frame(read.table("/Users/josephbulbulia/Dropbox/0R/DATA/NZAVSdata/NZLong/Data/NZAVSLongitudinalSPSSData.csv",
                                     header=TRUE,
                                     sep=",",
                                     na = -9999)))
                                     



##################################################################################################
# Note, the coding of these data is explicit and a bit inelegant: I've aimed for clarity.
# After model five below, there is a record of all the analyses that were done, including those
# not reported in the article. There are no surprises there and nothing that would challenge our 
# central (and rather modest) inferences in the published chapter.
##################################################################################################


l <- as.data.frame(longi4)

# Get territory
l$Territorial.Authority.T3 <- as.factor(as.character(l$Territorial.Authority.T3))

# Remove any dataset with this name.
rm(sgot)

# Create new data set for id variables
sgot <- data.frame(longi4)

# need to count the rows to recover observations
ssgot$Id <- 1:nrow(ssgot)


ssgot$Urban <- as.factor(as.character(l$Territorial.Authority.T3))

# Check data
head(ssgot)


## Now we write a function for coding urban as 0 if rural and 1 if urban


ssgot$Urban <- ifelse(   ssgot$Urban ==  5 |
                         ssgot$Urban ==  6 |
                         ssgot$Urban ==  7 |
                         ssgot$Urban ==  8 |
                         ssgot$Urban == 16 |
                         ssgot$Urban == 23 |
                         ssgot$Urban == 31 |
                         ssgot$Urban == 40 |
                         ssgot$Urban == 44 |
                         ssgot$Urban == 45 | 
                         ssgot$Urban == 46 |
                         ssgot$Urban == 47 |
                         ssgot$Urban == 52 |
                         ssgot$Urban == 60 | 
                         ssgot$Urban == 71, 
                       "Urban", "Rural")


## Make sure it is factor
ssgot$Urban <- as.factor(ssgot$Urban)
str(ssgot$Urban)


## now add religious id for T3
ssgot$rels3 <- as.factor(l$Religious.T3)
table(ssgot$rels3)
str(sgot$rels3)

## now add partner for T3
ssgot$partner3 <- as.factor(l$Partner.T3)

#Check 
table(ssgot$partner3)
str(ssgot$partner3)

# Employment Time 3
ssgot$employ <- as.factor(l$Employed.T3)
table(ssgot$employ)
str(ssgot$employ)

# Income Time 3
ssgot$income3 <- l$Household.INC.T3



# Gender Time 3
ssgot$gen3 <- as.factor(l$Gender.T3)
table(l$gen3)


# Children Time 3
ssgot$children3 <-as.numeric(l$ChildrenNum.T3)
ssgot$children3c <- scale(l$ChildrenNum.T3, scale=FALSE, center=TRUE)
table(ssgot$children3)
#hist(sgot$children3)

#  Children at home Time 3
ssgot$childrenHome3 <-as.numeric(l$ChildrenHome.T3)
table(sgot$childrenHome3)

#  Age Time 3

ssgot$age3 <- as.numeric(l$Age.T3)
summary(ssgot$age3)


# Ethnicity Time 3
ssgot$ethn3 <- as.factor(l$EthL3.NZEuro.T3)
table(ssgot$ethn3)


# Politics Time 3 (scale and center)
ssgot$pol3 <- as.numeric(l$Pol.Orient.T3)
table(ssgot$pol3)
#hist(sgot$pol3)
ssgot$pol3c <- scale(ssgot$pol3, center=TRUE, scale=TRUE)

###################################
###################################
###################################

#TIME 3.5


library(car)

# Earlier the data had a mistaken 2 coded, this should be sorted but if not:

ssgot$rels35 <- recode(l$Religious.T35, "2 = 1")
table(ssgot$rels35)
str(ssgot$rels35)


# Check how many NAs
p <- 2782 + 1675  
4514 - p 

# there are 57 NAs


#create gender variable 3.5 
ssgot$gen35 <- l$Gender.T35
ssgot$gen3 <- as.numeric(l$Gender.T35)


# Check
table(ssgot$gen35)
str(ssgot$gen35[110:120])
ssgot$gen35<-factor(ssgot$gen35)

# Time 3.5. denomination large
ssgot$denom<- as.factor(l$Religion.Denoms.T35)
table(ssgot$denom)

# Time 3.5 age
ssgot$age35 <- as.numeric(as.character(l$Age.T35))
summary(ssgot$age35)
ssgot$age35c <- scale(ssgot$age35, center=TRUE, scale=FALSE)


##############################################
##Make data frame 
##############################################

dim(ssgot)
head(ssgot)

ssgot$lincome3c  <-  scale(log((l$Household.INC.T3)+1), center = TRUE, scale = TRUE)
ssgot$incomeoverten <- (l$Household.INC.T3)/(10000)


ssgot$rels35 <- factor(l$Religious.T35)
str(ssgot$rels35)
ssgot$rels35[1:100]

ssgot$edu1 <- as.factor(l$EduOrdinal.T1)
str(ssgot$edu1)


# church
ssgot$church<- l$Religion.Church.T35 
mean(ssgot$church, na.rm=TRUE)


sd(ssgot$church, na.rm=TRUE)
table(ssgot$church)

## Scripture
table(l$Religion.Scripture.T35)


# Prayer
ssgot$pray <-l$Religion.Prayer.T35
mean(ssgot$pray , na.rm=TRUE)
sd(ssgot$pray , na.rm=TRUE)
#boxplot(ssgot$pray)


## Center
ssgot$prayC <- scale(ssgot$pray, center=TRUE, scale=FALSE)
ssgot$churchC <- scale(ssgot$church, center=TRUE, scale=TRUE)
summary(ssgot$edu1)


## NZ DEP
ssgot$NZDEP3 <- l$NZDep.T3
ssgot$NZDep3C <- scale(ssgot$NZDEP3, center =TRUE, scale =FALSE)
summary(ssgot$NZDep3C)


library(data.table)

# scripture

ssgot$script <- l$Religion.Scripture.T35
mean(ssgot$script, na.rm=TRUE)
table(ssgot$script)
hist(ssgot$script)

## Large denoms

# Time3.5 denominations censu  (blows up models!)
ssgot$ldenom<- as.factor(l$Religion.L4.Cats.T35)
table(ssgot$ldenom)

str(ssgot$ldenom)

## Small denoms
sgot$denom<- as.factor(l$Religion.Denoms.T35)
table(sgot$denom)


# Religious Identification 

ssgot$rimp3 <- l$Religion.Identification.T35
summary(ssgot$rimp3)
ssgot$rimp3C <- scale(ssgot$rimp3, center = TRUE, scale = FALSE)
ssgot$rimp35 <- l$Religion.Identification.T35
ssgot$rimp35C <- scale(ssgot$rimp35, center = TRUE, scale = FALSE)
summary(ssgot$rimp35C)


## Belief in God

ssgot$belief <- l$Believe.God.T35
summary(ssgot$belief)
#hist(ssgot$belief)
ssgot$beliefS <- scale(l$Believe.God.T35, center=TRUE, scale=TRUE)

# No longer exists
#head(l$REL.COLLNARC.T35)

# Create Religious ingroup ID variable
ta<- cbind(l$RelCollectNarc01.T35, l$RelCollectNarc02.T35, l$RelCollectNarc03.T35)


# Check

head(ta)
rmean <-rowMeans(ta, dims=1)
head(rmean)
ssgot$relnar <- as.numeric(rmean)# 
ssgot$relnar 
head(ssgot$relnar)

ssgot$relnarc <- scale(ssgot$relnar, center=TRUE, scale=TRUE)
summary(ssgot$relnar)

# Parents religious
ssgot$parRel35 <- as.factor(l$ParentsReligious.T35)
ssgot$gen35 <-  as.factor(l$Gender.T35)
ssgot$gen3 <-  as.factor(l$Gender.T3)
ssgot$gen35S <-  as.factor(l$Gender.T35)



## Contrast code ethn  (For some reason these get stripped later, so we have to put them back)
head(ssgot)
contrasts(ssgot$gen35S) <- contr.sum(2)


## Contrast code ethn  (For some reason these get stripped later, so we have to put them back)
contrasts(ssgot$ethn3) <- contr.sum(2)
str(ssgot$gen35S)


#########  Key subset   RELIGION ONLY GORUP  #####################
head(ssgot)

relgot <-subset(ssgot, rels35 == 1)
head(relgot)

######################################################################



head(relgot)



#########################################
# descriptive 
#########################################
head(ssgot)

library(xtable)
library(gmodels)

mytable <- table(ssgot$rels35, ssgot$gen35)
mytable
#Prop.Table(mytable)
prop.test(mytable)




tabs <- with(ssgot,CrossTable(children3, Urban))
table(tabs)
mytable<- table(ssgot$children3,ssgot$Urban)
dim(ssgot)


prop.test(mytable)




########################################
## Initial exploration of the data
### Check differences for urban v. rural in expected children, first in the general population

summary(glm(data=ssgot, children3 ~ factor(Urban), na.action=na.omit,family="quasipoisson"))

## Next in the religious subsample
summary(glm(data=relgot, children3 ~ factor(Urban), na.action=na.omit, family ="quasipoisson"))


## We use a poisson distribution  because this is a rate parameter, 
## and we exponentiate to recover the data scale.

#General sample
#Rural general
exp(.84220)

## Urban
exp(.84220-.25275)

## Difference
exp(.84220)- exp(.84220-.25275)

#Religion religion subset
## Urban
exp(.89785)
## Rural
exp(.89785-.23068)

## Difference
exp(.89785)-
exp(.89785-.23068)

chisq.test(mytable)
#	Pearson's Chi-squared test with Yates' continuity correction
#data:  mytable
#X-squared = 10.3172, df = 1, p-value = 0.001318

xtabs(~rels35 + gen35, data = ssgot)
table(relgot$partner)
head(ssgot)


### Note about employment and education (we do not generally use these varaibles in the NZAVS because it is not as inforamtive as Deprivation)

## Exploration of the data
str(relgot)
library(xtable)
library(gmodels)


summary(ssgot$ethn3)
with(relgot, CrossTable(gen35, children3))
with(relgot, CrossTable(employ, gen35))

with(relgot, CrossTable(church, gen35))
plot <- with(relgot, CrossTable(children3, gen35))
boxplot(relgot$gen35, relgot$children3, na.action = na.omit)
head(relgot)



## Summaries
library(ggplot2)
library(plyr)

ddply(relgot, "gen3", summarise, children3.mean=mean(children3, na.rm =TRUE))

ddply(relgot, "partner3", summarise, children3.mean=mean(children3, na.rm =TRUE))

ddply(relgot, "ethn3", summarise, children3.mean=mean(children3, na.rm =TRUE))



## Here is a useful function for extracting datacolumns and excluding missing data

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
%
nrow(relgot)

roo<- completeFun(relgot, c("Urban"))
nrow(relgot)- nrow(roo)



## Again we look at the effects of "urban which appear large"
summary(lm(data=roo, children3 ~ Urban, na.action=na.omit))
summary(aov(data=roo, children3 ~ Urban, na.action=na.omit))

ddply(ssgot, "Urban", summarise, children3.mean=mean(children3, na.rm =TRUE))

ddply(relgot, "Urban", summarise, children3.mean=mean(children3, na.rm =TRUE))



############# Some Exploratory graphs

pr <- ggplot(relgot, aes(x=children3, fill = gen3)) + geom_histogram(binwidth=1, position = "dodge", alpha=.5)  + facet_grid(~gen3)
pr


pr <- ggplot(relgot, aes(x=children3, fill = gen3)) + geom_density(position = "dodge", alpha=.1)   + facet_grid(~gen3)
pr

pr + geom_vline(data=relgottn, aes(xintercept=children3.mean,  colour=gen3),
               linetype="dashed", size=1)

pr <- ggplot(relgot, aes(x=partner3, fill = gen3)) + geom_density(position = "dodge", alpha=.1) 
pr + geom_vline(data=relgottn, aes(xintercept=children3.mean,  colour=gen3),
                linetype="dashed", size=1)


head(ssgot)
sr <- ggplot(ssgot, aes(x=children3, fill = rels3)) + geom_histogram(binwidth=1, position = "dodge", alpha=.5, frequency =TRUE)
sr + facet_grid(~gen3)
 ggplot(ssgot, aes(x=children3, fill = rels3)) + geom_histogram(position = "dodge", alpha=.5, frequency =TRUE)
 ggplot(ssgot, aes(x=rels3, y = children3, fill = rels3)) + geom_boxplot(notch=TRUE, na.rm =TRUE)


head(relgot)


#####################
## More exploration
#####################


library(ggplot2)
ggplot(relgot, aes(x=children3)) + geom_density()

hist(men$children3, freq = TRUE)
with(ssgot, CrossTable(rels35, gen35E))


with(ssgot, CrossTable(children3,rels35))
with(ssgot, chisq.test(children3,rels35))
with(roo, chisq.test(Urban,children3))


numb <- subset(relgot, children3 != 0 & gen35 == 1 & partner3 == 1)
nrow(numb)

       head(relgot)


head(relgot)

with(ssgot, CrossTable(ethn3, rels35))
with(ssgot, (ethn3, rels35))


mean(relgot$children3, na.rm=TRUE)
mean(ssgot$children3, na.rm=TRUE)


library(gmodels)
with(relgot, CrossTable(church, gen35))



### Political orientation, summaries

with(relgot, CrossTable(pol3, rels35))
mean(ssgot$pol3, na.rm=TRUE)
sd(ssgot$pol3, na.rm=TRUE)
mean(relgot$pol3, na.rm=TRUE)
sd(relgot$pol3, na.rm=TRUE)

##################################

head(relgot)

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}



head(relgot)
telgot <-  completeFun(relgot, c("gen35S","ldenom")
)
nrow(telgot)
head(telgot)

## Code Denom
str(telgot$Urban)
contrasts(telgot$Urban) <- contr.sum(2)


## Now we get data for imputing missing values. Note: will not impute "children" or "children at home". Children is an outcome variable and missingess can be estimated during MCMC. Children at Home is a random effect, and will be estiamted during MCMC. 


nrow(telgot)
head(telgot)
amdata <- telgot[ , c("Id", 
                         "age35",
                         "gen35S", 
                         "ethn3", 
                         "partner3",
                         "pol3", 
                         "rimp35",
                         "children3",
                         "belief",
                         "relnar", 
                         #"income3",
                         "pray",
                         "church",
                         #"denom",
                         "ldenom",
                          "NZDEP3",
                          "Urban",
                          "childrenHome3"
                         # "script"
                         )]

head(amdata)
nrow(amdata)
str(amdata)



library(mice)
head(amdata)
str(amdata)

# leave out Children as these values will be imputed during MCMC
head(amdata)
dim(amdata)
# We do not estimate children because it is a key variable of theortetical interest
mi <- mice(amdata[,-c(8,16)])
#head(mi)


#plot(mi)
now<- complete(mi)
head(now)
#add children back (we do not estimate, but let MCMC do this)
head(amdata)
cnow <- cbind(now[],amdata$childrenHome3,amdata$children3)
head(cnow)

head(cnow)
nrow(cnow)

## Time to get some tables for the data (including the newly imputed data).

library(stargazer)
head(now)
starnow <- cnow[,c(2,6,7,8,9,10,11,13,16)]
head(starnow)
colnames(starnow) <- c("Age","Conservativism", "Religious Identification","Belief", "Religious In-group Values", "Pray","Church", "NZDep","Children")
stargazer(starnow,summary=TRUE, omit.summary.stat=c("max","min"), type="text")

stargazer(starnow,summary=TRUE, omit.summary.stat=c("max","min","n"))
str(cnow)

head(cnow)
names(cnow)[16] <- "children"
names(cnow)[15] <- "childrenHome"
names(cnow)[3] <- "gender"

head(cnow)
ddply(cnow, "gender", summarise, children.mean=mean(children, na.rm =TRUE))

ddply(cnow, "partner3", summarise, children.mean=mean(children, na.rm =TRUE))

ddply(cnow, "ethn3", summarise, children.mean=mean(children, na.rm =TRUE))

ddply(cnow, "gender", summarise, partner3.mean=mean(childrenHome, na.rm =TRUE))


library(gmodels)

with(cnow, CrossTable(ethn3))
with(cnow, CrossTable(children, gender))


### Pols
with(relgot, CrossTable(pol3, rels35))


#######################################################################################################################################
# cor matrix: Next we need a correlation matrix. Note that in the publication the only correlation matrix is the simple table

########################################################################################################################################

######## Imputation



#install.packages("lubridate")

require(Hmisc)
require(xtable)


corstarsl <- function(x){ 
  require(Hmisc) 
  x <- as.matrix(x) 
  R <- rcorr(x)$r 
  p <- rcorr(x)$P 
  
  ## define notions for significance levels; spacing is important.
  mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " ")))
  
  ## trunctuate the matrix that holds the correlations to two decimal
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1] 
  
  ## build a new matrix that includes the correlations with their apropriate stars
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x)) 
  diag(Rnew) <- paste(diag(R), " ", sep="") 
  rownames(Rnew) <- colnames(x) 
  colnames(Rnew) <- paste(colnames(x), "", sep="") 
  
  
  ## remove upper triangle
  Rnew <- as.matrix(Rnew)
  Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
  Rnew <- as.data.frame(Rnew) 
  
  ## remove last column and return the matrix (which is now a data frame)
  Rnew <- cbind(Rnew[1:length(Rnew)-1])
  return(Rnew) 
}

head(cnow)
cornow <- cnow
head(cornow)
colnames(cornow) <- c("id", "Age","Gender","Eth","Partn","Consrv","RelImp","Belief",
                      "RelGrp","Pray","Church","Denom","NZDep", "Urban","ChildHome","Child")
head(cornow)


forcor <- cornow[,c(2,6:11,13)]
head(forcor)
str(forcor)
sforcor <- scale(forcor)
str(sforcor)
head(sforcor)
cor(sforcor)

############################################
##  USE
############################################
xtable(corstarsl(sforcor))
########################################

## Here we check variance inflation factors. The script below is from:

# 
# @book{zuur2013beginner,
#   Author = {Zuur, AF and Hilbe, JM and Ieno, EN},
# 	Publisher = {Highland Statistics},
# 	Series = {Newburgh, UK},
# 	Title = {A beginner's guide to GLM and GLMM with R},
# 	Year = {2013}}
# They should be cited.




corvif <- function(dataz) {
  dataz <- as.data.frame(dataz)
  #correlation part
  #cat("Correlations of the variables\n\n")
  #tmp_cor <- cor(dataz,use="complete.obs")
  #print(tmp_cor)
  #vif part
  form    <- formula(paste("fooy ~ ",paste(strsplit(names(dataz)," "),collapse=" + ")))
  dataz   <- data.frame(fooy=1,dataz)
  lm_mod  <- lm(form,dataz)
  
  cat("\n\nVariance inflation factors\n\n")
  print(myvif(lm_mod))
}



#Support function for corvif. Will not be called by the user
myvif <- function(mod) {
  v <- vcov(mod)
  assign <- attributes(model.matrix(mod))$assign
  if (names(coefficients(mod)[1]) == "(Intercept)") {
    v <- v[-1, -1]
    assign <- assign[-1]
  } else warning("No intercept: vifs may not be sensible.")
  terms <- labels(terms(mod))
  n.terms <- length(terms)
  if (n.terms < 2) stop("The model contains fewer than 2 terms")
  if (length(assign) > dim(v)[1] ) {
    diag(tmp_cor)<-0
    if (any(tmp_cor==1.0)){
      return("Sample size is too small, 100% collinearity is present")
    } else {
      return("Sample size is too small")
    }
  }
  R <- cov2cor(v)
  detR <- det(R)
  result <- matrix(0, n.terms, 3)
  rownames(result) <- terms
  colnames(result) <- c("GVIF", "Df", "GVIF^(1/2Df)")
  for (term in 1:n.terms) {
    subs <- which(assign == term)
    result[term, 1] <- det(as.matrix(R[subs, subs])) * det(as.matrix(R[-subs, -subs])) / detR
    result[term, 2] <- length(subs)
  }
  if (all(result[, 2] == 1)) {
    result <- data.frame(GVIF=result[, 1])
  } else {
    result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
  }
  invisible(result)
}
#END VIF FUNCTIONS
# Exlcude NA's

stargazer(test)
# Create Myvars 
MyVar <- c("Age","Gender", "Eth", "Partn", "Consrv", "RelImp", "Belief", "RelGrp", "Pray", "Church", "NZDep","Urban")
dim(now)

head(now)
cordat <- now
head(cordat)
corvarb <- cordat[,c(2,6,7,8,9,10,11,13)]
head(corvarb)
colnames(corvarb) <- c("Age",  "Consrv", "RelImp", "Belief", "RelGrp", "Pray", "Church", "NZDep")

MyVar <- c("Age",  "Consrv", "RelImp", "Belief", "RelGrp", "Pray", "Church", "NZDep")


  
###USE
stargazer(corvarb, summary=TRUE, title="VIF test: inflation of predictors", type="text")


# Check VIF

stargazer(corvif(corvarb[,MyVar]), summary =FALSE,title="VIF test: inflation of predictors")

library(stargazer)
head(cornow)
starnow <- cornow[,c(2,6,7,8,9,10,11,13,16)]
head(starnow)
names <- c("Age","Conservativism","Belief", "Rel_Group_Vals","Pray","Church","NZDep","Children")
colnames(starnow) <- c("Age","Conservativism","Religious_Ident","Belief", "Rel_Group_Vals","Pray","Church","NZDep","Children")
stargazer(starnow,summary=TRUE, omit.summary.stat=c("max","min","n"), type="text")

corvif(starnow[,])




# 
#            GVIF
# Age    1.160948
# Gender 1.110915
# Eth    1.019747
# Partn  1.156299
# Consrv 1.099262
# RelImp 1.782577
# Belief 1.380340
# RelGrp 1.359308
# Pray   1.172733
# Church 1.258573
# NZDep  1.071983
# Urban  1.034901
# 






######################################################################################
# Is a multi-level model needed?  We can assess this queston using the ICC
#######################################################################################
check if multilevel model needed
library(multilevel)
library(car)

head(cornow)
testml <- aov(Child ~ as.factor(Denom), data=cornow)
testm2 <- aov(Pray ~ as.factor(Denom), data=cornow)
testm3 <- aov(Church ~ as.factor(Denom), data=cornow)
testm4 <- aov(Belief ~ as.factor(Denom), data=cornow)
testm5 <- aov(Consrv ~ as.factor(Denom), data=cornow)
testm6 <- aov(RelGrp ~ as.factor(Denom), data=cornow)



round(ICC1(testml), 2) # .05   how much of the variance is explained by group membership
round(ICC2(testml), 2) #  .63  0.6318678  group mean reliability


round(ICC1(testm2), 2)  # .07 how much of the variance is explained by group membership
round(ICC2(testm2) , 2)  #  .77 group mean reliability

round(ICC1(testm3) , 2) # .11 how much of the variance is explained by group membership
round(ICC2(testm3) , 2)  # .84  group mean reliability

round(ICC1(testm4), 2)  #  .22 how much of the variance is explained by group membership
round(ICC2(testm4)  , 2) # .92  group mean reliability


round(ICC1(testm5), 2)  #  .9 how much of the variance is explained by group membership
round(ICC2(testm5)  , 2) # .82  group mean reliability


round(ICC1(testm6), 2)  #  .1 how much of the variance is explained by group membership
round(ICC2(testm6)  , 2) # .82  group mean reliability

stargazer(cbind(ta,ma), summary=F)

head(cornow)
corn <- cornow[,c(2,6:11,13)]
head(corn)
# 
# library(elipse)
# R = cor(corn)
# R <- round(R, 3)
# 
# plotcorr(R,col = colorRampPalette(c("firebrick3", "white", "navy"))(10))
# plotcorr(R, col = colorRampPalette(c("#E08214", "white", "#8073AC"))(10), type = "lower")

corrn <- corn[,-c(3)]

head(cnow)


### Correlation matrices (pretty but not used)
library(corrgram)
corrgram(cornow, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.pie, text.panel=panel.txt,
         main="Correlation of Modelled Variables")
%

corrgram(corn, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.shade, text.panel=panel.txt, dir="left",
         main="Correlation of Religion Variables")



corrgram(cornow, order=NULL, lower.panel=panel.shade,
         upper.panel=NULL, text.panel=panel.txt,
         main="Religion Variables (unsorted)")

corrgram(cornow, order=TRUE, lower.panel=panel.ellipse,
         upper.panel=panel.pts, text.panel=panel.txt,
         diag.panel=panel.minmax, 
         main="Religion Variables")


#### TEST
corrgram(cornow, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.pie, text.panel=panel.txt,
         main="Correlation of Religion Variables")

# 
 head(cornow)
colnames(cornow) <- c("Id", "Age","Gender","Eth","Prtnr","Consrv","RelImp","Belief",
                        "RelGrp" ,"Pray" ,"Church","Denom","NZDep","Urban","ChildHome", "Child")
head(cornow)
str(cornow)
%
#scale varaibles


## Note, NZAVS guidlines -- generally we don't use Income.  Education not needed
cornow$AgeC <- as.integer(scale(cornow$Age, center = TRUE, scale = FALSE))

## We want to center zero a 20, to get the full range fertility 

## 



# Get standard deviation of agey, we want to use a more sensible zero that the sample mean
mage<- mean(cornow$Age)
sdage <- sd(cornow$Age)

# mage-sdage  New Centered age is 34 years old



cornow$AgeCy <-  cornow$AgeC + sdage
mean(cornow$AgeCy)  # 15.8361
  
str(cornow)


### Now we make sure all the data are in shape for modelling.


contrasts(cornow$Urban) <- contr.sum(2)
contrasts(cornow$Urban)

cornow$PrayC <- as.integer(scale(cornow$Pray, center = TRUE, scale = FALSE))
cornow$ChurchC <- as.integer(scale(cornow$Church, center = TRUE, scale = FALSE))
cornow$RelImpC <- as.numeric(scale(cornow$RelImp, center = TRUE, scale = FALSE))
cornow$BeliefC <- scale(cornow$Belief, center = TRUE, scale = FALSE)
cornow$ConsrvC <- scale(cornow$Consrv, center =TRUE, scale = FALSE)
cornow$NZDepC <- scale(cornow$NZDep, center =TRUE, scale = FALSE)
cornow$RelGrpC <- scale(cornow$RelGrp, center =TRUE, scale = FALSE)

# Scales as well as center to improve MCMC Mixing

cornow$PrayS <- (scale(cornow$Pray, center = TRUE, scale = TRUE))
cornow$ChurchS <- (scale(cornow$Church, center = TRUE, scale = TRUE))
cornow$RelImpS <- scale(cornow$RelImp, center = TRUE, scale = TRUE)
cornow$BeliefS <- scale(cornow$Belief, center = TRUE, scale = TRUE)

### We don't use reprot these exogenous variables -- though when trying them, they make no difference to inference(and not significant.)
#cornow$Edu <- as.factor(cornow$Edu)
#cornow$InC <- scale(cornow$Inc, center =TRUE, scale = FALSE)
#cornow$ChildC <- scale(cornow$Children, center =TRUE, scale = FALSE)  # Note we estiamte children as multivariate normal during MCMC so no need to impute using a non-bayesian algorthim. 

cornow$ConservS <- scale(cornow$Consrv, center =TRUE, scale = TRUE)
cornow$NZDepS <- scale(cornow$NZDep, center =TRUE, scale = TRUE)
cornow$RelGrpS <- scale(cornow$RelGrp, center =TRUE, scale = TRUE)
cornow$ConsrvS <- scale(cornow$Consrv, center =TRUE, scale = TRUE)

head(cornow)

contrasts(cornow$Prtnr) <- contr.sum(2)
contrasts(cornow$Eth) <- contr.sum(2)
contrasts(cornow$Urban) <- contr.sum(2)
contrasts(cornow$Gender) <- contr.sum(2)
cornow$Id <- as.factor(cornow$Id)


##  Do correlation plot
library(corrgram)
corrgram(corn, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.pie, text.panel=panel.txt,
         main="Correlation of Modelled Variables")
#
corrgram(corn, order=TRUE, lower.panel=panel.ellipse,
upper.panel=panel.pts, text.panel=panel.txt,
diag.panel=panel.minmax, main="Correlation of Modelled Variables")




###########################
### Bayesian estimation 
###########################




## P.108 Jarrod Hadfield Course Notes:  This is the key methods bit
# We can then test for zero-flation by constraining the over-dispersion to be the same 
# for both process using a trait by units interaction in the R-structure, and by setting up 
# the contrasts so that the zero-altered regression coecients are expressed as di↵erences 
# from the Poisson regression coecients. 
# When this di↵erence is zero the variable causes no zero-flation, when it is negative it 
# causes zero-inflation and when it is positive it causes zero-deflation:


##p.101 Setting V=diag(2) and nu=0.002 we have the inverse-gamma prior with shape=scale=0.001 #for the residual component of the Poisson pro- cess which captures over-dispersion:

library(MCMCglmm)
library(coefplot2)

### make correlation matrix
#correlation.matrix <- cor(attitude[,c("rating","complaints","privileges")])
head(data)
require(Hmisc) 


this <- data[,c("Age","Consrv","Belief","RelGrp","Pray","Church","NZDep","Child","ChildHome")]
this <- na.exclude(this)

# Function defined above
xtable(corstarsl(this))
#stargazer(this, title="Correlation Matrix of Continuous Model Variables", type="text", align=FALSE, summary=FALSE)


## Note that in the end we did not use freq estimators, in keeping with the NZAVS 14 April 2014 protocols. 





# Check data
head(cornow)

# Fix bug in MCMCglmm
mod$Residual.nrt <- 2




# Note this is overlill, just doesn't help with mixing 
# prior5 =list(B= list(mu = c(rep(0,10),rep(0,10)), V= diag(20)*100), #*(1 + pi^2/3), 
#                            R = list(V = diag(2), nu = 2.002, fix =2),
#                            G = list( G1 = list(V = diag(2)*.02, nu = 3)))


# #  P11
# MCMCglmm uses an inverse Wishart prior for the (co)variances and a normal prior for the fixed e↵ects. 
# 
# inverse Wishart takes two scalar param- eters, V and nu. The distribution tends to a point mass on V as the degree of belief parameter, nu goes to infinity. The distribution tends to be right skewed when nu is not very large, with a mode of V⇤nu but a mean of V⇤nu (which is not defined for nu < 2).1 nu+2
# # 
# NOt a lot of information in the data for the marginal distribution for children, hence we use stronger priors
##
##
##  prior1.1<-list(G=list(
#+            G1=list(V=matrix(p.var/2),n=1)),
#+            R=list(V=matrix(p.var/2),n=1))

data  <- cornow
nittA <- 55000
nittC <- 505000
burnin <- 5000
thin <- 10




############################################
## Models
############################################

####################################################################################
## rel ingroup and prayer, what predicts it? 
####################################################################################

# Prior for Poiss, Zap, and Gauss modesl  # We use "~us variance to model how belief and in-group commiment might covary"
priorPsych = list( 
             R = list(V = diag(2)*.02, nu = 3),
             G = list( G1 = list(V = diag(2) *.02, nu = 3)))


# Note we tested censored regression because religious belief looks to be subject to ceiling effects (raw data)# Makes no diff.
# modtest <-  MCMCglmm(RelImp ~ 1, data=data)
# head(data, 18)
# head(data)
# data$Beliefmax <- data$Belief #ifelse(data$Belief <= 7, 7, data$Belief)
# data$Beliefmin <- rep(1, nrow(data)) #ifelse(data$Belief <= 1, 1, data$Belief)
# head(data)
# hist(data$Belief)

                    
###We use different

# No evidence of covariation so we fit the simpler IDH model
mod1 <-  MCMCglmm(cbind(cbind(Beliefmin, Beliefmax),  RelGrp)  ~ 
                    trait -1  + 
                    trait:PrayS + 
                    trait:ChurchS +  
                    trait:Gender +   
                    trait:Prtnr +  
                    trait:AgeC  +   
                    trait:Eth +  
                    trait:ConsrvS + 
                    trait:NZDepS +   
                    trait:Urban,   
                  random = ~us(trait):Denom, 
                  family=c("cengaussian","gaussian"),
                  rcov = ~ us(trait):units,  
                  prior=priorPsych,
                  nitt=nittA,
                #  pl=TRUE,
                #  pr=TRUE,
                  data=data)


# Prior for Poiss, Zap, and Gauss modesl  # We use "~us variance to model how belief and in-group commiment might covary"
# 
# 
# testmod1 <-  MCMCglmm(cbind(cbind(BeliefCenpoint, Belief))  ~ PrayS +
#                         ChurchS +  
#                         Gender +   
#                         Prtnr +  
#                         AgeC  +   
#                         Eth +  
#                         ConsrvS + 
#                         NZDepS +   
#                         Urban,   
#                         random = ~Denom, 
#                   family=c("cengaussian"),
#                   rcov = ~ units,  
#               #    prior=priorPsych,
#                   nitt=nitt,
#                   pl=TRUE,
#                   pr=TRUE,
#                   data=data)
# 
# plot(testmod1)

#### FIT EASIER MODEL TO SEE THAT cengauss does not effect inferene

  

mod1_simp <-  MCMCglmm(cbind(Belief,  RelGrp)  ~ 
                    trait -1  + 
                    trait:PrayS + 
                    trait:ChurchS +  
                    trait:Gender +   
                    trait:Prtnr +  
                    trait:AgeC  +   
                    trait:Eth +  
                    trait:ConsrvS + 
                    trait:NZDepS +   
                    trait:Urban,   
                  random = ~us(trait):Denom, 
                  family= rep("gaussian",2),
                  rcov = ~ us(trait):units,  
                  prior=priorPsych,
                  nitt=nittA,
                # pl=TRUE,
                # pr=TRUE,
                  data=data)



summary(mod1_simp)


summary(mod1_simp)

coefplot2(mod1_simp$Sol, intercept=TRUE)
coefplot2(mod1$VCV, var.inx=c(1:3))

plot(mod1$Sol)

  
abbfun <- function(x) {
   gsub("(Sol\\.)*(trait|at.level\\(trait, 1\\):)*","",
        gsub("trait_za","za",x))
}   

vn1 <- abbfun(rownames(coeftab(mod1)))
vnn1 <- vn1[1:20]
# insrt<- rep(c("Prayer(S)","Church(S)","Gender","Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban"),2)

vars1  <- c(("Intercept:Belief"), "Intercept:Grp", 
            "Bel:Prayer(S)",
            "Grp:Prayer(S)",
            "Bel:Church(S)",
            "Grp:Church(S)",
            "Bel:Gender",
            "Grp:Gender",
            "Bel:Partner",
            "Grp:Partner",
            "Bel:Age(C)",
            "Grp:Age(C)",
            "Bel:NZ Euro",
            "Grp:NZ Euro",
            "Bel:Conservative(S)",
            "Grp:Conservative(S)",
            "Bel:NZDep(S)",
            "GrpNZDep(S)",
            "Bel:Urban",
            "Grp:Urban")




# varss  <- c( "Intercept:Grp" ,"Prayer (S)","Church (S)","Gender","Partner","Age (C)","NZ Euro", "Conservative (S)","NZDep (S)","Urban","Z_Prayer (S)","Z_Church (S)","Z_Gender","Z_Partner","Z_Age (C)","Z_NZ Euro", "Z_Conservative (S)","Z_NZDep (S)","Z_Urban")


######USEE THIS #############
coefplot2(mod1_simp$Sol,
           intercept =TRUE,
          # merge.names=FALSE,
           varnames=vars1,
          col=c("gray","gray", "red","blue", "red","blue","red","blue", "green","green",rep("gray",10)))

######USEE THIS #############

# NO Diff between cengaussian and Gaussian
coefplot2(list("Gaussian"=mod1_simp$Sol,"CenGaussian"= mod1$Sol),
           intercept =TRUE,
           merge.names=FALSE,
           #varnames=vars,
           legend=TRUE,legend.x="right")


summary(mod1_simp)
coefplot2(mod1_simp)
#



# Diagnostics: non normality (latent normal) but lower dispersion.
 mcmc.units <- colMeans(mod1$Liab) - predict(mod1,   marginal = NULL, type = "terms")
 par(mfrow=c(1,2))
 qqnorm(mcmc.units) ; qqline(mcmc.units)
 hist(mcmc.units, col = grey(.6), main = "Over-dispersion")
%
head(mod1_simp$Liab)

# Diagnostics # normality but greater dispersion 
 mcmc.units <- colMeans(mod1_simp$Liab) - predict(mod1_simp,   marginal = NULL, type = "terms")
 par(mfrow=c(1,2))
 qqnorm(mcmc.units) ; qqline(mcmc.units)
 hist(mcmc.units, col = grey(.6), main = "Over-dispersion")
%

# enom  0.51509  0.7521  0.9243  1.1348  1.7038
# RelGrp:Belief.Denom -0.53617 -0.3243 -0.2436 -0.1725 -0.0541
# Belief:RelGrp.Denom -0.53617 -0.3243 -0.2436 -0.1725 -0.0541
# RelGrp:RelGrp.Denom  0.08301  0.1432  0.1879  0.2497  0.4274
# Belief:Belief.units  1.83827  1.9239  1.9715  2.0191  2.1152
# RelGrp:Belief.units  0.30745  0.3625  0.3937  0.4234  0.4830
# Belief:RelGrp.units  0.30745  0.3625  0.3937  0.4234  0.4830
# RelGrp:RelGrp.units  1.49955  1.5707  1.6103  1.6495  1.7267


# 
summary(mod1_simp$Sol)
#Belief effects vary across denominations
mean((mod1_simp$VCV[, 8])/(rowSums(mod1_simp$VCV)) )  # Expected correlation for belief by denominations = 002  

#RelGR/belief covariance 
# #Belief effects vary across denominations
# mean((mod1_simp$VCV[, 1])/(rowSums(mod1$VCV)) )  # Expected correlation for denominations = 09  
# 
# #RelGR effects vary across denominations
# mean((mod1_simp$VCV[, 2])/(rowSums(mod1$VCV)) )  # Expected correlation for denominations = 19 


 
# #########  Extracting means + standard errors of means for the table
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
tab <- summary(mod1_simp)
reg <- as.data.frame(tab$solutions)
reg
rreg <- round(reg,4)
rreg
rreg <- rreg[,-c(4)]
rreg

# mod1se <- mod1se[c(1, 3:11, 2, 12:20),]
# mod1se

# put BELIEF AND GROUP TOGETHER ODD Numbering
col1.1 <- rreg[c(1,3,5,7,9,11,13,15,17,19), ]
col2.1 <- rreg[c(2,4,6,8,10,12,14,16,18,20), ]
col1.1
col2.1


# 
# library(stargazer)
# head(now)
# starnow <- now[,c(2,6,7,8,9,10,13)]
# head(starnow)
# colnames(starnow) <- c("Age","Conservativism","Religious Ident","Rel Group Vals","NZDep")
# stargazer(starnow,summary=TRUE, omit.summary.stat=c("max","min","n"), type="text")






varss  <- c( "Intercept:Grp" ,"Prayer (S)","Church (S)","Gender","Partner","Age (C)","NZ Euro", "Conservative (S)","NZDep (S)","Urban","Z_Prayer (S)","Z_Church (S)","Z_Gender","Z_Partner","Z_Age (C)","Z_NZ Euro", "Z_Conservative (S)","Z_NZDep (S)","Z_Urban")


table1 <- cbind(col1.1,col2.1)
table1 
rownames(table1) <- c("Intercept","Pray(S)","Church(S)","Gender","Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban")
table1


colnames(table1)[1] <- "Belief"
colnames(table1)[5] <- "RelGrp"
table1

star1 <- stargazer(table1,summary=FALSE, type="text")
star1 <- stargazer(cbind(table1),summary=FALSE)




coefplot2(mod1$Sol, var.idx=c(1:20))





###############################################################
## Do Women Pray More than Men 
#################################################################
priorPray = list( 
             R = list(V = diag(1)*.02, nu = 2),
             G = list( G1 = list(V = diag(1) *.02, nu =3),
                       G2 = list(V = diag(1) *.02, nu= 3)))

mod2 <-  MCMCglmm(Pray ~  trait *(Gender +
                                    BeliefS + 
                                    Prtnr + 
                    AgeC +  
                    Eth + 
                    ConsrvS + 
                    NZDepS +  
                    Urban), 
                  random = ~ trait:Denom + trait:ChildHome,
                  family = "zapoisson",
                  rcov = ~ trait:units, 
                  prior=priorPray ,
                  nitt=nittA,
                  thin=thin,
                  burnin=burnin,
                #  pr = TRUE,
               #   pl = TRUE,
               #  saveX =TRUE,
               #  saveZ = TRUE,
                  data=data)




## ICC  family/family + units


summary(mod2)



mean((mod2$VCV[, 1])/(rowSums(mod2$VCV)) )  # Expected correlation for denominations = .11
mean((mod2$VCV[, 3]/(rowSums(mod2$VCV)) ))  # Expected correlation for denominations = .02

# 
#  names2<- c(("Intercept"),"ZA_Intercept","Men","Belief(S)", "Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban", "Men X Belief(S)","ZA_Men","ZA_Belief(S)", "ZA_Partner","ZA_Age(C)","ZA_NZ Euro", "ZA_Conservative(S)","ZA_NZDep(S)","ZA_Urban", "ZA_Men X Belief(S)")


names2<- c(("Intercept"),"ZA_Intercept","Men","Belief(S)", "Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban","ZA_Men","ZA_Belief(S)", "ZA_Partner","ZA_Age(C)","ZA_NZ Euro", "ZA_Conservative(S)","ZA_NZDep(S)","ZA_Urban")

coefplot2(mod2$Sol, intercept = TRUE,
          varnames= names2,
          col=c("gray","gray","red","red", rep("gray",6),"red","red", rep("gray",6)))
)

#
coefplot2(mod2,  var.idx=c(1:18))
coefplot2(mod2, var.idx=c(1:18))
summary(mod2)
plot(mod2)
#exp(.35)


# Diagnostics
 mcmc.units <- colMeans(mod2$Liab) - predict(mod2,   marginal = NULL, type = "terms")
 par(mfrow=c(1,2))
 qqnorm(mcmc.units) ; qqline(mcmc.units)
 hist(mcmc.units, col = grey(.6), main = "Over-dispersion")
%

 
# #########  Extracting means + standard errors of means for the table
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
tab2 <- summary(mod2)
reg2 <- as.data.frame(tab2$solutions)
reg2
rreg2 <- round(reg2,4)
rreg2
rreg2 <- rreg2[,-c(4)]
rreg2

# mod1se <- mod1se[c(1, 3:11, 2, 12:20),]
# mod1se

# put BELIEF AND GROUP TOGETHER ODD Numbering
col1.2 <- rreg2[c(1,3:10), ]
col2.2 <- rreg2[c(2,11:18), ]
col1.2
col2.2


# 
# library(stargazer)
# head(now)
# starnow <- now[,c(2,6,7,8,9,10,13)]
# head(starnow)
# colnames(starnow) <- c("Age","Conservativism","Religious Ident","Rel Group Vals","NZDep")
# stargazer(starnow,summary=TRUE, omit.summary.stat=c("max","min","n"), type="text")



table2 <- cbind(col1.2,col2.2)
table2 
# rownames(table2) <- c("Intercept","Men","Belief(S)", "Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban", "Men X Belief(S)")

rownames(table2) <- c("Intercept","Men","Belief(S)", "Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban")
table2


table1
colnames(table2)[1] <- "Pray Means"
colnames(table2)[5] <- "Zero_Pray"
table2

star2 <- stargazer(table2,summary=FALSE, type="text")
star2

star2.a<- stargazer(table2,summary=FALSE, align=FALSE)











####################################################################################
## rel ingroup and prayer, what predicts it? 
####################################################################################

####################################################################
## Do men church more?   No Singles
####################################################################
head(data)

cor
%%
  
priorChurch = list( 
             R = list(V = diag(1)*.02, nu = 2),
             G = list( G1 = list(V = diag(1) *.02, nu =3),
                       G2 = list(V = diag(1) *.02, nu= 3)))

#We just have to be careful about how we express the results. Stating that the family variance is 0.749 is meaningless without putting it in the context of the assumed residual variance. It is therefore more appropriate to report the intraclass correlation which in this context is the expected correlation between the state Pupated/Not Pupated, for members of the same family. It can be calculated as: p.50
%
head(data)
mod3 <-  MCMCglmm(Church ~  trait *(Gender +
                                      RelGrpS + 
                                     Prtnr + 
                    AgeC +  
                    Eth + 
                    ConsrvS+ 
                    NZDepS +  
                    Urban), 
                  random = ~ trait:Denom + trait:ChildHome,
                  family = "zapoisson",
                  rcov = ~ trait:units, 
                  prior=priorChurch ,
                  nitt=nittA,
                  thin=thin,
                  burnin=burnin,
                #  pr = TRUE,
              #    pl = TRUE,
               #  saveX =TRUE,
               #  saveZ = TRUE,
                  data=data)




coefplot2(mod3$VCV, var.idx =c(1:3))
summary(mod3t$VCV)
summary(mod3)


## ICC  family/family + units

## Intraclass Correlation
mean((mod3$VCV[, 1])/(rowSums(mod3$VCV)) )  # Expected correlation for denominations = 40   0.31
mean((mod3$VCV[, 2]/(rowSums(mod3$VCV)) ))  # Expected correlation for children home = 40   0.04


names3<- c(("Intercept"),"ZA_Intercept","Men","RelGrp(S)", "Partner","Age(C)","NZ Euro", "Conservative(S)","NZDep(S)","Urban","ZA_Men","ZA_RelGrp(S)", "ZA_Partner","ZA_Age(C)","ZA_NZ Euro", "ZA_Conservative(S)","ZA_NZDep(S)","ZA_Urban")

coefplot2(mod3$Sol,intercept=TRUE,
          varnames=names3,
          col=c("gray","gray","red","red", rep("gray",6),"red","red", rep("gray",6)))


coefplot2(mod3$VCV, var.idx=c(1:16))
summary(mod3)
plot(mod3)


# #########  Extracting means + standard errors of means for the table
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
tab3 <- summary(mod3)
reg3 <- as.data.frame(tab3$solutions)
reg3
rreg3 <- round(reg3,4)
rreg3
rreg3 <- rreg3[,-c(4)]
rreg3


# put BELIEF AND GROUP TOGETHER ODD Numbering
col1.3 <- rreg3[c(1,3:10), ]
col2.3 <- rreg3[c(2,11:18), ]
col1.3
col2.3





table3 <- cbind(col1.3,col2.3)
table3 
rownames(table3) <- c("Intercept","Gender","RelGroup(S)", "Partner","Age","NZ Euro", "Conservative(S)","NZDep(S)","Urban")
table3



colnames(table3)[1] <- "Church Means"
colnames(table3)[5] <- "Zero_Church"
table3

star3 <- stargazer(table3,summary=FALSE, type="text")


star3.a<- stargazer(table3,summary=FALSE)







####################################################################################
## Children Generally within Religion
####################################################################################
head(ssgot)
dim(ssgot)
head(ssgot)
datafull <- ssgot[,c(2,3,5,9,11,13,14,16,17,19,28)]
nrow(datafull     )
head(datafull)
str(datafull)

colnames(datafull) <- c("Id","Urban","Partner","Children","ChildrenHome", "Ethnic","Conservative","Religious","Gender", "Age","NZDep")


library(mice)


completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

dfull <- completeFun(datafull,c( "Gender","Religious"))
nrow(dfull)
head(dfull)

hist(dfull$Children)

sum(nrow(subset(dfull, Children == 0)))/sum(nrow(dfull)) #.24

# leave out Children as these values will be imputed during MCMC
mifull <- mice(dfull[,-c(4)])
#head(mi)


#plot(mi)
f.now<- complete(mifull)

head(f.now)
f.cnow <- f.now

f.cnow$Children <- dfull$Children
head(f.cnow)

with(f.cnow, CrossTable(Children,Religious))
with(data, CrossTable(Child,Gender))

with(f.cnow, CrossTable(Religious, Children))
with(data, CrossTable(Eth))

with(data, CrossTable(Pray,Gender))

head(data)

tmean <- subset(f.cnow, Religious ==1)
mean(tmean$Children, na.rm=TRUE)
mean(data$Child, na.rm=TRUE)


f.cnowS <- f.cnow[,-c(1,4)]
library(stargazer)
stargazer(f.cnow, type ="text",omit.summary.stat=c("max","min","n"), digits=2)
stargazer(f.cnowS, omit.summary.stat=c("max","min","n"), digits=2)
stargazer(f.cnowS, omit.summary.stat=c("max","min","n"), digits=2, type="text")

stargazer(f.cnowS, omit.summary.stat=c("max","min","n"), digits=2, type="text")



head(cnow)
ddply(f.cnow, "Gender", summarise, Children.mean=mean(Children, na.rm =TRUE))

ddply(f.cnow, "Partner", summarise, children.mean=mean(Children, na.rm =TRUE))

ddply(f.cnow, "Ethnic", summarise, children.mean=mean(Children, na.rm =TRUE))

ddply(f.cnow, "Religious", summarise, Children.Mean=mean(Children, na.rm =TRUE), SD=sd(Children, na.rm =TRUE))

ddply(f.cnow, "Religious", summarise, Age.Mean=mean(Age), SD=sd(Age))



head(data)

staruse <- data[,c(2,3,4,5,6,7,9,10,11,12,14,15,16,17)]
stargazer(staruse,  omit.summary.stat=c("max","min","n"), digits=2, type="text")
stargazer(staruse,  omit.summary.stat=c("max","min","n"), digits=2)

# transform from integer
data$NZDep3 <- as.numeric(as.character(data$NZDep3))
str(cnow)


## Religon only
ddply(data, "Gender", summarise, P.Mean=mean(Pray), SD=sd(Pray))
ddply(data, "Gender", summarise, C.Mean=mean(Church), SD=sd(Church))
ddply(data, "Gender", summarise, C.Mean=mean(Child,na.rm=TRUE), SD=sd(Child,na.rm=TRUE))
ddply(data, "Gender", summarise, B.Mean=mean(Belief), SD=sd(Belief))
ddply(data, "Gender", summarise, ch.Mean=mean(ChildHome, na.rm=TRUE), SD=sd(ChildHome, na.rm=TRUE))

ddply(f.cnow, "Religious", summarise,  NZDep.mean=mean(NZDep), SD=sd(NZDep))
1647-568

ddply(cnow, summarise, NZDep.Mean=mean(NZDEP3), SD=sd(NZDEP3))


head(data)
stard <- data[,c(2,3,4,5,6,8,9,10,11,13,14,15,16)]

stargazer(stard, omit.summary.stat=c("max","min","n"), type ="text", digits=2) 
stargazer(stard, omit.summary.stat=c("max","min","n"),  digits=2) 

library(gmodels)

with(data, CrossTable(ChildHome))

with(data, CrossTable(Urban))
with(f.cnow, CrossTable(Urban))

with(f.cnow, CrossTable(Ethnic))
with(f.cnow, CrossTable(Gender))

with(f.cnow, CrossTable(Religious))

with(f.cnow, CrossTable(Partner))
with(data, CrossTable(Prtnr))

stargazer(f.cnowS, digits=2)


head(cnow)



colnames(cnow) <- c("id", "Age","Gender","Eth","Partn","Consrv","RelImp","Belief",
                      "RelGrp","Pray","Church","Denom","NZDep", "Urban","ChildHome","Child")
head(cnow)

ccnow <- cnow[,-c(1,3,4,5,7,12,14)]
head(ccnow)


############################################
##  USE
############################################
xtable(corstarsl(ccnow))
ts <- corstarsl(ccnow)

ts
stargazer(ts, summary=FALSE, type = "text")
##############################

library(corrgram)
corrgram(f.cnow, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.pie, text.panel=panel.txt,
         main="Correlation of Modelled Variables")

corrgram(ccnow, order=TRUE, lower.panel=panel.shade,
         upper.panel=panel.pie, text.panel=panel.txt,
         main="Correlation of Modelled Variables")



 f.cnow$Partner<- as.factor(f.cnow$Partner)
 f.cnow$Ethnic<- as.factor(f.cnow$Ethnic)  
# sum(nrow(subset(f.cnow, Children ==0)))/sum(nrow(f.cnow)) #.24
 str(f.cnow)


#add children back (we do not estimate, but let MCMC do this)



contrasts(f.cnow$Urban  )  <- contr.sum(2)

contrasts(f.cnow$Partner ) <- contr.sum(2)


contrasts(f.cnow$Ethnic  )  <- contr.sum(2)


contrasts(f.cnow$Religious)  <- contr.sum(2)

contrasts(f.cnow$Gender)  <- contr.sum(2)


f.cnow$AgeC    <- scale(f.cnow$Age  , center =TRUE, scale =FALSE)

f.cnow$ConservativeS     <- scale(f.cnow$Conservative , center =TRUE, scale =TRUE)
f.cnow$NZDepS    <- scale(f.cnow$NZDep  , center =TRUE, scale =TRUE)

head(f.cnow)
str(f.cnow)
# Descriptive Data
with(f.cnow, CrossTable(Children, Religious))
head(f.cnow)
f.cnowS <- f.cnow[,-c(1,4)]
library(stargazer)
stargazer(f.cnowS, type="text")

stargazer(f.cnowS, type ="text", digits=2)




priorFull = list( 
             R = list(V = diag(1), nu = 3))

modFpois <-  MCMCglmm(Children ~  1, 
                  family = "poisson", 
                  rcov = ~ trait:units, 
                  prior=priorFull, 
                  nitt=nitt,
                  thin=thin,
                  burnin=burnin,
                  pr = TRUE,
                 pl = TRUE,
                 saveX =TRUE,
                 saveZ = TRUE,
                 data=f.cnow)

str(f.cnow)


modFul <-  MCMCglmm(Children ~  trait * (Religious + 
                                      Gender +  
                                      AgeC +  
                                      Ethnic + 
                                      Partner +
                                      ConservativeS+ 
                                      NZDepS +  
                                      Urban), 
                  family = "zapoisson", 
                  rcov = ~ trait:units, 
                  prior=priorFull ,
                  nitt=50000,
                  thin=thin,
                  burnin=burnin,
                  #pr = TRUE,
                  #pl = TRUE,
                  #saveX =TRUE,
                  #saveZ = TRUE,
                  data=f.cnow)


nrow(f.cnow)



# Diagnostics
 mcmc.units <- colMeans(modFul$Liab) - predict(modFul,   marginal = NULL, type = "terms")
 par(mfrow=c(1,2))
 qqnorm(mcmc.units) ; qqline(mcmc.units)
 hist(mcmc.units, col = grey(.6), main = "Over-dispersion")

plot(modFpois)
#
summary(modFul)
summary(modFul)

namesF <- c(("Intercept"),"ZA_Intercept","Religious","Gender","Age(C)","NZ Euro", "Partner", "Conservative(S)","NZDep(S)","Urban","ZA_Religious","ZA_Gender","ZA_Age(C)","ZA_NZ Euro", "ZA_Partner", "ZA_Conservative(S)","ZA_NZDep(S)","ZA_Urban")

coefplot2(modFul$Sol, intercept=TRUE, 
          varnames = namesF,   col=c("gray","gray","red",rep("gray",7),"red",rep("gray",7))
)

coefplot2(modFul$VCV, intercept=TRUE)



## ICC  family/family + units

## Intraclass Correlation
head(data)
mean(data$Church)
sd(data$Church)
mean(exp(rowSums(modFul$VCV)))  # 1.068  dispersion

#

# Check contrasts in religion coding to make sure it is sume (-1,1)
str(f.cnow)
%%
# Magnitiude of difference in children comparing religious and non religious people
summary(modFul)
## not religious expected offspring   1.94

## Note that you vary the coding scheme to recover expected outcomes for different effects.

round(mean(exp(
    modFul$Sol[,1] * 1 + #intercept
    modFul$Sol[,3] * 1 +  #rel
    modFul$Sol[,4] * 0  + # gen
    modFul$Sol[,5] * 0 +  #age
    modFul$Sol[,6] * 0  + #eth
    modFul$Sol[,7] * 1 + #part
    modFul$Sol[,8] * 0 + #cons
    modFul$Sol[,9] * 0 +  #dep
    modFul$Sol[,10] * 0  #urban
)),2)

%

#Religous expected offspring  2.42
round(mean(exp(
    modFul$Sol[,1] * 1 + #intercept
    modFul$Sol[,3] * 0 +  #rel
    modFul$Sol[,4] * 0  + # gen
    modFul$Sol[,5] * 0 +  #age
    modFul$Sol[,6] * 0  + #eth
    modFul$Sol[,7] * 1 + #part
    modFul$Sol[,8] * 0 + #cons
    modFul$Sol[,9] * 0 +  #dep
    modFul$Sol[,10] *0  #urban
)),2)



## nonreligious and religious expected zeros (no diff)  =.48
round(mean(exp(-exp(
    modFul$Sol[,2] * 1 + #intercept
    modFul$Sol[,11] * 1 +  #rel
    modFul$Sol[,12] * 0  + # gen
    modFul$Sol[,13] * 0 +  #age
    modFul$Sol[,14] * 0  + #eth
    modFul$Sol[,15] * 1 + #part
    modFul$Sol[,16] * 0 + #cons
    modFul$Sol[,17] * 0 +  #dep
    modFul$Sol[,18] * 0  #urban
))),2)


 
# #########  Extracting means + standard errors of means for the table
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
tabF <- summary(modFul)
regF <- as.data.frame(tabF$solutions)
regF
rregF <- round(regF,4)
rregF
rregF <- rregF[,-c(4)]
rregF

# mod1se <- mod1se[c(1, 3:11, 3, 13:20),]
# mod1se

# put BELIEF AND GROUP TOGETHER ODD Numbering
col1.F <- rregF[c(1,3:10), ]
col2.F <- rregF[c(2,11:18), ]
col1.F
col2.F


tableF <- cbind(col1.F,col2.F)
tableF
rownames(tableF) <- c("Intercept","Religious","Gender","Age(C)","NZ Euro", "Partner", "Conservative(S)","NZDep(S)","Urban")
tableF


tableF
colnames(tableF)[1] <- "Children"
colnames(tableF)[5] <- "Zero_Children"
tableF

starF <- stargazer(tableF,summary=FALSE)
starF <- stargazer(tableF,summary=FALSE, type="text")


starF.a<- stargazer(tableF,summary=FALSE)



summary(mod5)

#######################################
## Children within Religion
#######################################


data$Gender)

prior1 = list( 
             R = list(V = diag(1)*.02, nu = 3),
             G = list( G1 = list(V = diag(1) *.02, nu = 3)))

prior2 = list( 
             R = list(V = diag(1)*.02, nu = 3),
             G = list( G1 = list(V = diag(2) *.02, nu = 3)))




mod5pois <-  MCMCglmm(Child ~  1, 
                  random = ~ trait:Denom,
                  family = "poisson", 
                  rcov = ~ trait:units, 
                  prior=prior, 
                  nitt=nitt,
                  thin=thin,
                  burnin=burnin,
                  pr = TRUE,
                  pl = TRUE,
                 # saveX =TRUE,
                  #saveZ = TRUE,
                  data=data)


str(data)


  
head(amdata)
mitest <- mice(amdata[,-c(16)])
ddat<- complete(mitest)
head(ddat)  

colnames(ddat) <- c("Id","Age","Gender","Eth","Prtnr", "Consrv", "Child", "Relimp","Belief","RelGrp","Pray","Church","Denom","NZDep","Urban")

library(arm)
Child,
%
dota <- transform(subset(ddat,select=c(Id,Age,Gender,Eth,Prtnr,Consrv,Child,Belief,RelGrp,Pray,Church,Denom,NZDep,Urban)),
                 Id=factor(Id),
                 Age=rescale(Age,"center"),
                 Gender=rescale(Gender, binary.inputs= "center"),
                 Eth=rescale(Eth, "center"),
                 Prtnr=rescale(Prtnr,"center"),
                 Belief=rescale(Belief, "full"),
                 RelGrp=rescale(RelGrp, "full"),
                 Pray=rescale(Pray,"full"),
                 Church=rescale(Church,"full"),
                 Denom=factor(Denom),
                 Consrv = rescale(Consrv, "full"),
                 NZDep=rescale(NZDep,"full"))
                 
                 
str(doata)
head(dota)
nrow(dota)
sum(is.na(dota))

prior1 = list( 
             R = list(V = diag(1)*.02, nu = 3),
             G = list( G1 = list(V = diag(1) *.02, nu = 3)))

head(dota)

mod5tt <-  MCMCglmm(Child ~            Pray + 
                                      Church +  
                                      Gender + 
                                      Age +  
                                      Eth + 
                                      Prtnr +
                                      Consrv + 
                                      NZDep +  
                                      Urban, 
                  random = ~ Denom,
                  family = "poisson", 
                  rcov = ~ trait:units, 
                  prior=prior1,
                 data=dota)


  
  
mod5 <-  MCMCglmm(Child ~  trait * (  PrayS + 
                                      ChurchS +  
                                      Gender + 
                                      AgeC +  
                                      Eth + 
                                      Prtnr +
                                      ConsrvS + 
                                      NZDepS+  
                                      Urban), 
                  random = ~ trait:Denom,
                  family = "zapoisson", 
                  rcov = ~ trait:units, 
                  prior=prior,
                  #nitt= 200000,  # Mixing is an issue
                  #thin=thin,
                  #burnin=burnin,
                  slice =FALSE,
                  #pr = TRUE,
                 # pl = TRUE,
                  #saveX =TRUE,
                  #saveZ = TRUE,
                 data=data)

rm(mod5t)

head(telgot)
%

                 
%mod$Residual.nrt <- 2
mod5ttt$Residual.nrt  <- 2
summary(mod5ttt)
coefplot2(mod5ttt$VCV)
summary(mod5)
coefplot2(list(mod5t,mod5tt), intercept=TRUE,merge.names=FALSE)



mod5 <-  MCMCglmm(Child ~  trait * (  PrayS + 
                                      ChurchS +  
                                      Gender + 
                                      AgeC +  
                                      Eth + 
                                      Prtnr +
                                      ConsrvS + 
                                      NZDepS+  
                                      Urban), 
                  random = ~ trait:Denom,
                  family = "zapoisson", 
                  rcov = ~ trait:units, 
                  prior=prior1,
                  nitt= 500000,  # Mixing is an issue
                  #thin=thin,
                  #burnin=burnin,
                  slice =FALSE,
                  #pr = TRUE,
                 # pl = TRUE,
                  #saveX =TRUE,
                  #saveZ = TRUE,
                 data=data)

str(data)


summary(mod5)

str(data)
mean(data$Church) 
sd(data$Pray)
## not religious expected offspring   1.94
round(mean(exp(
    mod5$Sol[,1] * 1 + #intercept
    mod5$Sol[,3] * 1 +  #pray
    mod5$Sol[,4] * 0  + # church
    mod5$Sol[,5] * 0 +  #men
    mod5$Sol[,6] * 0  + #age
    mod5$Sol[,7] * 0 + #eth
    mod5$Sol[,8] * 1 + #part
    mod5$Sol[,9] * 0 +  #cons
    mod5$Sol[,10] * 0 + #nzdep
    mod5$Sol[,11] * 1  #urban
)),2)
%
round(mean(exp(-exp(
    mod5$Sol[,2] * 1 + #intercept
    mod5$Sol[,12] *0 +  #pray
    mod5$Sol[,13] * 1  + # church
    mod5$Sol[,14] * 0 +  #men
    mod5$Sol[,15] * 0  + #age
    mod5$Sol[,16] * 0 + #eth
    mod5$Sol[,17] * 1 + #part
    mod5$Sol[,18] * 0 +  #cons
    mod5$Sol[,19] * 0 + #nzdep
    mod5$Sol[,20] * 1  #urban
))),2)

sd(data$Church)






## Prior for Poiss, Zap, and Gauss modesl

head(data)
str(data)
ddata <- data
ddata$Pray <- as.numeric(as.character(ddata$Pray))
ddata$Church <- as.numeric(as.character(ddata$Church))


coefplot2(mod5$Sol, intercept =TRUE)
summary(mod5)

data)

names5 <- c(("Intercept"),"ZA_Intercept", "Prayer(C)","Church(C)","Gender","Age(C)","NZEuro","Partner", "Conservative (C)","NZDep(C)","Urban", "ZA_Prayer(C)","ZA_Church(C)","ZA_Gender","ZA_Age(C)","ZA_NZEuro","ZA_Partner", "ZA_Conservative (C)","ZA_NZDep(C)","ZA_Urban")

coefplot2(mod5$Sol, intercept=TRUE, 
          varnames = names5, 
          col = c("gray","gray","red","red","green","gray","gray","gray","gray","gray","gray",'red','red',"green","gray","gray","gray","gray","gray","gray"))
          
          


## not religious expected offspring   1.94
round(mean(exp(
    mod5$Sol[,1] * 1 + #intercept
    mod5$Sol[,3] * 0 +  #pray
    mod5$Sol[,4] * 1  + # church
    mod5$Sol[,5] * 0 +  #men
    mod5$Sol[,6] * 0  + #age
    mod5$Sol[,7] * 0 + #eth
    mod5$Sol[,8] * 1 + #part
    mod5$Sol[,9] * 0 +  #cons
    mod5$Sol[,10] * 0 + #nzdep
    mod5$Sol[,11] * 1  #urban
)),2)
%
round(mean(exp(-exp(
    mod5$Sol[,2] * 1 + #intercept
    mod5$Sol[,12] * 0 +  #pray
    mod5$Sol[,13] * 0  + # church
    mod5$Sol[,14] * 0 +  #men
    mod5$Sol[,15] * 0  + #age
    mod5$Sol[,16] * 0 + #eth
    mod5$Sol[,17] * 1 + #part
    mod5$Sol[,18] * 0 +  #cons
    mod5$Sol[,19] * 0 + #nzdep
    mod5$Sol[,20] * 1  #urban
))),2)

          mean(data$Church)  
  sd(data$Church)
mean(data$Pray) 
sd(data$Pray)


library(coefplot2)

summary(mod5)
coefplot2(mod5, var.idx=c(1:20))
coefplot2(mod5$VCV, var.idx=c(1:2))


coefplot2(mod5$Sol, col=(c("red","blue")),var.idx=c(3,4,12,13))

# Denomination effects
ic.5 <- (mod5$VCV[, 2])/(rowSums(mod5$VCV)) # Expected correlation for denominations = .13
round(mean(ic.5),2) # 56
plot(ic.5)
#plot(mod5)


zBack  <- fun{x,mu,sd}(x-mu/sd)
zBack(1.2,8,7)


# #########  Extracting means + standard errors of means for the table
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
tab <- summary(mod5)
reg <- as.data.frame(tab$solutions)
reg
rreg <- round(reg,4)
rreg
rreg <- rreg[,-c(4)]
rreg

# mod1se <- mod1se[c(1, 3:11, 2, 12:20),]
# mod1se

mod5tab <- rreg[c(1, 3:11, 2, 12:20),]
col1.5 <- mod5tab[c(1:10),]
col2.5 <- mod5tab[c(11:20),]




table5 <- cbind(col1.5,col2.5)
table5
rownames(table5) <- c("Intercept","Prayer(S)","Church(S)","Gender","Age(S)","NZ Euro","Partner", "Conservative (S)","NZDep(S)","Urban")
table5


colnames(table5)[1] <- "Children"
colnames(table5)[5] <- "Zero_Children"
table5

stargazer(table5, summary =FALSE, type="text")

stargazer(table5, summary =FALSE)

star1 <- stargazer(table5,summary=FALSE, type="text")




# Diagnostics
 mcmc.units <- colMeans(mod5C$Liab) - predict(mod5C, 	marginal = NULL, type = "terms")
 par(mfrow=c(1,2))
 qqnorm(mcmc.units) ; qqline(mcmc.units)
 hist(mcmc.units, col = grey(.6), main = "Over-dispersion")
%


# Check plots
plot(mod5C)
#
#mod5$Residual$nrt <- 2  # this is from a bug in mcmcglmm at the time. (2012)
summary(mod5C$Sol)
summary(mod5$VCV)
summary(mod5)
#




sd(cornow$Pray)
rm(mm)
summary(mod5)
rm(lm)
lm <-  with(cornow, expand.grid(
        Child = 0, 
        PrayC = seq(0, 2, by =.1),
        ChurchC = 0, 
        #PrayC= 0,
        #ChurchC = seq(-2,0, by=.25),# Someone who goes to church every week
        Gender = 0,
        AgeCy = 0, 
        Eth = 0,
        Prtnr =1,
        ConsrvC = 0,
        NZDepC = 0,
        Urban  =1   
 ))

#check
#lm

mm<- model.matrix(Child ~ PrayC + ChurchC + Gender + AgeCy + Eth + Prtnr + ConsrvC+ NZDepC + Urban, data=lm)
mm



#residual variance

# od.mcmc5 <- mean(mod5C$VCV[,1])
# head(od.mcmc5)


wm <- as.data.frame(mm)
mm

dim(wm)


#########  Extracting means + standard errors of means for the plot
library(MCMCpack)
library(Hmisc)


## Here we yank the posterior means (above we used the modes.)
rm(mod5se)
qreg <- summary(mod5$Sol)
str(qreg)
mod5se <-as.data.frame(qreg$statistics)


# put ZA's with each other
mod5se <- mod5se[c(1, 3:11, 2, 12:20),]

# Yank Mean + Time-series SE
ymod5se  <- mod5se[, c(1,4)]
ymod5se
rm(meanprP)
meanprP

#Get col means
meanprP  <-  matrix(mod5se[1:10,1])

#Get TIme-serioes
diffP
mod5se
diffP <- matrix(mod5se[1:10,4])
lop <-  meanprP - diffP
hip <-  meanprP + diffP

head(meanprPZ)
meanprPZ  <-  matrix(mod5se[11:20,1])
diffPZ <- matrix(mod5se[11:20,4])

## create interval, not we will be exponentiating these numbers, so lo and hi don't have the same meaning.

lopZ <-  meanprPZ - diffP
hipZ <-  meanprPZ + diffP

# print to check
lopZ
hipZ

## Note you can get the same values using the coefplot2 package.  Just create a df, and yank, but remember that for matrix math you need to put the data into matrices.
##coeftab(mod5$Sol,quantiles = TRUE, cmult=1.96) # quantiles =FALSE,clevel=c(.5,0,.95),cmult=1.96)


rm(W)
W  <- as.matrix(mm)
W


## Recover data scale for the Poisson process. 
exp(1)


childrenPm <- exp(((W %*% meanprP)) )
plo <-  (exp(W %*% lop ) )
phi <-  (exp(W %*% hip  ))

## Recover data scale for zero-inflation
zchildrenPm <-  (exp(-exp((W %*% meanprPZ) )))   
zplo <-  (exp(-exp((W %*% lopZ) )))
zphi <- (exp(-exp((W %*% hipZ) ))) 

## Data frame
wm <- as.data.frame(mm)

## Check
head(wm)

## Chuck predited values into the data frame
wm$children <- as.numeric(childrenPm)

# Check
head(wm)

# Works, next...
wm$upper  <- as.numeric(plo)
wm$lower <-as.numeric(phi)
wm$zerochildren <- as.numeric(zchildrenPm)
wm$zeroupper  <- as.numeric(zplo)
wm$zerolower <-as.numeric(zphi)


head(wm)
colnames(wm) <- c("Intcpt","Pray","Church","Gender","Age","Eth","Part","Cons","NZDep",'Urban',"Children","Upper","Lower", 
                  "ZeroChildren","ZeroUpper","ZeroLower")
dim(wm)
str(wm)
head(wm)



## Ok onto the plots

library(ggplot2)
%%
p <- ggplot(wm, aes(x= Pray,y= Children)) + geom_point()
p1 <- p + geom_errorbar(aes(ymin = Lower, ymax = Upper) ,alpha=.3, colour="blue") + scale_y_continuous(limits=c(2.2,2.65)) + theme_bw() + labs(title = "Weekly prayer and expected fertility")
p1


zp <- ggplot(wm, aes(x= Pray,y= ZeroChildren)) + geom_point()
zp1 <- zp + geom_errorbar(aes(ymin = ZeroLower, ymax = ZeroUpper) ,alpha=.3, colour=" blue")  + scale_y_continuous(limits=c(.3,.45))+ theme_bw() + 
  labs(title = "Weekly prayer and expected prob of zero children")

zp1



library(gridExtra)

mygrid <- grid.arrange(p1, zp1, nrow=2) 



summary(mod5)
rm(clm)
clm <-  with(cornow, expand.grid(
        Child = 0, 
        #PrayC = seq(0, 50, by =1),
       # ChurchC = 0, 
        PrayC= 0,
        ChurchC = seq(0,2, by=.1),# Someone who goes to church every week
        Gender = 0,
        AgeCy = 0, 
        Eth = 0,
        Prtnr =1,
        ConsrvC = 0,
        NZDepC = 0,
        Urban  =1   
 ))

#check
#lm

cmm<- model.matrix(Child ~ PrayC + ChurchC + Gender + AgeCy + Eth + Prtnr + ConsrvC+ NZDepC + Urban, data=clm)
cmm


#residual variance

#od.mcmc5 <- mean(mod5C$VCV[,1])
#head(od.mcmc5)


cwm <- as.data.frame(cmm)
cmm

dim(cwm)



qreg <- summary(mod5C$Sol)
str(qreg)
mod5se <-as.data.frame(qreg$statistics)

# put ZA's with each other
mod5se <- mod5se[c(1, 3:11, 2, 12:20),]

# Yank Mean + Time-series SE
ymod5se  <- mod5se[, c(1,4)]
ymod5se
meanprP

#Get col means
meanprP  <-  matrix(mod5se[1:10,1])

#Get TIme-serioes
diffP
diffP <- matrix(mod5se[1:10,4])
lop <-  meanprP - diffP
hip <-  meanprP + diffP


meanprPZ  <-  matrix(mod5se[11:20,1])
diffPZ <- matrix(mod5se[11:20,4])

## create interval, not we will be exponentiating these numbers, so lo and hi don't have the same meaning.

lopZ <-  meanprPZ - diffP
hipZ <-  meanprPZ + diffP

# print to check
lopZ
hipZ

## Note you can get the same values using the coefplot2 package.  Just create a df, and yank, but remember that for matrix math you need to put the data into matrices.

##coeftab(mod5$Sol,quantiles = TRUE, cmult=1.96) # quantiles =FALSE,clevel=c(.5,0,.95),cmult=1.96)


rm(cW)
cW  <- as.matrix(cmm)
cW

cchildrenPm <- exp(((cW %*% meanprP)) )
cplo <-  (exp(cW %*% lop ) )
cphi <-  (exp(cW %*% hip) )


## Recover data scale for zero-inflation
czchildrenPm <-  (exp(-exp((cW %*% meanprPZ) )))   
czplo <-  (exp(-exp((cW %*% lopZ) )))
czphi <- (exp(-exp((cW %*% hipZ) ))) 

## Data frame
cwm <- as.data.frame(cmm)

## Check
head(cwm)

## Chuck predited values into the data frame
cwm$children <- as.numeric(cchildrenPm)

# Check
head(cwm)

# Works, next...
cwm$upper  <- as.numeric(cplo)
cwm$lower <-as.numeric(cphi)
cwm$zerochildren <- as.numeric(czchildrenPm)
cwm$zeroupper  <- as.numeric(czplo)
cwm$zerolower <-as.numeric(czphi)


head(cwm)
colnames(cwm) <- c("Intcpt","Pray","Church","Gender","Age","Eth","Part","Cons","NZDep",'Urban',"Children","Upper","Lower", 
                  "ZeroChildren","ZeroUpper","ZeroLower")
dim(cwm)
str(cwm)
head(cwm)



## Ok onto the plots

library(ggplot2)



library(ggplot2)
p <- ggplot(wm, aes(x= Pray,y= Children)) + geom_point()
p1 <- p + geom_errorbar(aes(ymin = Lower, ymax = Upper) ,alpha=.3, colour="blue") + scale_y_continuous(limits=c(2.2,2.6)) + theme_bw() + labs(title = "Weekly prayer rate and expected fertility")
p1


zp <- ggplot(wm, aes(x= Pray,y= ZeroChildren)) + geom_point()
zp1 <- zp + geom_errorbar(aes(ymin = ZeroLower, ymax = ZeroUpper) ,alpha=.3, colour=" blue")  + scale_y_continuous(limits=c(.3,.55))+ theme_bw() + 
  labs(title = "Weekly prayer rate and probability of no children")

zp1

%%
cp <- ggplot(cwm, aes(x= Church,y= Children)) + geom_point()
cp
cp1 <- cp + geom_errorbar(aes(ymin = Lower, ymax = Upper) ,alpha=.3, colour="red") + theme_bw() + scale_y_continuous(limits=c(2.2,2.6)) +labs(title = "Monthly church rate and expected fertility")
cp1


czp <- ggplot(cwm, aes(x= Church,y= ZeroChildren)) + geom_point()
czp1 <- czp + geom_errorbar(aes(ymin = ZeroLower, ymax = ZeroUpper) ,alpha=.3, colour="red") + theme_bw()  +scale_y_continuous(limits=c(.3,.55))+ labs(title = "Monthly church rate and probability of no children")
czp1




library(gridExtra)
mygrid <- grid.arrange(p1,cp1,zp1,czp1, clip=TRUE, nrow=2) 




####################################################################################################################
## BELOW CAN BE IGNORED  THIS IS ALL OLD ANLAYSIS. I leave it so  that you can follow the process of our thinking.
######################################################################################################################



#############################################
# compare religious and not religious people
#############################################
head(ssgot)

subg  <- subset(ssgot, is.na(gen3) ==FALSE & is.na(gen35)==FALSE & is.na(rels3) ==FALSE)

nrow(subg)
head(subg)

fdata <- subg[ , c("id",
                      "age3", 
                      "gen3",
                      "ethn3", 
                      "partner3",
                      "pol3", 
                       "rels3",
                      "children3",                       
                      "NZDEP3"
)]

head(ssgot)

nrow(fdata)
head(fdata)



#####################################
# AMELIA MODELS
#####################################




## make prediction data frame
#newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female"))
## design matrix (fixed effects)
mm <- model.matrix(delete.response(terms(fm2)),newdat)
## linear predictor (for GLMMs, back-transform this with the
##  inverse link function (e.g. plogis() for binomial, beta;
##  exp() for Poisson, negative binomial
newdat$children <- mm %*% fixef(fm2)

predvar <- diag(mm %*% vcov(fm2) %*% t(mm))
newdat$SE <- sqrt(predvar) 
newdat$SE2 <- sqrt(predvar+fm2$alpha^2)
#(Probably overly complicated) ggplot code:
head(dat)

library(ggplot2) 
pd <- position_dodge(width=0.4)
theme_update(theme_bw())
g0 <- ggplot(dat,aes(x=pray,y=children3,colour=gen35))+
  stat_sum(alpha=0.2,aes(size=..n..))+
  scale_size_continuous(breaks=1:4,range=c(2,5))
g0
g1 <- g0+geom_line(data=newdat,position=pd)+
  geom_point(data=newdat,shape=17,size=3,position=pd)
g1 + geom_linerange(data=newdat,
                    aes(ymin=children-2*SE,ymax=children+2*SE), position=pd)



## prediction intervals 
g1 + geom_linerange(data=newdat,
                    aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)


data <- na.omit(relgot)

newdat <-
  with(data, expand.grid(
    pray = seq(min(pray), max(pray), length = 10),
    Gen35 = levels(factor(gen35)),
    church= seq(min(church),max(church), length=5),
    children3= seq(min(children3), max(children3), length = 10)
  ))

#NZDepT35
fm2 <- glmmadmb(children3 ~ gen35 + prayS + churchS, random = ~ gen35|denom,  data = dat, family="nbinom1",zeroInflation = TRUE)
fm2 <- lmer(children3 ~ gen35 + prayS + churchS + 1|gen35,  data = dat, family="poisson")
summary(fm2)

rm(newdat)
newdat <-
  with(data, expand.grid(
    prayS = seq(0,2, by = .1),
    #seq(min(pray), max(pray), length = 10),
    gen35 = levels(factor(gen35)),
    churchS= seq(0,2, by =.1),
    children3=0,
    age35c=0
    #denom= levels(factor(denom))    
    #children3= seq(min(children3), max(children3), length = 10)
  ))


head(newdat)

newdat$childrenborn <- stats::predict(fm2, newdata=newdat, type = "response")



library(ggplot2)
qplot(prayS,childrenborn, data=newdat, colour=gen35) + geom_point()+theme_bw()
qplot(churchS,childrenborn, data=newdat, colour=gen35) + geom_point()+theme_bw()
#qplot(time, arousal, data=grid, colour=factor(role)) + geom_line(data=grid)

#Aerr <- stats::predict(fm2, newdata=newdat,  se = TRUE)
#Aerr <- stats::predict(fm2, newdata=newdat,  se = TRUE)
#head(Aerr)

library(MCMCglmm)

prior.1 =
  list(
    R = list(
      V = diag(2), nu = 0.002, fix =2),
    G = list(
      G1 = list(V = diag(4)*.002, nu = 4)))
%
# Subcripts too long to print so we relable them

#head(cornow)
#head(cornow)

# Not mixing, adust constrasts


k<-rpois(1000,3)
hist(k)



test7a<-  MCMCglmm(Child ~   PrayS  + ChurchS + Gender + Age + Eth + Prtnr + Consrv + NZDepS), random= ~idh(trait:Gender):Denom,
                  family = "zapoisson", rcov = ~idh(trait):units, prior=prior.1, data = cornow, pl=FALSE,pr=FALSE, saveX=FALSE)

test7<-  MCMCglmm(Church ~  trait * ( RelImpS  + RelGrpS + Gender + Age + Eth + Prtnr + Consrv + NZDepS), random= ~idh(trait:Gender):Denom,
                  family = "zapoisson", rcov = ~idh(trait):units, prior=prior.1, data = cornow, pl=FALSE,pr=FALSE, saveX=FALSE)


test7$Residual$nrt <- 2  
coefplot2(test7$Sol, intercept =TRUE)
#
#summary(test7$Sol)
summary(test7)
summary(test7$Sol)[1:3]
coefplot2(test7,intercept =TRUE)
#
plot(test7)
#
# qreg <- summary(test7)
# str(qreg)
# 
# 
# #List of 6 
# # $ statistics: num [1:3, 1:4] 0.3479 0.9029 0.3413 0.2602 0.0391 ... 
# #  ..- attr(*, "dimnames")=List of 2 
# #  .. ..$ : chr [1:3] "(Intercept)" "x" "sigma" 
# #  .. ..$ : chr [1:4] "Mean" "SD" "Naive SE" "Time-series SE" 
# # $ quantiles : num [1:3, 1:5] -0.139 0.825 0.257 0.168 0.878 ... 
# #  ..- attr(*, "dimnames")=List of 2 
# #  .. ..$ : chr [1:3] "(Intercept)" "x" "sigma" 
# #  .. ..$ : chr [1:5] "2.5%" "25%" "50%" "75%" ... 
# ... 
# # - attr(*, "class")= chr "summary.mcmc" 
# 
# toLatex(qreg$solutions,file="")
# 
# Table(qreg$quantiles,file="") 


library(xtable)
library(Hmisc)
library(texreg)
library(coefplot2)




coefplot2(test2, intercept=TRUE,)

library(coefplot)
ploto <- coefplot(test2, lwdInner=1, lwdOuter=.5, pointsize = .5, fillColor= "empty", secret.weapon=F)
multiplot(glm.nb, glm.out, lwdInner=1, lwdOuter=.5, pointSize = 2, )
ploto
rm(plot)
plotcoef(test2)

extract.coef(test2)

tab <- round(coeftab(test2, ptype="fixef", digits=3, quantiles=FALSE, ),2)
tab

library(gridExtra)
grid <- tableGrob(tab,  nrow=1)

gridd <- grid.arrange(plot, grid,  nrow=1)
griddd <- grid + labs(title = "Children")

gridd


xtable(tab)

tsum <- summary(test2)
coeftab(mod5$Sol,quantiles = TRUE, cmult=1.96) # quantiles =FALSE,clevel=c(.5,0,.95),cmult=1.96)

latex(obsum)
screenreg(list(obsum))
exp(.69 )


# look at people under 42

mean(relgot$age35, na.rm=T)

test3 <-  glmmadmb(children3 ~  gen35 + age35c + ethn3 + pol3c + factor(edu3) + lincome3c + prayS + churchS, random = ~gen35|denom,  
                   family="nbinom1",zeroInflation = TRUE, data = now)


## Check young people
now2 <- subset(now, age35c <= -8)
nrow(now2) # 453


test3 <-  glmmadmb(children3 ~  gen35 + age35c + ethn3 + pol3c + factor(edu3) + lincome3c + prayS + churchS, random = ~gen35|denom,  
                   family="nbinom1",zeroInflation = TRUE, data = now2)

 gap <- summary(test2)
tableGrob(gap, nrow=1, as.table=TRUE, main="test main",
sub=textGrob("test sub", gp=gpar(font=2)))


plogis(-4)

summary(test2)

#warnings()
#remove(combined.results)
#combined.results <- mi.meld(q = b.out, se = se.out)
#combined.results




#####################################
# MODESL REL NOT REL
#####################################

library(MCMCglmm)
library(coefplot2)
#intercpet only


completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}




#~cor(trait):units

### p.108 course notesWe can then test for zero-flation by constraining the over-dispersion to be the same for both process using a trait by units interaction in the R-structure, and by setting up the contrasts so that the zero-altered regression coecients are expressed as di↵erences from the Poisson regression coecients. When this di↵erence is zero the variable causes no zero-flation, when it is negative it causes zero-inflation and when it is positive it causes zero-deflation:

#
#
#
#mean(exp(modF.1$Sol[,1] + modF.1$Sol[,3]*1 + modF.1$Sol[,4] +  modF.1$Sol[,5]*0  +  modF.1$Sol[,6]*0 + modF.1$Sol[,7]*0 ) )


#(exp(modF.1$Sol[,1] + modF.1$Sol[,3] + modF.1$Sol[,4]*1 +  modF.1$Sol[,5]*0  +  modF.1$Sol[,6]*0 + modF.1$Sol[,7]*0 ) )


#HURDLE
#summary(modF.1 <- (MCMCglmm(children3 ~ trait - 1 +  at.level(trait, 1):gen35  + at.level(trait, 1):age35c +  at.level(trait, 1):ethn3 +  at.level(trait, 1):rels35 + at.level(trait, 1):pol3c, rcov = ~idh(trait):units, data = cagot, family = "hupoisson",  nitt=13000, prior = prior.01) ) )

#mean(exp(modF.1$Sol[,1] + modF.1$Sol[,3] + modF.1$Sol[,4]*1 +  modF.1$Sol[,5]*0  +  modF.1$Sol[,6]*0 + modF.1$Sol[,7]*0 ) )

# zero inflation in hurdle model
#c2 <- (16 * sqrt(3)/(15 * pi))^2

#HPDinterval(plogis(modF.1$Sol[, 2]/sqrt(1 + c2)))


### THEORETICAL MODEL ###   USE 

# prior for poisson
prior.ch  =
    list(R = list(V = diag(2), nu = 0.02, fix = 2))

table(sagot$children3)

summary(modF.1I <- (MCMCglmm(children3 ~  trait, rcov = ~idh(trait):units,  data = sagot, family = "zapoisson",   prior = prior.ch) ) )



### USE

summary(modF.1 <- (MCMCglmm(children3 ~  trait *( gen35  + age35c +  ethn3 +   pol3c+ rels35), rcov = ~idh(trait):units,  data = sagot, family = "zapoisson",   prior = prior.ch) ) )


exp(mean((modF.1$Sol[,1] + modF.1$Sol[,3]*1 + modF.1$Sol[,4] +  modF.1$Sol[,5]*0  +  modF.1$Sol[,6]*1 + modF.1$Sol[,7]*1 +  modF.1$Sol[,8]*1 ) ))

plot(modF.1)
#
#)


head(relgot)

#######################################################
# religious importance   NO -- EXPLORATION
##############

exp(.13)



######################################################################
##  Does religion predict CHILDREN  FULL MODEL USE  (no interaction)  EXPLORE
########################################################

# Cautious here

#########################################################################
## Do religious women have more children ?  ignore
#########################################################################



#########################################################women#####
#  Do religious men have more children?   ingnore 
#######

#Says that religious importance predicts fertility in men but not women





###############################################################################
###############################################################################
##############################################################################


################################################################################
######  Analysis of religion  Do religious men and women express different signals
###############################################################################




# GENDER OF SUBSAMPLE
table(relgot$gen35)
exp(.63 + .13)

#head(relgot)
dim(relgot)
rre#GENDER
head(relgot)

a <- relgot$gen3

b <- subset(a, a==1)#male
c <- subset(a, a==0)#female
d <- subset(a, is.na(a))
#length(d)
length(b) # male is 1568
length(c) # female is 2813
length(b)/(length(b)+length(c)) #ratio of men/women [1] 33
length(c)/(length(b)+length(c)) #ratio of women/men [1] 67



## obtain cases
#corgot <- completeFun(relgot, c("parentsRel", "rels"))



#summary(glm(rels ~ parentsRel, data = corgot, family = "binomial"))

#plogis(1.84) # Simple effect (.86)
#plogis(-1.933) # probablity of being religious at all .12


######################
# Tests
#####################

mean(relgot$belief, na.rm =TRUE)

gen0 <- subset(relgot, gen35 ==0)
gen1 <- subset(relgot, gen35 ==1)



# yes
t.test(as.numeric(gen0$belief),as.numeric(gen1$belief), na.rm=TRUE)

# Not significant pray
t.test(as.numeric(gen0$pray),as.numeric(gen1$pray), na.rm=TRUE)

# yes
t.test(as.numeric(gen0$church),as.numeric(gen1$church), na.rm=TRUE)

# Not significant script -- no
t.test(as.numeric(gen0$script),as.numeric(gen1$script), na.rm=TRUE)

# Importance  No 
t.test(as.numeric(gen0$rimp35),as.numeric(gen1$rimp35), na.rm=TRUE)

# #relnar
t.test(as.numeric(relmen$relnar),as.numeric(relfem$relnar), na.rm=TRUE)

data:  as.numeric(relmen$relnar) and as.numeric(relfem$relnar)
t = 1.8753, df = 979.32, p-value = 0.06104
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.006792002  0.299442177
sample estimates:
mean of x mean of y 
 3.714579  3.568254



boxplot(gen0$belief, gen1$belief)





###################################
# Function for creating dataframe
# ALL cases for predictor variables      
###################################

#NOTE WE DO NOT INCLUDE INCOME BECAUSE IT COMPROMISES DATA



completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}





#completeFun(DF, "y")


head(relgot)


rrelwork  <- completeFun(relgot, c("id", "partner3", "employ",  "gen35", "children3", "age35c", "ethn3", "pol3",  "beliefc",  "relnarc", "rimp35c", "denom"))





nrow(rrelwork)
 

mens <- subset(rrelwork, gen35==1)
fems <- subset(rrelwork, gen35==0)


nrow(rrelwork) # 1138


library(MCMCglmm)
library(coefplot2)
library(MuMIn)


prior.11 =
    list(R = list(V = 1, nu = 0.002),
         G = list(
         G1 = list( V= 1, nu = .002)))

#zero poiss
prior.1 =
    list(R = list(V = diag(1), nu = 0.002, fix = 1),
         G = list(
         G1 = list( V= diag(2)*.002, nu = 3)))





#head(rrelwork)
 #     summary(relwork$church)


## Sample without NAs in the predictors

margot <- completeFun(relgot, c("id", "beliefc", "age35c","gen35", "rels35", "partner3", "gen35", "ethn3", "pol3c", "rimp3", "relnarc","prayS", "churchS", "beliefc", "children3c", "denom","Religion.L2.Cats.T35","Religion.L3.Cats.T35","Religion.L4.Cats.T35"))






head(relgot)

nrow(margot)

head(margot)

#full    PRAYER IS ABOUT CHILDREN

#summary(modR <- (MCMCglmm(children3 ~  gen35 + age35c +  ethn3  + beliefc  + rimp3 + pol3c + relnarc + prayS + churchS,  family = "gaussian",  nitt=13000, prior = prior.11, data = margot ) ) )


summary(modN.I <- (MCMCglmm(belief ~ 1, random =  ~idh(gen35):denom, rcov = ~units,  family = "ordinal",  data = margot, prior = prior.1)))

summary(modN.I)



summary(modN <- (MCMCglmm(belief ~  gen35 + age35c +  ethn3  + partner3   + pol3c + children3c + prayS + churchS, ~idh(gen35):denom, rcov = ~units,  family = "ordinal", prior = prior.1, data = margot)))


summary(modN)

coefplot2(modN, intercept=TRUE)


#summary(modN <- (MCMCglmm(rimp35 ~  gen35 + age35c +  ethn3  + partner3   + pol3c + children3c + prayS + churchS, ~idh(gen35):denom, rcov = ~units,  family = "gaussian", prior = prior.1, data = margot )))


summary(modN)





str(margot)

## Church is about ingroup commitment


round((eta<-mean(modN$Sol[,1] * 1 + modN$Sol[,2] * 1  + modN$Sol[,3]* 0 + modN$Sol[,4] * 0  + modN$Sol[,5] * 0  + modN$Sol[,6] * 0 + modN$Sol[,7] * 0  + modN$Sol[,8] * 1  + modN$Sol[,9] * 1)+ 1),2)
HPDinterval(((modN$Sol[,1] * 1 + modN$Sol[,2] * 1  + modN$Sol[,3]* 0 + modN$Sol[,4] * 0  + modN$Sol[,5] * 0  + modN$Sol[,6] * 0 + modN$Sol[,7] * 0  + modN$Sol[,8] * 1 + modN$Sol[,9] * 1 )+ 1))


2head(margot)
summary(modN)


head(margot)


#zero poiss
prior.1R =
    list(R = list(V = diag(1), nu = 0.002),
         G = list(
         G1 = list( V= diag(2)*.002, nu = 3)))



summary(modR.I <- (MCMCglmm(relnar ~ 1, random =  ~idh(gen35):denom, rcov = ~units,  family = "gaussian", prior = prior.1R, data = margot )))






summary(modR <- (MCMCglmm(relnar ~   gen35 + age35c +  ethn3  + partner3   + pol3c + children3c + prayS + churchS, ~idh(gen35):denom, rcov = ~units,  family = "gaussian", prior = prior.1R, data = margot, nitt=20000 )))

summary(modR)


### PREDICTIVE MODEL
print((eta<-mean(modR$Sol[,1] * 1 + modR$Sol[,2] * 0  + modR$Sol[,3]* 0 + modR$Sol[,4] * 0  + modR$Sol[,5] * 0  + modR$Sol[,6] * 0 + modR$Sol[,7] * 0  + modR$Sol[,8] * 1 + modR$Sol[,9] * 0 )))


HPDinterval((modR$Sol[,1] * 1 + modR$Sol[,2] * 0  + modR$Sol[,3]* 0 + modR$Sol[,4] * 0  + modR$Sol[,5] * 0  + modR$Sol[,6] * 0 + modR$Sol[,7] * 0  + modR$Sol[,8] * 1 + modR$Sol[,9] * 0 ))


summary(modR)


3

summary(modR <- (MCMCglmm(rimp3 ~  gen35 + age35c +  ethn3  + partner3 +  beliefc  + relnarc + pol3c + children3c + prayS + churchS, random =  ~idh(gen35):denom, rcov = ~units,  family = "gaussian", prior = prior.1, data = margot )))


#### GT

summary(modR <- (MCMCglmm(rimp3 ~  gen35 + age35c +  ethn3  + partner3 +   pol3c + children3c + churchS), rcov = ~units, random = ~ as.factor(Religion.L2.Cats.T35),  family = "gaussian",  data = margot )))




summary(modR <- (MCMCglmm(children3 ~ as.factor(Religion.L3.Cats.T35) , rcov = ~units,  family = "gaussian",  data = margot )))
head(margot)



summary(modR <- (MCMCglmm(children3 ~ as.factor(Religion.L4.Cats.T35) , rcov = ~units,  family = "gaussian",  data = margot )))

coefplot2(modR, intercept=TRUE)


head(margot)


sgot$Religion.L2.Cats.T35




coefplot2(modR, intercept=TRUE)




set.seed(1)
nn <- 100
placements <- seq(0,10)
foo <- data.frame(disorder=sample(c(0,1),size=nn,replace=TRUE),
placement=sample(placements,size=nn,replace=TRUE),
ethnic=sample(c("White","Black","Other"),size=nn,replace=TRUE))
M <- glm(disorder~placement*ethnic, family=binomial, data=foo)
#Plot the fits:

plot(placements,
predict(M,newdata=data.frame(ethnic="White",placement=placements),type="response"),
type="o",pch=21,bg="white",ylab="",ylim=c(0,1))
points(placements,
predict(M,newdata=data.frame(ethnic="Black",placement=placements),type="response"),
type="o",pch=21,bg="black")
points(placements,
predict(M,newdata=data.frame(ethnic="Other",placement=placements),type="response"),
type="o",pch=21,col="red",bg="red")
legend(x="topleft",pch=21,col=c("black","black","red"),pt.bg=c("white","black","red"),
legend=c("White","Black","Other"))




##############################################################################################################################################################################################################################!#######################

say <- lm(church ~age35, data = relgot, na.action=na.omit)
summary(say)


rrelwork  <- completeFun(relgot, c("id", "partner3",  "gen35", "children3c", "age35c", "ethn3", "pol3c",  "beliefc",  "relnarc",  "denom", "rimp35c"))

#nrow(rrelwork)


 #gen35 + age35c +  ethn3  + partner3   + pol3c + children3c + prayS + churchS,


prior.8 =
    list(R = list(V = diag(1), nu = 0.002),
         G = list(
         G1 = list( V= diag(4)*.002, nu = 3)))






summary(fit1.I <- MCMCglmm(church ~  trait, random =  ~idh(trait:gen35):denom, rcov = ~trait:units,  family = "zapoisson", prior = prior.8, data =rrelwork))

summary(fit1.I <- MCMCglmm(church ~  trait*denom, rcov = ~trait:units,  family = "zapoisson", data =rrelwork))


coefplot2(fit1.I$Sol, intercept=TRUE)
        


ggplot(childata1, aes(rimp35, fill = denom,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(denom ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Denom"))


test <- MCMCglmm(church ~  age35c + rimp35 + as.factor(denom),  rcov = ~units,  family = "gaussian", data =rrelwork, nitt= 20000)


coefplot2(test, intercept =TRUE)

head(margot)

Fit1.1 <- MCMCglmm(church ~  trait * (gen35 + age35c +  ethn3  + partner3  + pol3c + children3c + relnarc + beliefc), random =  ~idh(trait:gen35):denom, rcov = ~trait:units,  family = "zapoisson", prior = prior.8, data =rrelwork, nitt= 20000)


coefplot2(fit1.1, intercept = TRUE)



summary(fit1.1)
plot(fit1.1$Sol)

## USE
round(mean(exp(fit1.1$Sol[,1] * 1 + fit1.1$Sol[,3] * 0  + fit1.1$Sol[,4]*0 + fit1.1$Sol[,5] * 0  + fit1.1$Sol[,6] * 0  + fit1.1$Sol[,7] * 0 + fit1.1$Sol[,8] * 0  + fit1.1$Sol[,9] * 0 +  fit1.1$Sol[,10] * 0)),2)



round(HPDinterval(exp((fit1.1$Sol[,1] * 1+ fit1.1$Sol[,3] * 0  + fit1.1$Sol[,4]*0 + fit1.1$Sol[,5] * 0  + fit1.1$Sol[,6] * 0  + fit1.1$Sol[,7] * 0 + fit1.1$Sol[,8] * 0  + fit1.1$Sol[,9] * 0 +  fit1.1$Sol[,10] * 0))),2)



## ZERO
round(mean(exp(-exp((fit1.1$Sol[,2] * 1    + fit1.1$Sol[,11]* 0 + fit1.1$Sol[,12] * 0  + fit1.1$Sol[,13] * 0  + fit1.1$Sol[,14] * 0 + fit1.1$Sol[,15] * 0  + fit1.1$Sol[,16] * 0 +fit1.1$Sol[,17] * 0 + fit1.1$Sol[,18] * 0 )))),2)

round(mean(exp(-exp((fit1.1$Sol[,2] * 1    + fit1.1$Sol[,11]*0  + fit1.1$Sol[,12] * 0 + fit1.1$Sol[,13] * 0  + fit1.1$Sol[,14] * 0 + fit1.1$Sol[,15] * 0  + fit1.1$Sol[,16] * 0 +fit1.1$Sol[,17] * 0 + fit1.1$Sol[,18] * 0)))),2)

round(HPDinterval(exp(-exp((fit1.1$Sol[,2] * 1    + fit1.1$Sol[,11]*0 + fit1.1$Sol[,12] * 0  + fit1.1$Sol[,13] * 0  + fit1.1$Sol[,14] * 0 + fit1.1$Sol[,15] * 0  + fit1.1$Sol[,16] * 0 +fit1.1$Sol[,17] * 0 + fit1.1$Sol[,18] * 0)))),2)




prior.8 =
    list(R = list(V = diag(1), nu = 0.002),
         G = list(
         G1 = list( V= diag(4)*.002, nu = 3)))




        


fit1.2 <- MCMCglmm(pray ~  trait * (gen35 + age35c +  ethn3  + partner3  + pol3c + children3c + relnarc + beliefc), random =  ~idh(trait:gen35):denom, rcov = ~trait:units,  family = "zapoisson", prior = prior.8, data =rrelwork, nitt= 50000)



fit1.2$Residual$nrt <- 2
summary(fit1.2)
#fix bug
#fit1.2$Residual$nrt <- 2

exp(.04)

# 
# women abv
# 
# is 5.70, 95% HPDinterval from 4.25 to  7.34; probability of zero prayer = .71, 95% HPDinterval from .60 to .81.
# 
# men average
# is 4.07, 95% HPDinterval from 2.85  to 5.54; probability of zero prayer = .69, 95% HPDinterval from .56 to .82.
# 
# 
# 
# women sd
# is  8.10, 95% HPDinterval from 5.91 to 10.49 ; probability of zero prayer = .46, 95% HPDinterval from .29 to .64.
# 
# 
# men sd
# is  5.78, 95% HPDinterval from 4.05  to 7.86 ; probability of zero prayer = .44, 95% HPDinterval from .25 to .64.


## USE
round(exp(mean(
    fit1.2$Sol[,1] * 1 + #intercept
    fit1.2$Sol[,3] * 1  +  #gen
    fit1.2$Sol[,4] * 0 + # age
    fit1.2$Sol[,5] * 0  + #ethn
    fit1.2$Sol[,6] * 0  +# partner
    fit1.2$Sol[,7] * 0 + #Pol
    fit1.2$Sol[,8] * 0  + #children
    fit1.2$Sol[,9] * 0 + #narc
    fit1.2$Sol[,10] *0 #belief
)
),2)



HPDinterval(exp(
    fit1.2$Sol[,1] * 1 + #intercep
    fit1.2$Sol[,3] * 1  +  #gen
    fit1.2$Sol[,4]* 0 + # age
    fit1.2$Sol[,5] * 0  + #ethn
    fit1.2$Sol[,6] * 0  +# partner
    fit1.2$Sol[,7] * 0 + #Pol
    fit1.2$Sol[,8] * 0  + #children
    fit1.2$Sol[,9] * 0 + #narc
    fit1.2$Sol[,10] *0 #belief
    ))

## USE
round(mean(exp(-exp(
    fit1.2$Sol[,2] * 1 + #intercept
    fit1.2$Sol[,11] * 1  +  #gen
    fit1.2$Sol[,12]*0 + # age
    fit1.2$Sol[,13] * 0  + #ethn
    fit1.2$Sol[,14] * 0  +# partner
    fit1.2$Sol[,15] * 0 + #Pol
    fit1.2$Sol[,16] * 0  + #children
    fit1.2$Sol[,17] * 0 + #narc
    fit1.2$Sol[,18] *0 #belief
))),2)


HPDinterval(exp(-exp(
    fit1.2$Sol[,2] * 1 + #intercept
    fit1.2$Sol[,11] * 1  +  #gen
    fit1.2$Sol[,12]* 0 + # age
    fit1.2$Sol[,13] * 0  + #ethn
    fit1.2$Sol[,14] * 0  +# partner
    fit1.2$Sol[,15] * 0 + #Pol
    fit1.2$Sol[,16] * 0  + #children
    fit1.2$Sol[,17] * 0 + #narc
    fit1.2$Sol[,18] *0 #belief
    )))





#mean(exp(-exp(fit1.2$Sol[,2] * 1   + fit1.2$Sol[,11]*0 + fit1.2$Sol[,12] * 0  + fit1.2$Sol[,13] * 0  + fit1.2$Sol[,14] * 0 + fit1.2$Sol[,15] * 0  + fit1.2$Sol[,16] * 0 +fit1.2$Sol[,17] * 0 + fit1.2$Sol[,18] * 0 )))



summary(fit1.2)


library(gmodels)

summary(ssgot$ethn3)


with(relgot, CrossTable(relnar, gen35))

with(relgot, CrossTable(church, gen35))




womsummary(fit1.1)


summary(fit1.1)
        
# SUGGESTS RELIGION IS A SIGNAL    Men go to church but less
summary(fit1.2I <- MCMCglmm(pray ~  trait, random =  ~idh(trait):denom, rcov = ~idh(trait):units,  family = "zapoisson", prior = prior.1, data =rrelwork))



summary(fit1.2 <- MCMCglmm(pray ~  trait * (gen35 + partner3 + age35c + children3   + relnarc + beliefc + rimp35c +  pol3c), random =  ~idh(trait):denom, rcov = ~idh(trait):units,  family = "zapoisson", prior = prior.1, data =rrelwork))

coefplot2(fit1.2, intercept = TRUE)


summary(fit1.2)        

#women
print(exp(eta<-mean(fit1.2$Sol[,1]+ fit1.2$Sol[,3] *1 + fit1.2$Sol[,4]*0 + fit1.2$Sol[,5] + fit1.2$Sol[,6] + fit1.2$Sol[,7] + fit1.2$Sol[,8]+  fit1.2$Sol[,9] +  fit1.2$Sol[,10])))

HPDinterval((exp(fit1.2$Sol[,1]+ fit1.2$Sol[,3] *1 + fit1.2$Sol[,4]*0 + fit1.2$Sol[,5] + fit1.2$Sol[,6] + fit1.2$Sol[,7] + fit1.2$Sol[,8]+  fit1.2$Sol[,9] +  fit1.2$Sol[,10])))



mean((exp(-exp(fit1.2$Sol[,2]+ fit1.2$Sol[,11] *1 + fit1.2$Sol[,12]*0 + fit1.2$Sol[,13]*0 + fit1.2$Sol[,14]*0 + fit1.2$Sol[,15]*0 + fit1.2$Sol[,16]*0+  fit1.2$Sol[,17]*0 +  fit1.2$Sol[,18]*0))))


HPDinterval((exp(-exp(fit1.2$Sol[,2]+ fit1.2$Sol[,11] *1 + fit1.2$Sol[,12]*0 + fit1.2$Sol[,13]*0 + fit1.2$Sol[,14]*0 + fit1.2$Sol[,15]*0 + fit1.2$Sol[,16]*0+  fit1.2$Sol[,17]*0 +  fit1.2$Sol[,18]*0))))




########################
########################CHILDREND 

##### # USSEE ######## !!!  1


# prior for poisson
#prior.02  =
   # list(R = list(V = diag(2), nu = 1.02, fix = 2))

prior.02  =
    list(R = list(V = diag(1), nu = 0.002, fix = 2),
             G = list(
                 G1 = list(V = diag(2)*.002, nu = 3)))


############################################################################################### #############################################################

summary(modF.1 <- (MCMCglmm(children3 ~ trait * ( gen35  + age35c + ethn3  + rels35 + pol3c), random = ~idh(trait):gen35,  rcov = ~trait:units, data = sagot, family = "zapoisson",  nitt=13000, prior=  prior.02,) ) )


summary(modF.1 <- (MCMCglmm(children3 ~ trait * ( gen35  + age35c + ethn3  + rels35 + pol3c), random = ~idh(trait):gen35,  rcov = ~trait:units, data = sagot, family = "zapoisson",  nitt=13000, prior=  prior.02,) ) )




coefplot2(modF.1, intercept =TRUE)

1+1

summary(modF.1)



mean(eta <- exp(modF.1$Sol[,1] + modF.1$Sol[,3]*1 + modF.1$Sol[,4] +  modF.1$Sol[,5]*0  +  modF.1$Sol[,6]*1 + modF.1$Sol[,7]*0 ) )

HPDinterval(eta)

mean(eta <- exp(-exp(modF.1$Sol[,2] + modF.1$Sol[,8]*1 + modF.1$Sol[,9] +  modF.1$Sol[,10]  +  modF.1$Sol[,11]*1 + modF.1$Sol[,12]*0) ))

HPDinterval(eta)

#####


#zero poiss
prior.1 =
    list(R = list(V = diag(2), nu = 0.002, fix = 1),
         G = list(
         G1 = list( V= diag(2)*.002, nu = 3)))

prior.3 =
    list(R = list(V = diag(2), nu = 0.002, fix = 1),
         G = list(
         G1 = list( V= diag(2)*.002, nu = 3),
         G2 = list( V= diag(2) * .002, nu = 3)))



### US structure
#zero poiss
prior.4 =
    list(R = list(V = diag(1), nu = 0.002, fix = 1),
         G = list(
         G1 = list( V= diag(2)*.5, nu = 3)))



## Binary
prior.12 =
    list(R = list(V = 1, nu = 0.02, fix = 1),
         G = list(
         G1 = list( V= 1, nu = .02)))

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}


#####USE

#childata  <- completeFun(relgot, c("id", "partner3", "employ",  "gen35", "prayS", "churchS",  "age35c", "ethn3", "pol3c",  "beliefc",  "denom"))

childata  <- completeFun(relgot, c("id",  "gen35", "prayS", "churchS",  "age35c",  "denom"))


nrow(childata)
dim(childata)
# "relnarc", "rimp35c",


### WOW
#interceptonly
library(MCMCglmm)
library(coefplot2)

#summary(fit.callrel.Int<- MCMCglmm(children3 ~  trait, random =  ~idh(trait:gen35):denom, rcov = ~idh(trait):units,  family = "zapoisson", prior = prior.2, data = childata))

#summary(fit.callrel<- MCMCglmm(children3 ~  trait * (age35c + ethn3 + church +  relnarc + beliefc + rimp35c +  pol3c +gen35*prayS +gen35*churchS), random =  ~idh(trait):denom, rcov = ~ idh(trait):units, prior = prior.1,  family = "zapoisson",  data = childata))


#################################################################################
#################################################################################


#### Use
rcov=~us(trait:at(sex, "M")):units+us(trait:at(sex, "F")):units

# strong prior
prior.3a =
    list(R = list(V = diag(1), nu =.05),
         G = list(
         G1 = list( V= diag(4)*.002, nu = 5)))

prior.3 =
    list(R = list(V = diag(1), nu =.02),
         G = list(
         G1 = list( V= diag(4)*.002, nu = 5)))



prior.4 =
    list(R = list(V = diag(1), nu =.05),
         G = list(
         G1 = list( V= diag(2)*.02, nu =3)))



#summary(modFC <- (MCMCglmm(children3 ~  trait * (gen35 + rim35 + beliefc + relnarc + prayS + churchS)), random = ~idh(trait:gen35):denom,  rcov = ~idh(trait):units,  data = childata, family = "poisson",   prior = prior.3,) ) )


hist(childata2$prayS)
childata2 <- subset(childata, prayS <2.5)



summary(glm(children3~ age35c  + ethn3 + pol3c +  gen35 * prayS  + gen35 * churchS, data = childata, na.action=na.omit, family="poisson"))


completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}


#####USE

#childata  <- completeFun(relgot, c("id", "partner3", "gen35", "prayS", "churchS",  "age35c", "ethn3",  "denom"))


nrow(childata)


#summary(modP <- (MCMCglmm(children3 ~ age35c +  ethn3 +  gen35 * prayS  + gen35 * churchS, random = ~idh(gen35):denom,  rcov=~trait:units,  data = childata, family = "poisson", prior=prior.4,   nitt = 50000) ) )





M <- glm(children3 ~  age35c  + gen35 * prayS  + gen35 * churchS, family ="poisson", data = childata, na.action=na.omit) 

                                        #test1 <- glm(rimp ~  age35 * denom, family ="poisson", data = childata, na.action=na.omit) 
summary(test1)

t1 <- predict(test1, type ="response")
plot(t1)








#summary(modF.1 <- (MCMCglmm(children3 ~ trait * ( gen35  + age35c + ethn3  + rels35 + pol3c), random = ~idh(trait):gen35,  rcov = ~trait:units, data = sagot, family = "zapoisson",  nitt=13000, prior=  prior.02,) ) )


library(glmmADMB)
#cAge +  lPincome + Euro + cEdu + Partner + Parent + cPolo + cSocdes * Reli, random = ~ 1|NZDep, data = smatap, family = "nbinom")

childata3  <- completeFun(relgot, c("id", "partner3", "gen35", "prayS", "churchS",  "age35c", "ethn3",   "denom","children3"))

help(glmm
     
summary(M4)

M3 <- glmmadmb(as.integer(children3) ~ age35c  + gen35*prayS  + gen35*churchS, random = ~ 1|denom, data = childata3, family = "nbinom")

     
M4 <- glmmadmb(as.integer(children3) ~ age35c  + gen35*prayS  + gen35*churchS, random = ~ 1|denom, data = childata3, family = "poisson")

summary(M4)
coefplot2(M3, intercept=TRUE)


hist(childata$children3)  

summary(modFC)
coefplot2(modF$Sol, intercept =TRUE)




plot(modF)
#

round(mean(exp((modF$Sol[,1] * 1 + modF$Sol[,3] * 0  + modFC$Sol[,4]* 0 + modF$Sol[,5] * 0  + modF$Sol[,6] * 0 +  modF$Sol[,7] * 0  +  modF$Sol[,8] *0  +  modF$Sol[,9] * 0 +  modF$Sol[,10] *0))),2)

HPDinterval(exp((modF$Sol[,1] * 1 + modF$Sol[,3] * 0  + modFC$Sol[,4]* 0 + modF$Sol[,5] * 0  + modF$Sol[,6] * 0 +  modF$Sol[,7] * 0  +  modF$Sol[,8] *0  +  modF$Sol[,9] * 0 +  modF$Sol[,10]  * 0)))


round(mean(exp(-exp(modF$Sol[,2] * 1 +  modF$Sol[,11]* 0  + modF$Sol[,12]* 0 + modF$Sol[,13] * 0  + modF$Sol[,14] * 1 +  modF$Sol[,15] * 1  +  modF$Sol[,16] * 0  +  modF$Sol[,17] * 1  +  modF$Sol[,18] * 0))),2)

HPDinterval((exp(-exp(modF$Sol[,2] * 1  + modF$Sol[,11]* 0  + modF$Sol[,12]* 0 + modFC$Sol[,13] * 0  + modF$Sol[,14] * 0 +  modF$Sol[,15] * 0   +  modF$Sol[,16] *   +  modF$Sol[,17] * 1 +  modF$Sol[,18] * 0))))



#
childata$cgen35 <- childata$gen35
contrasts(childata$cgen35) <- contr.sum(2)
str(childata$cgen35)
head(childata$cgen35)

#str(childata$gen35)
#head(childata)
exp(.04)

prior.3 =
    list(R = list(V = diag(1), nu =.02),
         G = list(
         G1 = list( V= diag(4)*.02, nu = 5)))

# No th

summary(modFa <- (MCMCglmm(children3 ~ -1 +  trait *  (gen35 + age35c + prayS  + churchS), random = ~idh(trait:gen35):denom,  rcov=~trait:units,  data = childata, family = "zapoisson",  prior = prior.3, nitt =20000) ) )




# ######################################################################



#################################################################################
#################################################################################

#summary(modFC <- (MCMCglmm(children3 ~   trait *  ( age35c + gen35*prayS  +gen35* churchS), random = ~idh(trait:gen35):denom,  rcov=~trait:units,  data = childata, family = "zapoisson",  prior = prior.3, nitt =500000) ) )

summary(modFC)
round(mean(exp((modFC$Sol[,1] * 1 + modFC$Sol[,3] * 0  + modFC$Sol[,4]* 1 + modFC$Sol[,5] * 0 + modFC$Sol[,6] * 0 +  modFC$Sol[,7] * 0  +  modFC$Sol[,8] *0))),2)

HPDinterval(exp((modFC$Sol[,1] * 1 + modFC$Sol[,3] * 0  + modFC$Sol[,4]*1  + modFC$Sol[,5] * 0  + modFC$Sol[,6] * 0 +  modFC$Sol[,7] * 0  +  modFC$Sol[,8] *0)))


round(mean(exp(-exp(modFC$Sol[,2] * 1  +  modFC$Sol[,9] * 0  +  modFC$Sol[,10] * 1 +  modFC$Sol[,11]* 0  + modFC$Sol[,12]* 0 + modFC$Sol[,13] * 0  + modFC$Sol[,14]*1 ))),2)

HPDinterval((exp(-exp(modFC$Sol[,2] * 1  +  modFC$Sol[,9] * 0  +  modFC$Sol[,10] * 1 +  modFC$Sol[,11]* 0  + modFC$Sol[,12]* 0 + modFC$Sol[,13] * 0  + modFC$Sol[,14]*0))))


newdata1 <- expand.grid(0:3, factor(0:1), 1:4)
colnames(newdata1) <- c("child", "camper", "persons")
newdata1$phat <- predict(m1, newdata1)

ggplot(newdata1, aes(x = child, y = phat, colour = factor(persons))) +
  geom_point() +
  geom_line() +
  facet_wrap(~camper) +
  labs(x = "Number of Children", y = "Predicted Fish Caught")

p <- predict(modFC)

head(p)

newdata1 <- expand.grid(factor(0:1), 1:20, 1:10)
colnames(newdata1) <- c("gen35", "pray", "church")

newdata1$children <- predict(modFC, newdata=NULL)



ggplot(newdata1, aes(x = child, y = phat, colour = factor(persons))) +
  geom_point() +
  geom_line() +
  facet_wrap(~camper) +
  labs(x = "Number of Children", y = "Predicted Fish Caught")





#################################################################################
#################################################################################
#
#!!!USE USE USE!!!
###############################################################################

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

datt <- transform(relgot,pray=scale(pray,center=TRUE,scale=FALSE), church=scale(church,center=TRUE,scale=FALSE))



datt  <- completeFun(datt, c("id",  "gen35", "pray", "church",  "age35c",  "denom", "pol3c", "ethn3"))


head(datt)


prior.4 =
    list(R = list(V = diag(1), nu =.02),
         G = list(
         G1 = list( V= diag(2)*.02, nu = 5),
         G2 = list( V= diag(4)*.02, nu = 5)))


mmodFC$Residual$nrt <- 2


odat <- datt
contrasts(odat$gen35) <- contr.sum(2)

str(odat$gen35)
head(odat$gen35)



summary(odFC <- (MCMCglmm(children3 ~   trait *  ( age35c + ethn3 + pol3c + gen35*prayS  + gen35* churchS), random = ~ idh(trait):gen35 + idh(trait:gen35):denom,  rcov= ~trait:units,  data = datt, family = "zapoisson",  prior = prior.4) ) )



summary(odFC <- (MCMCglmm(children3 ~   trait *  ( age35c + ethn3 + pol3c + at.level(gen35,1)*prayS  + at.level(gen35,1)* churchS), random = ~ idh(trait):gen35 + idh(trait:gen35):denom,  rcov= ~trait:units,  data = odat, family = "zapoisson",  prior = prior.4) ) )


summary(odFCa <- (MCMCglmm(children3 ~   trait-1 *  ( age35c + ethn3 + pol3c + at.level(gen35,1)*prayS  + at.level(gen35,1)* churchS), random = ~ idh(trait):gen35 + idh(trait:gen35):denom,  rcov= ~trait:units,  data = odat, family = "zapoisson",  prior = prior.4) ) )








coefplot2(odFC, intercept=TRUE)



summary(mmodFC <- (MCMCglmm(children3 ~   trait *  ( age35c + ethn3 + pol3c + gen35*prayS  +gen35* churchS), 
                            random = ~idh(trait:gen35):denom,  rcov= ~trait:units,  data = datt, family = "zapoisson",  prior = prior.3, nitt =100000) ) )

summary(mmodFC)
summary(mmodFC)


library(glmmADMB)
help(glmmADMB)
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

cdat <- completeFun(odat, c("children3"))

head(cdat)
nrow(cdat)                    
library(glmmADMB)
help(glmmADMB)
#completeFun <- function(data, desiredCols) {
 # completeVec <- complete.cases(data[, desiredCols])
  #return(data[completeVec, ])
#}

#fm2 <- glmmadmb(children3 ~  gen35*pray + gen35*church, random = ~ 1|denom,  data = dat, family="poisson",zeroInflation = TRUE)

fm2 <- glmmadmb(children3 ~  gen35 + age35c + ethn3 +pol3c + pray + church, random = ~ gen35|denom,  data = cdat, family="nbinom1",zeroInflation = TRUE)

            

#fm2 <- glmmadmb(children3 ~  gen35*pray + gen35*church, random = ~ 1|denom,  data = dat, family="poisson",zeroInflation = TRUE)

fm3 <- glmmadmb(children3 ~  + at.level(gen35,1) + age35c + ethn3 +pol3c + pray + church, random = ~ gen35|denom,  data = cdat, family="nbinom1",zeroInflation = TRUE)



summary(fm2)
coefplot2(fm2, intercept=TRUE)


summary(fm3)
coefplot2(fm3, intercept=TRUE)


#praysd <-  sd(relgot$pray, na.rm=T) - mean(relgot$pray, na.rm=T)
#churchsd<- sd(relgot$church, na.rm=T)-mean(relgot$church, na.rm=T)

churchsd
#praysd

# Women mean 
## Z
women average = 2.22 ; 95% HPDInterval from 1.89   to 2.56   
women probability of zero  = .55 ; 95% HPDInterval from .46  to .63

 # women sd prayer
##Z
women average =  2.36 ; 95% HPDInterval from 1.98   to  2.72
women probability of zero  =.50  ; 95% HPDInterval from .40  to .60
 

 # women sd church
##Z
women average =  2.37 ; 95% HPDInterval from 1.93   to  2.77
women probability of zero  =.61  ; 95% HPDInterval from .52  to .71


# men mean 
## Z
men average =  2.42 ; 95% HPDInterval from 2.08 to  2.72
men probability of zero  = .60 ; 95% HPDInterval from .51  to .69

 # men sd prayer
##Z
men average = 2.45  ; 95% HPDInterval from 2.02  to 2.89
men probability of zero  =.55  ; 95% HPDInterval from .39  to .69


 # men sd church
##Z
men average = 2.58  ; 95% HPDInterval from 2.11   to 2.94
men probability of zero  =.65  ; 95% HPDInterval from .57 to .73


summary(mmodFC)


coefplot2(mmodFC$Sol)

round(mean(exp((
      mmodFC$Sol[,1] * 1 + # intercept
      mmodFC$Sol[,3]* 0 + #age
      mmodFC$Sol[,4] * 0  +#ethn
                    mmodFC$Sol[,5] * 0 +  #pol
                    mmodFC$Sol[ ,6] *  1 + #gen
                    mmodFC$Sol[,7] * 0 +#pray
                    mmodFC$Sol[,8] * 0 +#ch
                    mmodFC$Sol[,9] *1  +#pray
                     mmodFC$Sol[,10] *0  #chu
                    ))),2)



HPDinterval(exp(
      mmodFC$Sol[,1] * 1 + # intercept
      mmodFC$Sol[,3]* 0 + #age
      mmodFC$Sol[,4] * 0  +#ethn
                    mmodFC$Sol[,5] * 0 +  #pol
                    mmodFC$Sol[ ,6] *  1 + #gen
                    mmodFC$Sol[,7] * 0 +#pray
                    mmodFC$Sol[,8] * 1 +#ch
                    mmodFC$Sol[,9] *0  +#pray
                     mmodFC$Sol[,10] *1  #chu
                    ))




round(mean(exp(-exp(
      mmodFC$Sol[,2] * 1 + # intercept
      mmodFC$Sol[,11]* 0 + #age
      mmodFC$Sol[,12] * 0  +#ethn
                    mmodFC$Sol[,13] * 0 +  #pol
                    mmodFC$Sol[ ,14] * 1 + #gen
                    mmodFC$Sol[,15] * 0 +#pray
                    mmodFC$Sol[,16] * 1 +#ch
                    mmodFC$Sol[,17] *0  +#pray
                     mmodFC$Sol[,18] *1  #chu
                    ))),2)




HPDinterval(exp(-exp(
      mmodFC$Sol[,2] * 1 + # intercept
      mmodFC$Sol[,11]* 0 + #age
      mmodFC$Sol[,12] * 0  +#ethn
                    mmodFC$Sol[,13] * 0 +  #pol
                    mmodFC$Sol[ ,14] * 1 + #gen
                    mmodFC$Sol[,15] * 0 +#pray
                    mmodFC$Sol[,16] * 1 +#ch
                    mmodFC$Sol[,17] * 0  +#pray
                     mmodFC$Sol[,18] *1  #chu
                    )))

#

coefplot2(mmodFC$Sol, intercept =TRUE)

library(MCMCglmm)
prior.3 =
    list(R = list(V = diag(2), nu =.02, fix =2),
         G = list(
         G1 = list( V= diag(4)*.02, nu = 5))
         )
                                        )#,
        # G2 = list( V= diag(4)*.02, nu = 5)))



summary(mmodF <- (MCMCglmm(children3 ~   
                             trait *  (gen35 + age35c + ethn3 + pol3c + prayS +  churchS),
                             random = ~idh(trait:gen35):denom,  rcov= ~idh(trait):units,  
                             data = now, 
                            family = "zapoisson",  prior = prior.3, 
                            pr=TRUE)))
#                          , nitt =100000)
) )


%
mmodF$Residual$nrt<-2

summary(mmodF)

coefplot2(test1$Sol,intercept=TRUE)

remove(newdat)
rm(newdat)
newdat <-
  with(data, expand.grid(
  prayS = seq(0,2, by = .1),
      #seq(min(pray), max(pray), length = 10),
  gen35 = levels(factor(gen35)),
  churchS= seq(0,2, by =.1),
  children3=0,
  age35c=0,
ethn3=0
  #denom= levels(factor(denom))    
  #children3= seq(min(children3), max(children3), length = 10)
))

head(newdat)
str(newdat)


nrow(newdat)
?predict.MCMCglmm

prior.3 =
    list(R = list(V = diag(1), nu =.02),
         G = list(
         G1 = list( V= diag(2)*.02, nu=3)))



test1$Residual$nrt<-2
#??mcmcglmm


X <- model.matrix(children3 ~ gen35 + age35c + ethn3 + pol3c + prayS +  churchS, newdat)

head(X)
Z <- predict(mmodF, type = 'terms', interval="confidence", level=0.95, marginal=mmodF$Random$X, newdata=newdat)

head(Z)
nrow(Z)
tail(Z)
nrow(Z)
max(Z)


nrow(Z)
nrow(X)



#X <- model.matrix(children3 ~ gen35 + age35c + ethn3 + pol3c + prayS +  churchS, newdat)
head(X)
head(Z)
nrow((X))

Z <- as.data.frame((exp(Z)))
colnames(Z) <- c("children","lwr","upr")
colnames(Z) <- "children"
dim(Z)
head(Z)
tail(Z)

dim(X)
X <- as.data.frame(X)
colnames(X)[2] <- "Gender"
# Rename Men and Women
X$Gender <- as.factor(as.character(X$Gender))
levels(X$Gender)[levels(X$Gender)=="0"] <- "Women"
levels(X$Gender)[levels(X$Gender)=="1"] <- "Men"

head(X)
head(Z)


le <- cbind(X,Z)






head(le)
#mm <- model.matrix(terms(test$Sol),newdat)
#head(X)
#newdat$children <- X#exp(X %*% test$Sol)
#head(newdat)

#guass only
#remove(newdat)
#newdat <- data.frame(
 #   newdat
  #  , plo = newdat$children-2*sqrt(pvar1)
   # , phi = newdat$children+2*sqrt(pvar1)
   # , tlo = newdat$children-2*sqrt(tvar1)
   # , thi = newdat$children+2*sqrt(tvar1)
#)


#pvar1 <- diag(mm %*% tcrossprod(vcov(fm4),mm))
#tvar1 <- pvar1+VarCorr(fm1)$denom[1]  ## must be adapted for more complex models


#plot confidence
head(le)
g2 <- ggplot(le, aes(x=prayS, y=children, color=Gender))+ geom_point()


g2
(position="identity", alpha=.3)+ geom_pointrange(aes(ymin = lwr, ymax = upr) ,alpha=.3, position="jitter")+ theme_bw()


g2

#g0 + scale_x_continuous(limits = c(0, 20)) +scale_y_continuous(limits=c(0,10))

g0

g1 <- g0+stat_sum(aes(group=pray,color=Gender), alpha=.5)#http://docs.ggplot2.org/current/stat_sum.html
head(datt)
g1

#g1 <- ggplot(datme, aes(x=pray, y=children, color=Gender))+   geom_line( aes(ymin = lwr, ymax = upr)) + theme_bw()+scale_x_continuous(limits = c(0, 90))



                                        #+ opts(title="CI based on fixed-effects uncertainty ONLY")#

#g0 + geom_errorbar(aes(ymin = lwr, ymax = upr) ,alpha=.3)  
o





summary(mmodF)


churchsd
#praysd

# Women mean 
## Z
women average =  2.20 ; 95% HPDInterval from 1.90  to   2.49
women probability of zero  =  ; 95% HPDInterval from to

 # women sd prayer
##Z
women average = 2.31 ; 95% HPDInterval from 1.98 to  2.61
women probability of zero  = .49  ; 95% HPDInterval from .39  to .58
 

 # women sd church
##Z
women average = 2.34 ; 95% HPDInterval from 1.99 to 2.68  
women probability of zero  = .60 ; 95% HPDInterval from .51 to .67


# men mean 
## Z
men average = 2.39 ; 95% HPDInterval from 2.04 to 2.76  
men probability of zero  = .59 ; 95% HPDInterval from .50 to .66

 # men sd prayer
##Z
men average = 2.51 ; 95% HPDInterval from 2.13  to  2.89  
men probability of zero  = .54 ;  95% HPDInterval from .44 to .62



 # men sd church
##Z
men average = 2.55 ; 95% HPDInterval from 2.15 to  2.99
men probability of zero  = .64  ; 95% HPDInterval from .56 to .71


summary(mmodF)
coefplot(modF, intercept =TRUE)


coefplot2(odF$Sol)



round(mean(exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 1 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0#ch
                    )), 2)



HPDinterval(exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 1 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0 #ch
                    ))




round(mean(exp(-exp(
      mmodF$Sol[,2] * 1 + # intercept
      mmodF$Sol[,9]* 1 + # gen
      mmodF$Sol[,10] * 0  +#age
                    mmodF$Sol[,11] * 0 +  #ethn
                    mmodF$Sol[ ,12] *  0 + #pol
                    mmodF$Sol[,13] * 0 +#pray
                    mmodF$Sol[,14] * 0 #ch
                    ))),2)




HPDinterval(exp(-exp(
      mmodF$Sol[,2] * 1 + # intercept
      mmodF$Sol[,9]* 1 + # gen
      mmodF$Sol[,10] * 0  +#age
                    mmodF$Sol[,11] * 0 +  #ethn
                    mmodF$Sol[ ,12] *  0 + #pol
                    mmodF$Sol[,13] * 0 +#pray
                    mmodF$Sol[,14] * 0 #chch
                    )))





coefplot2(odF$Sol)


#### DIFFS 
wbase <- (exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 0 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0#ch
                    ))

mbase <- (exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 1 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0#ch
                    ))

dbase <- wbase - mbase

 dbas<- mean(round(dbase,2))


whp <- HPDinterval(exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 0 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0 #ch
                    ))

summary(mmodF$Sol)


mhp <- HPDinterval(exp(
      mmodF$Sol[,1] * 1 + # intercept
      mmodF$Sol[,3]* 1 + # gen
      mmodF$Sol[,4] * 0 + #age
                    mmodF$Sol[,5] * 0 +  #ethn
                    mmodF$Sol[ ,6] *  0 + #pol
                    mmodF$Sol[,7] * 0 + #pray
                    mmodF$Sol[,8] * 0 #ch
                    ))

p <- whp-mhp
p
head(p)
dbas
summary(mmodF)
coefplot2(mmodF$VCV, intercept=TRUE)



round(mean(exp(-exp(
      mmodF$Sol[,2] * 1 + # intercept
      mmodF$Sol[,9]* 1 + # gen
      mmodF$Sol[,10] * 0  +#age
                    mmodF$Sol[,11] * 0 +  #ethn
                    mmodF$Sol[ ,12] *  0 + #pol
                    mmodF$Sol[,13] * 0 +#pray
                    mmodF$Sol[,14] * 0 #ch
                    ))),2)




HPDinterval(exp(-exp(
      mmodF$Sol[,2] * 1 + # intercept
      mmodF$Sol[,9]* 1 + # gen
      mmodF$Sol[,10] * 0  +#age
                    mmodF$Sol[,11] * 0 +  #ethn
                    mmodF$Sol[ ,12] *  0 + #pol
                    mmodF$Sol[,13] * 1 +#pray
                    mmodF$Sol[,14] * 0 #chch
                    )))




summary(fm4)


remove(newdat)
newdat <- expand.grid(
    pray=c(0,1,2,3),
    gen35 = c("m","f"),
    church=c(1:8),
    age35c=0,
    pol
    children3 =0)





###############



#plot(modFC$Sol)

logit(3.16)

round(mean(exp(-exp(mmodFC$Sol[,2]*1 + mmodFC$[,15] * 1)))


prior.3 =
    list(R = list(V = diag(1), nu =.02),
         G = list(
         G1 = list( V= diag(4)*.002, nu = 5)))


summary(modFCa <- (MCMCglmm(children3 ~  age35c  +  ethn3 +   gen35*prayS  + gen35*churchS, random = ~idh(gen35):denom,  rcov= ~units,  data = childata, family = "threshold",   prior = prior.3,  nitt =13000) ) )



predicted <- predict(mmodF, data=newdata, type="terms", interval="confidence", level=0.90)
%

#hist(predicted)



#childata$csgen35 <- contr.sum(childata$gen35)







## Check zeros expected
table(relgot$children3 ==0)

nrow(relgot)
e
312/1003



table(relgot$church ==0)
774/1630

ch <- as.data.frame(na.exclude(relgot$church))
nrow(ch)


ppois(0, mean(relgot$church, na.rm=TRUE))


                                        # expect 12%

#summary(modF <- (MCMCglmm(children3 ~ trait *( gen35 *(prayS + church)), random = ~idh(trait:gen35):denom,  rcov = ~idh(trait):units,  data = childata, family = "zipoisson",   prior = prior.zc) ) )



#summary(modF <- (MCMCglmm(children3 ~ trait *( gen35 *(prayS + church)), random = ~idh(trait:gen35):denom,  rcov = ~idh(trait):units,  data = childata, family = "zipoisson",   prior = prior.zc) ) )



#exp(-.11)


summary(modF)


coefplot2(modFC, intercept=TRUE)







c2 <- (16 * sqrt(3)/(15 * pi))^2

div <- sqrt(1 + c2)



HPDinterval(eta <- exp(fit1.2$Sol[,1] * 1 + fit1.2$Sol[,3] * 1 + fit1.2$Sol[,4]*0 + fit1.2$Sol[,5] * 1  + fit1.2$Sol[,6] * 0  + fit1.2$Sol[,7] * 0 + fit1.2$Sol[,8] * 0  + fit1.2$Sol[,9] * 1 +  fit1.2$Sol[,10] * 0))

div

HPDinterval(ppois(0, ((fit1.2$Sol[,2] * 1   + fit1.2$Sol[,8]* 1 + fit1.2$Sol[,9] * 0  + fit1.2$Sol[,10] * 0  + fit1.2$Sol[,11] * 0 + fit1.2$Sol[,12] * 0/div))))











# not sensible:
#summary(fit.callrel3<- MCMCglmm(children3 ~  trait * (age35c + gen35 + ethn3 +  relnarc + beliefc + rimp35c +  pol3c + prayS + churchS), random = ~idh(trait):gen35 + idh(trait):denom, rcov = ~ idh(trait):units, prior = prior.3,  family = "zapoisson",  data = childata))



coefplot2(modF.rel$Sol, intercept = TRUE)
coefplot2(modF.rel$VCV, intercept = TRUE)



### ????


#summary(modF.relo <- (MCMCglmm(children3 ~  trait *( gen35 + age35c +  ethn3 +  pol3c  + rimp35c + relnarc + beliefc + prayS + churchS), random = ~idh(trait:gen35):denom,  rcov = ~idh(trait):units,  data = childata, family = "zapoisson",   prior = prior.3) ) )



#summary(fit.relo)

#coefplot2(modF.relo$Sol, intercept = TRUE)
#coefplot2(modF.relo$VCV, intercept = TRUE)

summary(test1)




###########
# PLOT
############################


#### GRAPH POSITIVE ONLY

## Scale
datt <- transform(relgot,pray=scale(pray,center=TRUE,scale=FALSE), church=scale(church,center=TRUE,scale=FALSE))



datt  <- completeFun(datt, c("id",  "gen35", "pray", "church",  "age35c",  "denom", "pol3c", "ethn3"))


dato <- subset(datt, children3 >0)

summary(test1)
%%

test1 <- MCMCglmm(children3 ~ trait*( age35c + ethn3 + pol3c + gen35 * pray + gen35* church), data=datt, rcov= ~idh(trait):units, family = "zapoisson", nitt=20000)


test1 <- MCMCglmm(children3 ~  age35c + gen35 + pray + church, data=datt, family = "poisson", nitt=100000)

head(relgot)

#fix bug in MCMCglmm
test1$Residual$nrt<-2

summary(test1)

coefplot2(test1$Sol,intercept=TRUE)

remove(newdat)
newdat <- expand.grid(
    prayS=seq(0:3),
    gen35 = c("m","f"),
    churchS=seq(0:3),
    age35c=0,
    children3 =0)

head(newdat)
str(newdat)


nrow(newdat)

Z <- predict.MCMCglmm(mmodF, type = 'terms', interval="confidence", level=0.95)
head(Z)
tail(Z)
nrow(Z)
max(Z)


nrow(Z)
nrow(X)

head(prdicted)


X <- model.matrix(children3 ~ age35c + gen35 + pray + church, datt)
head(X)





#X <- model.matrix((children3 ~  age35c + gen35 + pray + church), datt)
head(X)
head(Z)
nrow((X))

Z <- as.data.frame((exp(Z)))
colnames(Z) <- c("children","lwr","upr")
dim(Z)
head(Z)
tail(Z)

dim(X)
X <- as.data.frame(X)
colnames(X)[3] <- "Gender"
# Rename Men and Women
X$Gender <- as.factor(as.character(X$Gender))
levels(X$Gender)[levels(X$Gender)=="0"] <- "Women"
levels(X$Gender)[levels(X$Gender)=="1"] <- "Men"

head(X)
head(Z)


datme <- cbind(X,Z)





head(datme)

#mm <- model.matrix(terms(test$Sol),newdat)
#head(X)
#newdat$children <- X#exp(X %*% test$Sol)
#head(newdat)

#guass only
#remove(newdat)
#newdat <- data.frame(
 #   newdat
  #  , plo = newdat$children-2*sqrt(pvar1)
   # , phi = newdat$children+2*sqrt(pvar1)
   # , tlo = newdat$children-2*sqrt(tvar1)
   # , thi = newdat$children+2*sqrt(tvar1)
#)


#pvar1 <- diag(mm %*% tcrossprod(vcov(fm4),mm))
#tvar1 <- pvar1+VarCorr(fm1)$denom[1]  ## must be adapted for more complex models


#plot confidence

g0 <- ggplot(datme, aes(x=pray, y=children, color=Gender))+ geom_point(position="jitter", alpha=.3)+geom_pointrange(aes(ymin = lwr, ymax = upr) ,alpha=.3, position="jitter")+ theme_bw()



g0 + scale_x_continuous(limits = c(0, 20)) +scale_y_continuous(limits=c(0,10))

g1 <- g0+stat_sum(aes(group=pray,color=Gender), alpha=.5)#http://docs.ggplot2.org/current/stat_sum.html
head(datt)
g1

#g1 <- ggplot(datme, aes(x=pray, y=children, color=Gender))+   geom_line( aes(ymin = lwr, ymax = upr)) + theme_bw()+scale_x_continuous(limits = c(0, 90))



                                        #+ opts(title="CI based on fixed-effects uncertainty ONLY")#

#g0 + geom_errorbar(aes(ymin = lwr, ymax = upr) ,alpha=.3)  
o







 + guides(fill=guide_legend(title="Gender"))



childata1 <- childata
head(childata)

dpois(mean(childata$children3, na.rm=TRUE), (nrow(childata))
)

# Rename Men and Women
levels(childata1$gen35)[levels(childata1$gen35)=="0"] <- "Women"

levels(childata1$gen35)[levels(childata1$gen35)=="1"] <- "Men"


library(ggplot2)
p1 <- ggplot(childata1, aes(children3, fill = gen35,  na.rm=TRUE)) + geom_histogram() + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender")) +  theme(legend.position="none")


library(ggplot2)

p2 <- ggplot(childata1, aes(pray, fill = gen35,  na.rm=TRUE)) + geom_histogram(width = 1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))


P3 <- ggplot(childata1, aes(church, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))



library(ggplot2)
ggplot(childata1, aes(belief, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))


library(ggplot2)
ggplot(childata1, aes(rimp3, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))


library(ggplot2)
ggplot(childat1, aes(relnar, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))

### GRID
head(childata)

head(childata1)
colnames(childata1)[6] <- "Children"
colnames(childata1)[22] <- "Church"
colnames(childata1)[23] <- "Pray"


library(ggplot2)


p1 <- ggplot(childata1, aes(Children, fill = gen35,  na.rm=TRUE)) + geom_histogram() + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender")) +  theme(legend.position="none")


p1

p2 <- ggplot(childata1, aes(Pray, fill = gen35,  na.rm=TRUE)) + geom_histogram(width = 1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))  +  theme(legend.position="none")

p2
p3 <- ggplot(childata1, aes(Church, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))

ggplot(childata1, aes(Church, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(denom ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))
p3

ggplot(childata1, aes(rimp35, fill = denom,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(denom ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Denom"))


ggplot(childata1, aes(rimp35, fill = denom,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(denom ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Denom"))




grid <- grid.arrange(p1, p2, p3, nrow=1)
gridd <- grid + labs(title = "Zero inflation")



ggplot(childata1, aes(belief, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))

ggplot(childata1, aes(relnar, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))


ggplot(childata1, aes(rimp35, fill = gen35,  na.rm=TRUE)) + geom_histogram(width=1) + facet_grid(gen35 ~ ., margins = TRUE, scales = "free") + guides(fill=guide_legend(title="Gender"))


head(childata)
### What would a possion look like


a
datamat
hist(datamat$col2)

ggplot(datamat, aes(x=col1,y=col2)) + stat_bin)# + facet_grid(col1 ~ ., margins = TRUE, scales = "free") 


library(arm)
library(lme4)



data <- relgot

head(newdata)
str(newdata)

ddata <- transform(data,pray=scale(pray,center=TRUE,scale=FALSE), church=scale(church,center=TRUE,scale=FALSE))

mean(ddata$pray,na.rm=TRUE)



summary(mod1 <-  lm(children3 ~  gen35 * pray + gen35*church, data = newdata, na.action=na.omit))



head(relgot)


library(lme4)
library(ggplot2) # Plotting


# Create for modelling overdispersion
ddata$obs_effect<-1:nrow(ddata)

mdata <- (subset(relgot, children3>0))
mdata$obs_effect<-1:nrow(mdata)

nrow(mdata)

head(relgot)
dat <- completeFun(relgot,c("denom","gen35","prayS","pray","churchS","church", "children3"))

fm4 <- glmer(
    formula = children3 ~ age35c + gen35 + pray + church +(1|denom)
    , data = mdata, na.action=na.omit, family="poisson"
)


summary(fm4)

remove(newdat)
newdat <- expand.grid(
    pray=c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),
    gen35 = c("m","f"),
    church=c(1:20),
    age35c=0,
    children3 =0)

head(newdat)
str(newdat)

rm(mm)



mm <- model.matrix(terms(fm4),newdat)
head(mm)
newdat$children <- exp(mm %*% fixef(fm4))
head(newdat)

#guass only
#remove(newdat)
#newdat <- data.frame(
 #   newdat
  #  , plo = newdat$children-2*sqrt(pvar1)
   # , phi = newdat$children+2*sqrt(pvar1)
   # , tlo = newdat$children-2*sqrt(tvar1)
   # , thi = newdat$children+2*sqrt(tvar1)
#)


#pvar1 <- diag(mm %*% tcrossprod(vcov(fm4),mm))
#tvar1 <- pvar1+VarCorr(fm1)$denom[1]  ## must be adapted for more complex models


#plot confidence
g0 <- ggplot(newdat, aes(x=pray, y=children, colour=gen35))+
    geom_smooth() + geom_point(position="jitter") #+ opts(title="CI based on fixed-effects uncertainty ONLY")#
g0

#g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
#    opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
#g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
   # opts(title="CI based on FE uncertainty + RE variance")

#summary(mod1)


%
relgot$obs_effect<-1:nrow(relgot)
head(relgot)











dat <- completeFun(relgot,c("denom","gen35","prayS","pray","churchS","church", "children3", "gen35"))

rm(newdat)
newdat <- expand.grid(
    pray=seq(0,2,by=.1),
    gen35c=c("Male","Female"),
    church=seq(0,2,by=.1),
    children3=0,
    age35c=0,
    pol3c=0
)

mm <- model.matrix(terms(fm1),newdat)
head(mm)
nrow(mm)

newdat$children <-exp( mm %*% fixef(fm1))
head(newdat)

pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$denom[1]  ## must be adapted for more complex models
newdat <- data.frame(
    newdat
    , plo = newdat$children-2*sqrt(pvar1)
    , phi = newdat$children+2*sqrt(pvar1)
    , tlo = newdat$children-2*sqrt(tvar1)
    , thi = newdat$children+2*sqrt(tvar1)
)


#plot confidence
g0 <- ggplot(newdat, aes(x=pray, y=children, colour=gen35c))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
    opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
    opts(title="CI based on FE uncertainty + RE variance")

summary(fm2)


head(relgot)


library(glmmADMB)
help(glmmADMB)
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

head(datcom)
datm<- completeFun(dat, c("ethn3", "age35c", "pol3c", "lincome3c"))


#fm2 <- glmmadmb(children3 ~  gen35*pray + gen35*church, random = ~ 1|denom,  data = dat, family="poisson",zeroInflation = TRUE)


datm <- completeFun(datm, c(  "lincome3c", "age3c", "ethn3", "pol3c", "prayS", "churchS", "denom", "gen35"))


#datcom$ageF <-datcom$age35c
nrow(datm)

fm2 <- glmmadmb(children3 ~ gen35 + age3c + ethn3 + pol3c + lincome3c+  prayS + churchS, random = ~ gen35|denom,
                data = datm,
                #family="poisson",
                family="nbinom1",
                zeroInflation = TRUE)




# family="nbinom1"
summary(fm2)
coefplot2(fm2, intercept=TRUE)
nrow(dat)




#3
#model.matrix(terms(fm2))

rm(newdat)
newdat <- expand.grid(
    prayS=seq(0,2,by=.1),
    #prayS=0,
    #gen35=c(0,1),
    gen35=.5,
    ethn3 = 0,
    churchS=0,
    #churchS=seq(0,2,by=.1),
    children3=0,
    ageF=0,
    #age35c=0,
    pol3c=0
)

head(newdat)


#newdat$children <- exp(mm %*% fixef(fm2))
head(newdat)



#library(glmmADMB)


## design matrix (fixed effects)

mm <- model.matrix(delete.response(terms(fm2)),newdat)
head(mm)
nrow(mm)

## linear predictor (for GLMMs, back-transform this with the
##  inverse link function (e.g. plogis() for binomial, beta;
##  exp() for Poisson, negative binomial
newdat$children3 <- mm %*% exp(fixef(fm2))
rm(predvar)
head(newdat)

predvar <- diag(mm %*% vcov(fm2) %*% t(mm))
newdat$SE <- sqrt(predvar)
head(newdat)
head(predvar)
head(fm2$alpha^2)
newdat$SE2 <- sqrt(predvar+fm2$alpha^2)



pvar1 <- diag(mm %*% tcrossprod(vcov(fm2),mm))
tvar1 <- pvar1+VarCorr(fm2)$denom[1]  ## must be adapted for more complex models
#remove(newdata)
newdat <- data.frame(
    newdat
   , plo = newdat$children -2*sqrt(pvar1)
   , phi = newdat$children +2*sqrt(pvar1)
    , tlo = newdat$children -2*sqrt(tvar1)
    , thi = newdat$children +2*sqrt(tvar1)
)
head(newdat)

remove(g0)
rm(g1)
rm(g2)
#plot confidence
# colour=gen35
g0 <- ggplot(newdat, aes(x=prayS, y=children3)) + geom_point(position="", alpha=.3)+theme_bw()
g0

g2 <- ggplot(newdat, aes(x=prayS, y=children3, colour=gen35))+geom_smooth(se=FALSE) + geom_point(position="jitter", alpha=.3)+theme_bw()
g2


g3 <- ggplot(newdat, aes(x=churchS, y=children3, colour=gen35))+geom_smooth(se=FALSE) + geom_point(position="jitter", alpha=.3)+theme_bw()
g3


#plot confidence
g0 <- ggplot(newdat, aes(x=churchS, y=children3, colour=gen35)) + geom_point(position="", alpha=.3)+theme_bw()
g0

g2 <- ggplot(newdat, aes(x=churchS, y=children3, colour=gen35))+geom_smooth(se=FALSE) + geom_point(position="jitter", alpha=.3)+theme_bw()
g2


g3 <- ggplot(newdat, aes(x=churchS, y=children3, colour=gen35))+geom_smooth(se=FALSE) + geom_point(position="jitter", alpha=.3)+theme_bw()
g3



summary(fm2)



exp(.83)

library(ggplot2) 
#pd <- position_dodge(width=0.4)
theme_update(theme_bw())
rm(g0)

rm(g1)
rm(g2)

head(mm)
g0 <- ggplot(newdat,aes(x=prayS,y=children3))

g0 <- g0 + geom_point(alpha=.3) +geom_line(alpha=.3) +  geom_linerange(data=newdat,
                    aes(ymin=children3-2*SE,ymax=children3+2*SE), alpha=.9)+ theme_bw()

g0

## prediction intervals 
g2 <- g0 + geom_linerange(data=newdat,
                    aes(ymin=children3-2*SE2,ymax=children3+2*SE2, alpha=.2)) + theme_bw()


g2

head(mm)#,colour=factor(gen35
g0 <- ggplot(newdat,aes(x=churchS,y=children3))
             
g0 + geom_point(alpha=.3)+ geom_linerange(data=newdat,
                    aes(ymin=children3-2*SE,ymax=children3+2*SE), alpha=.9)+ theme_bw()


## prediction intervals 
g2 <- g1 + geom_linerange(data=newdat,
                    aes(ymin=children3-2*SE2,ymax=children3+2*SE2)) + theme_bw()

g2




summary(fm2)
stargazer(fixef(fm2))

#g0 <- ggplot(newdat, aes(x=church, y=children, colour=children))+geom_point(position= "jitter")

                                                 


library(stargazer)
library(texreg)

texreg(fm2)


#g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
    opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
#g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
 #   opts(title="CI based on FE uncertainty + RE variance")

summary(fm2)



remove(g0)



## make prediction data frame
#newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female"))
## design matrix (fixed effects)
mm <- model.matrix(delete.response(terms(fm2)),newdat)
## linear predictor (for GLMMs, back-transform this with the
##  inverse link function (e.g. plogis() for binomial, beta;
##  exp() for Poisson, negative binomial
newdat$children <- mm %*% fixef(fm2)

predvar <- diag(mm %*% vcov(fm2) %*% t(mm))
newdat$SE <- sqrt(predvar) 
newdat$SE2 <- sqrt(predvar+fm2$alpha^2)
#(Probably overly complicated) ggplot code:

    head(dat)
    
library(ggplot2) 
pd <- position_dodge(width=0.4)
theme_update(theme_bw())
g0 <- ggplot(dat,aes(x=pray,y=children3,colour=gen35))+
    stat_sum(alpha=0.2,aes(size=..n..))+
    scale_size_continuous(breaks=1:4,range=c(2,5))
g0
g1 <- g0+geom_line(data=newdat,position=pd)+
      geom_point(data=newdat,shape=17,size=3,position=pd)
g1 + geom_linerange(data=newdat,
                    aes(ymin=children-2*SE,ymax=children+2*SE), position=pd)



## prediction intervals 
g1 + geom_linerange(data=newdat,
                    aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)


data <- na.omit(relgot)

newdat <-
  with(data, expand.grid(
  pray = seq(min(pray), max(pray), length = 10),
  Gen35 = levels(factor(gen35)),
  church= seq(min(church),max(church), length=5),
  children3= seq(min(children3), max(children3), length = 10)
))

#NZDepT35
fm2 <- glmmadmb(children3 ~ gen35 + prayS + churchS, random = ~ gen35|denom,  data = dat, family="nbinom1",zeroInflation = TRUE)

fm2 <- lmer(children3 ~ gen35 + prayS + churchS + 1|gen35,  data = dat, family="poisson")
summary(fm2)

rm(newdat)
newdat <-
  with(data, expand.grid(
  prayS = seq(0,2, by = .1),
      #seq(min(pray), max(pray), length = 10),
  gen35 = levels(factor(gen35)),
  churchS= seq(0,2, by =.1),
  children3=0,
  age35c=0
  #denom= levels(factor(denom))    
  #children3= seq(min(children3), max(children3), length = 10)
))


head(newdat)

newdat$childrenborn <- stats::predict(fm2, newdata=newdat, type = "response")


#plot(fm2$residuals, fm2$fitted.values)

library(ggplot2)
qplot(prayS,childrenborn, data=newdat, colour=gen35) + geom_point()+theme_bw()
qplot(churchS,childrenborn, data=newdat, colour=gen35) + geom_point()+theme_bw()
#qplot(time, arousal, data=grid, colour=factor(role)) + geom_line(data=grid)

#Aerr <- stats::predict(fm2, newdata=newdat,  se = TRUE)
#Aerr <- stats::predict(fm2, newdata=newdat,  se = TRUE)
#head(Aerr)




%
#stat="identity"






####################################################################################################################################################################################################################################################################################################################################






### WOMEN CHILDREN WHO PRAY




print(fit.callrel)

#
#

######################################################################################
######################################################################################      
######################################################################################   
       
t3_0 <-  as.data.frame(read.csv(file = "/Users/josephbulbulia/Dropbox/0R/DATA/2013/NZAVS/T_3_O/NZAVST32011.csv", header = TRUE, sep =","))    #NZAVS T3 (2011) Mplus Base Dataset.txt  

       
head(t3_0)
t3dat<-data.frame(t3_0)
head(t3dat)
dim(t3dat)

#Replace SPSS -9999 with "NA" so that R performs properly
t3dat[t3dat == "-9999" ] = NA
head(t3dat)

# ReliIdT3  !Religion.Identification.T3
# ReliGdT3  !Religion.Believe.God.T3 (1 yes, 0 no)
# ReliSpT3  !Religion.Believe.Spirit.T3 (1 yes, 0 no)
# RelCatT3  !Religion.Believe.Cats.T3
#           !(1 = God & Spirit, 2 = God - No Spirit,
#           !3 = Spirit - No God, 4 = Atheists & Agnostics) 
rtdat3<-subset(t3dat, Religious.T3== 1)

dim(rtdat3)

t3church<- as.data.frame(as.numeric(rtdat3[ ,"Religion.Church.T3"]))
dim(t3church)

t3type<-as.data.frame(as.factor(rtdat3[ ,"Religion.L3.Cats.T3"]))
       summary(t3type)
       
       rtdat3$w.Gender.T3
       
t3male<-as.data.frame(as.factor(rtdat3[ ,"Gender.T3"])) 

       
t3eth<- as.data.frame(as.factor(rtdat3[ ,"EthL2.Euro.T3"]))       
       
t3id<-as.data.frame(as.factor(rtdat3[ ,"Questionnaire.Num"]))
       
t3age<-as.data.frame(as.numeric(rtdat3[,"Age.T3"]))
       
t3strength<-as.data.frame(as.numeric(rtdat3[,"Religion.Identification.T3"])) 
      colnames(t3strength)<-c("t3st")
       library(arm)

      t3strength$t3st<-rescale(t3strength$t3st)
       
       dim(t3strength)
       
       
t3age<-scale(t3age, center =TRUE, scale=FALSE)

t3reldata<-as.data.frame(cbind(t3id,t3church,t3type,t3male,t3eth,t3age,t3strength))     
       
       
dim(t3reldata)    
head(t3reldata)
       
t3work<-na.exclude(t3reldata)
dim(t3work)
colnames(t3work)<-c("id","church","denom","male","eth","age","strength")
       head(t3work)
       str(t3work)
       t3work$male
# ReliChT3  !Religion.Church.T3
# CharitT3  !CharityDonate.T3
# EthEurT3  !EthL2.Euro.T3 (1 NZ Euro, 0 other)
      
       
#################################################
#      model  t3
#################################################
#################################################
prior_fit_t3c = 
  list(R = list(V = diag(2), nu = 0.002, fix = 2),
       G = list(
         G1 = list( V= diag(2)*.002, nu = 3)))


#Note idh fixes covariance to zero

###LOOKS GOOD
library(MCMCglmm)
library(coefplot2)
       
       exp(-0.139121)
       
fit_7_t3c <- MCMCglmm(church ~ trait - 1  + at.level(trait, 1):male 
                      + at.level(trait, 1):strength + at.level(trait, 1):age + at.level(trait, 1):eth, random =  ~idh(trait):denom, data=t3work, rcov = ~idh(trait):units, prior =prior_fit_t3c, family = "zapoisson", verbose = TRUE, nitt=50000)
 t4

       
       
fit_7_t3cI <- MCMCglmm(church ~ trait - 1  + at.level(trait, 1):male + + at.level(trait, 1):strength + at.level(trait, 1):age + at.level(trait, 1):eth, random =  ~idh(trait):denom, data=t3work, rcov = ~idh(trait):units, prior =prior_fit_t3c, family = "zipoisson", verbose = TRUE)

       

summary(fit_7_t3cI)
coefplot2(fit_7_t3cI$Sol, intercept=TRUE)
coefplot2(fit_7_t3cI$VCV, intercept=TRUE)
plot(fit_7_t3cI)
       

fit_7_t3JH <- MCMCglmm(church ~ trait * (male + age + strength + eth), random = ~trait:denom, data=t3work, rcov = ~trait:units, family = "zapoisson", verbose = TRUE)
       

       
summary(fit_7_t3JH)
coefplot2(fit_7_t3JH$Sol, intercept=TRUE)
coefplot2(fit_7_t3JH$VCV, intercept=TRUE)
plot(fit_7_t3c)
       
exp(0.140778)

##################################################################
##################################################################



## ZERO ALTERED(SAME BUT PREDICTIN NON-ZEROS: GENERALLY BETTER)
# See Hadfield course notes p. 108, note use family "hupoisson"  # WORKS BETTER
#This time with prayer   DO NOT USE THE HURDLE ASSUMPTIOSN PROB NOT MET


prior_fit_8 = 
  list(R = list(V = diag(2), nu = 0.002, fix = 2),
       G = list(
         G1 = list( V= diag(2)*.002, nu = 3),
         G2 = list( V= diag(2)*.002, nu = 3)))


# note that we fit for both men and women

fit_8 <- MCMCglmm(pray ~ trait - 1 + at.level(trait, 1):rage + at.level(trait, 1:2):relgen + at.level(trait, 1):belief + at.level(trait, 1):reth, random = ~us(trait):relid + idh(trait):denom, data=relwork, rcov = ~idh(trait):units, prior = prior_fit_8, family = "hupoisson", verbose = TRUE)


DIC(fit_8)
summary(fit_8)
coefplot2(fit_8$Sol, intercept=TRUE)
coefplot2(fit_8$VCV, intercept=TRUE)
plot(fit_8)



# Checks p. 107 course notes

c2 <- (16 * sqrt(3)/(15 * pi))^2
HPDinterval(plogis(fit_8$Sol[, 2]/sqrt(1 + c2)))

HPDinterval(ppois(0, exp(fit_8$Sol[, 1] + 0.5 *  fit_8$VCV[, 1])))
c2 <- (16 * sqrt(3)/(15 * pi))^2
HPDinterval(ppois(0, exp(fit_8$Sol[, 1] + fit_8$Sol[,4] + 0.5 * fit_8$VCV[, 1])))




#
##
#
#



DIC(fit_2)
summary(fit_2)
coefplot2(fit_2$Sol, intercept=TRUE)
Scoefplot2(fit_2$VCV, intercept=TRUE)
plot(fit_2)
autocorr(fit_2)


# head(data)
# 
prior.
a = list(
  R = list(V =diag(4), nu=1.002),
  G = list(
    G1 = list( V= diag(4), nu = 5 ),
    G2 = list( V= diag(4), nu=5)))

 DIC(mod_A)
 summary(mod_A)
 coefplot2(mod_A, intercept=TRUE)
 coefplot2(mod_A$VCV, intercept=TRUE)




mod_A<- MCMCglmm(cbind(rimp,church,pray,script) ~ trait - 1, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.a)




mod_B<- MCMCglmm(cbind(rimp,church,pray,script) ~ trait - 1 + trait:rage + trait:relgen + trait:belief + trait:reth, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.a)


plot(mod_B$Sol)
#
DIC(mod_B)
summary(mod_B)
coefplot2(mod_B$Sol, intercept=TRUE)
coefplot2(mod_B$VCV, intercept=TRUE)


#SAME as b BUT WITH INTERACTION of belief and gender
mod_B_2<- MCMCglmm(cbind(rimp,church,pray,script) ~ trait - 1 + trait:rage + trait:relgen * trait:belief + trait:reth, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.a)




DIC(mod_B_2)
summary(mod_B_2)
coefplot2(mod_B_2, intercept=TRUE)
coefplot2(mod_B_2$VCV, intercept=TRUE)



#
# model but with poisson substitute relnar for rimp
mod_C<- MCMCglmm(cbind(relnar,church,pray,script) ~ trait - 1 + trait:rage + trait:relgen + trait:belief + trait:reth, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.a)

DIC(mod_C)
summary(mod_C)
coefplot2(mod_$Sol, intercept=TRUE)
coefplot2(mod_C$VCV, intercept=TRUE)


#SAME as C BUT WITH INTERACTION of belief and gender

mod_C_2<- MCMCglmm(cbind(relnar,church,pray,script) ~ trait - 1 + trait:rage + trait:relgen * trait:belief + trait:reth, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.a)

DIC(mod_C_2)
summary(mod_C_2)
coefplot2(mod_C_2$Sol, intercept=TRUE)
coefplot2(mod_C_2$VCV, intercept=TRUE)






# Same model as C but with close to god 

prior.b = list( 
  R = list(V = diag(5), nu=1.002),
  G = list(  
    G1 = list( V= diag(5)*.05, nu = 6 ),
    G2 = list( V= diag(5)*.05, nu = 6)))



mod_D<- MCMCglmm(cbind(relnar,rgolcon,church,pray,script) ~ trait - 1 + trait:rage + trait:relgen + trait:belief, random = ~us(trait):relid + idh(trait):denom,  family= c("gaussian","gaussian","poisson","poisson","poisson"), rcov= ~us(trait):units, data=relwork, prior=prior.b)


summary(mod_D)
coefplot2(mod_D$Sol, intercept=TRUE)
coefplot2(mod_D$VCV, intercept=TRUE)




prior.b = list( 
  R = list(V =diag(5), nu=1.002),
  G = list(  
  	G1 = list( V= diag(5)*.05, nu = 5 ),
    G2 = list( V= diag(5)*.05, nu=5)))


mod_A <- MCMCglmm(cbind(quest,intrinsic,pextrinsic,sextrinsic) ~ trait - 1 + gen, random = ~us(trait):id + idh(trait):denom,  family=rep("gaussian",4), rcov= ~us(trait):units, data=work, prior=prior.a)


mod_B <- MCMCglmm(cbind(harm,fairness,loyalty,authority,sanctity) ~ trait - 1 + gen, random = ~us(trait):id  + idh(trait):denom,  family=rep("gaussian",5), rcov= ~us(trait):units, data=work, prior = prior.b)


mod_C <- MCMCglmm(loyalty ~ gen, family = "gaussian", data=work)
mod_D <- MCMCglmm(cbind(pextrinsic,sextrinsic) ~ trait-1 + gen, rcov=~us(trait):units, family = rep("gaussian",2), data=work)


coefplot2(mod_C$VCV, intercept=T)
coefplot2(mod_C$Sol, intercept=T)

coefplot2(mod_D$VCV, intercept=T)
coefplot2(mod_D$Sol, intercept=T)




prior = list( 
  R = list(V =diag(4), nu=.02),
	G = list(	
		G1 = list( V= diag(4) * .5, nu = 5 ),
		G2 = list( V= diag(4) * .5, nu = 5 ), 
		G3 = list( V= diag(4) * .5, nu = 5 ),
    G4 = list( V= diag(4) * .5, nu = 5 ),
    G5 = list( V= diag(4) * .5, nu = 5 ) ) )


prior2 = list( 
  R = list(V =diag(4), nu=1, fix = 1),
  G = list(	
		G1 = list( V= diag(4) * .5, nu = 5 ),
		G2 = list( V= diag(4) * .5, nu = 5 ), 
		G3 = list( V= diag(4) * .5, nu = 5 ),
    G4 = list( V= diag(4) * .5, nu = 5 ),
    G5 = list( V= diag(4) * .5, nu = 5 ) ) )

prior3 = list( 
  R = list(V =diag(4), nu=5),
  G = list(  
		G1 = list( V= diag(4) * .5, nu = 5 ))
)



mod_1 <- MCMCglmm(cbind(quest,intrinsic,pextrinsic,sextrinsic) ~ trait - 1 + trait:harm + trait:fairness + trait:loyalty + trait:authority + trait:sanctity, random = ~us(trait):id, family=rep("gaussian",4), rcov= ~us(trait):units, data=work, prior = prior3)

coefplot2(mod_1$VCV, intercept=T)
coefplot2(mod_1$Sol, intercept=T)



#Just the fixed effects
mod_1 <- MCMCglmm(cbind(quest,intrinsic,pextrinsic,sextrinsic) ~ trait - 1 + trait:harm + trait:fairness + trait:loyalty + trait:authority + trait:sanctity,  family=rep("gaussian",4), rcov= ~us(trait):units, data=work)

summary(mod_1)

#
coefplot2(mod_1$VCV, intercept=T)
coefplot2(mod_1$Sol, intercept=T)


library(MCMCglmm)
library(coefplot2)
head(data)

prior.a = list( 
  R = list(V =diag(4), nu=4, fix =1),
  G = list(	
		G1 = list( V= diag(4) * .5, nu = 5 ))
  )


mod_0 <- MCMCglmm(cbind(quest,intrinsic,pextrinsic,sextrinsic) ~ trait - 1, random = ~us(trait):id,  family=rep("gaussian",4), rcov= ~us(trait):units, data=work, prior =prior.a)

coefplot2(mod_0$VCV, intercept=T)


prior = list( 
  R = list(V =diag(4), nu=.02),
	G = list(	
		G1 = list( V= diag(4) * .5, nu = 5 ),
		G2 = list( V= diag(4) * .5, nu = 5 ), 
		G3 = list( V= diag(4) * .5, nu = 5 ),
    G4 = list( V= diag(4) * .5, nu = 5 ),
    G5 = list( V= diag(4) * .5, nu = 5 ) ) )


prior2 = list( 
  R = list(V =diag(4), nu=1, fix = 1),
  G = list(	
		G1 = list( V= diag(4) * .5, nu = 5 ),
		G2 = list( V= diag(4) * .5, nu = 5 ), 
		G3 = list( V= diag(4) * .5, nu = 5 ),
    G4 = list( V= diag(4) * .5, nu = 5 ),
    G5 = list( V= diag(4) * .5, nu = 5 ) ) )



mod_1 <- MCMCglmm(cbind(quest,intrinsic,pextrinsic,sextrinsic) ~ trait - 1, random = ~us(trait):harm + us(trait):fairness + us(trait):loyalty + us(trait):authority+ us(trait):sanctity, family=rep("gaussian",4), rcov= ~us(trait):units, data=work, prior = prior)








summary(mod_1)

coefplot2(mod_1$Sol, intercept = TRUE)
coefplot2(mod_1$VCV, intercept = TRUE)
plot(mod_1$VCV)
#
plot(mod_1$Sol)


```


RQue1T35  !RelQuest01.T35
RQue2T35  !RelQuest02.T35
RQue3T35  !RelQuest03.T35
RInt1T35  !RelIntrin01.T35
RInt2T35  !RelIntrin02.T35
RInt3T35  !RelIntrin03.T35
RExP1T35  !RelExtrnPers01.T35
RExP2T35  !RelExtrnPers02.T35
RExP3T35  !RelExtrnPers03.T35
RExS1T35  !RelExtrnSoc01.T35
RExS2T35  !RelExtrnSoc02.T35
RExS3T35  !RelExtrnSoc03.T35

MFHa1T35  !MFQ.harm01.T35
MFHa2T35  !MFQ.harm02.T35
MFHa3T35  !MFQ.harm03.T35
MFHa4T35  !MFQ.harm04.T35
MFHa5T35  !MFQ.harm05.T35
MFHa6T35  !MFQ.harm06.T35
MFFa1T35  !MFQ.fairness01.T35
MFFa2T35  !MFQ.fairness02.T35
MFFa3T35  !MFQ.fairness03.T35
MFFa4T35  !MFQ.fairness04.T35
MFFa5T35  !MFQ.fairness05.T35
MFFa6T35  !MFQ.fairness06.T35
MFIn1T35  !MFQ.ingroup01.T35
MFIn2T35  !MFQ.ingroup02.T35
MFIn3T35  !MFQ.ingroup03.T35
MFIn4T35  !MFQ.ingroup04.T35
MFIn5T35  !MFQ.ingroup05.T35
MFIn6T35  !MFQ.ingroup06.T35
MFAu1T35  !MFQ.authority01.T35
MFAu2T35  !MFQ.authority02.T35
MFAu3T35  !MFQ.authority03.T35
MFAu4T35  !MFQ.authority04.T35
MFAu5T35  !MFQ.authority05.T35
MFAu6T35  !MFQ.authority06.T35
MFPu1T35  !MFQ.purity01.T35
MFPu2T35  !MFQ.purity02.T35
MFPu3T35  !MFQ.purity03.T35
MFPu4T35  !MFQ.purity04.T35
MFPu5T35  !MFQ.purity05.T35
MFPu6T35  !MFQ.purity06.T35

RelIDT35  !Religion.Identification.T35


```{r fig.width=7, fig.height=6}


```



### APPENDIX INFORATION ABOUT THE NZAVS
# ======================================================================
#   
# BASE TEMPLATE FOR T3.5
# 
# ======================================================================
#   
#   DATA:
#   
#   FILE IS NZAVS T3.5 (2012) Mplus Base Dataset.txt;            
# NOBSERVATIONS ARE 4515; 
# !LISTWISE IS ON; 
# !NGROUPS = XXX;                   !Number of groups in multi-group analysis 
# 
# VARIABLE:   
#   
#   IDVARIABLE IS subnum;
# MISSING ARE ALL (-9999);
# !TSCORES = TSCR035;               !defines TSCORES for time-varying growth models
# !AUXILIARY = (m) AuxT4;           !specifies missing data covariate in LGM 
# !AUXILIARY = GenT35 (DCAT)(DCON); !specifies auxiliary variables in a mixture model
# !WEIGHT IS STATE_VAR;             !sample weight correction factor (integer weight).
# !CLUSTER = CAU06T351;             !specifies Census Area Units as a Level 2 unit in MRCM
# !BETWEEN = STATE_VAR;             !specifies BETWEEN (Level 2) variables in MRCM
# !WITHIN = STATE_VAR;              !specifies WITHIN (Level 1) variables in MRCM
# !CENSORED ARE STATE_VAR(b);       !used to test a Censor-Inflated Regression Model (bi)
# !GROUPING IS STATE_VAR;           !(1 = Group1 2 = Group 2, etc) used in multi-group
# 
# NAMES ARE 
# 
# subnum    !Questionnaire.Num
# TSCR035   !TSCORE count from Day 0 (where 0 = 30th June 2009)
# MBU06T35  !MBU.2006.T35
# MBU10T35  !MBU.2010.T35
# CAU06T35  !CAU.2006.T35
# MMICET35  !Completed.MMMICE.T35 (1 yes, 0 no)
# PIWBST35  !Completed.PIWBS.T35 (1 yes, 0 no)
# 
# GenT35    !Gender.T35 (0 female, 1 male)
# AgeT35    !Age.T35
# NZDepT35  !NZDep Percentile Score T3 (1 affluent -- 10 deprived)
# EthEuT35  !EthL2.Euro.T3 (1 NZ Euro, 0 other)
# EthMaT35  !EthL2.Maori.T3 (1 Maori, 0 other)
# EthPaT35  !EthL2.Pacific.T3 (1 Pacific, 0 other)
# EthAsT35  !EthL2.Asian.T3 (1 Asian, 0 other)
# EthOtT35  !EthL2.AllOther.T3 (1 other eth, 0 previously defined)
# ReligT35  !Religious.T35 (1 yes, 0 no)
# RelDsT35  !Religion.Denoms.T35
# !(1 Anglican, 2 Catholic, 3 Presbyterian, Christian NFD, 
#   !5 Christian Other, 6 Other Religion).
# RelL2T35  !Religion.L2.Cats.T35 
# !(-1 Outside Scope/Unreported, 0 Non-Religious, 
#   !1 Christian, 2 Other Religion)
# RelL3T35  !Religion.L3.Cats.T35 
# !(-2 Outside Scope, -1 Unreported, 0  No Religion, 1  Anglican, 2  Catholic,
#   !3  Presbyterian, 4  Methodist, 5  Baptist, 6  Mormon, 7  Ratana/Ringatu, 
#   !8  Evan/Pentecostal, 9  Christian NFD, 10 Christian Other,  11 Buddhist, 
#   !12 Hindu, 13 Other Religion)
# RelL4T35  !Religion.L4.Cats.T35 (Coding of all distinct groups at Level 4)
# !(-2 Outside scope, -1 Unreported, 0 No religion, 1 Anglican, 
#   !2 Apostolic, 3 Assemblies of God, 4 Baha'i, 5 Baptist, 6 Born Again, 
#   !7 Brethren, 8 Buddhist, 9 Catholic, 10 Christian Other, 11 Christian NFD, 
#   !12 Church of Christ, 13 Congregationist, 14 Destiny Church, 15 Evan/BA/Fund, 
#   !16 Freemason, 17 Hindu, 18 Jehova's Witness, 19 Jewish, 20 Lutheran, 
#   !21 Maori Other Religion, 22 Methodist, 23 Mormon, 24 Muslim, 25 New Age, 
#   !26 New Life, 27 Orthodox, 28 Other Religion, 29 Pagan, 30 Pentecostal, 
#   !31 Presbyterian, 32 Quaker, 33 Ratana, 34 Reformed, 35 Ringatu, 
#   !36 Salvation Army 37 Scientology, 38 Seventh Day Adventist, 39 Sikh 
#   !40 Spiritual, 41 Union, 42 Wicca).
# PRelgT35  !ParentsReligious.T35 (1 yes, 0 no)
# NZDpRT35  !NZDep Raw Score T3
# CityT35   !CityRegion.T3 (City Region Codes for 2006 areas)
# !(1 Auckland; 2 Bay of Plenty; 3 Canterbury; 4 Gisborne; 5 Hawkes Bay
#   !6 Manawatu-Wanganui; 7 Marlborough; 8 Nelson; 9 Northland; 10 Otago
#   !11 Southland; 12 Tasman; 14 Waikato; 15 Wellington; 16 West Coast)
# TerrtT35  !Territorial Area units for 2006 territorial authorities
# !1 Far North District, 2 Whangarei District,3 Kaipara District, 
# !4 Rodney District 5 North Shore City, 6 Waitakere City, 7 Auckland City, 
# !8 Manukau City, 9 Papakura District, 10 Franklin District, 
# !11 Thames-Coramandel District, 12 Hauraki District, 13 Waikato District, 
# !15 Matamata-Piako District, 16 Hamilton, 17 Waipa District, 18 Otorohanga 
# !District, 19 South Waikato District, 20 Waitomo District, 21 Taupo District, 
# !22 Western Bay of Plenty District, 23 Tauranga City, 24 Rotorua District, 
# !25 Whakatane District, 26 Kawerau District, 27 Opotiki District, 28 Gisborne 
# !District, 29 Wairoa District, 30 Hastings District, 31 Napier City, 
# !32 Central Hawke's Bay District, 33 New Plymouth District, 34 Stratford 
# !District, 35 South Taranaki District, 36 Ruapehu District, 37 Wanganui 
# !District, 38 Rangitikei District, 39 Manawatu District 40 Palmerston North
# !City, 41 Tararu District, 42 Horowhenua District, 43 Kapiti Coast District, 
# !44 Porirua City, 45 Upper Hutt City, 46 Lower Hutt City, 47 Wellington, 
# !48 Masterton District, 49 Carterton District, 50 South Wairarapa District, 
# !51 Tasman District, 52 Nelson City, 53 Marlborough District, 54 Kaikoura 
# !District, 55 Buller District, 56 Grey District, 57 Westland District, 
# !58 Hurunui District, 59 Waimakariri District, 60 Christchurch, 62 Selwyn 
# !District, 63 Ashburton District, 64 Timaru District, 65 Mackenzie District, 
# !66 Waimate District, 67 Chatham Islands Territory, 68 Waitaki District,
# !69 Central Otago District, 70 Queenstown-Lakes District, 71 Dunedin City
# !72 Clutha District, 73 Southland District, 74 Gore District, 
# !75 Invercargill District, 999 Area Outside Territorial Authority
# WardT35   !WardCode.T3 (Ward Codes for 2006 Areas)
# EthCaT35  !Ethnic Categories (Priority coded) T3.5
# !(1 NZ European, 2 Maori, 3 Pacific Nations, 4 Asian)
# RelIDT35  !Religion.Identification.T35
# RelGdT35  !Believe.God.T35 
# RelSpT35  !Believe.Spirit.T35 
# RelChT35  !Religion.Church.T35
# RelPrT35  !Religion.Prayer.T35
# RelScT35  !Religion.Scripture.T35
# EtL01T35  !EthL3.NZEuro.T3
# EtL02T35  !EthL3.NZer.T3
# EtL03T35  !EthL3.EuroNFD.T3
# EtL04T35  !EthL3.EuroOther.T3
# EtL05T35  !EthL3.Maori.T3
# EtL06T35  !EthL3.PacificNFD.T3
# EtL07T35  !EthL3.Samoan.T3
# EtL08T35  !EthL3.CookIsland.T3
# EtL09T35  !EthL3.Tongan.T3
# EtL10T35  !EthL3.Niuean.T3
# EtL11T35  !EthL3.Tokelauan.T3
# EtL12T35  !EthL3.Fijian.T3
# EtL13T35  !EthL3.PacificOth.T3
# EtL14T35  !EthL3.AsianNFD.T3
# EtL15T35  !EthL3.SthEastAsn.T3
# EtL16T35  !EthL3.Chinese.T3
# EtL17T35  !EthL3.Indian.T3
# EtL18T35  !EthL3.AsianOther.T3
# EtL19T35  !EthL3.MiddleEastern.T3
# EtL20T35  !EthL3.LatinAmerican.T3
# EtL21T35  !EthL3.African.T3
# EtL22T35  !EthL3.Other.T3
# EtL23T35  !EthL3.Unreported.T3
# EMP01T35  !Start Filler Variables 
# EMP02T35
# EMP03T35
# EMP04T35
# EMP05T35
# EMP06T35
# EMP07T35
# EMP08T35
# EMP09T35
# EMP10T35
# EMP11T35
# EMP12T35
# EMP13T35
# EMP14T35
# EMP15T35
# EMP16T35
# EMP17T35
# EMP18T35
# EMP19T35  !End Filler Variables
# ElectT35  !2012 General Electorate Districts
# !(1 Auckland Central, 2 Bay of Plenty, 3 Botany, 4 Christchurch Central, 
# !5 Christchurch East, 6 Clutha-Southland, 7 Coramandel, 8 Dunedin North,
# !9 Dunedin South, 10 East Coast, 11 East Coast Bays, 12 Epsom,
# !13 Hamilton East, 14 Hamilton West, 15 Helensville, 16 Hunua, 17 Hutt South
# !18 Ilam, 19 Invercargill, 20 Kaikoura, 21 Mana, 22 Mangere, 23 Manukau East, 
# !24 Manurewa, 25 Maungakiekie, 26 Mt Albert, 27 Mt Roskill, 28 Napier, 
# !29 Nelson, 30 New Lynn, 31 New Plymouth, 32 North Shore, 33 Northcote, 
# !34 Northland, 35 Ohariu, 36 Otaki, 37 Pakuranga, 38 Palmerston North, 
# !39 Papakura, 40 Port Hills, 41 Rangitata, 42 Rangitei, 43 Rimutaka, 44 Rodney, 
# !45 Rongotai, 46 Rotorua, 47 Selwyn, 48 Tamaki, 49 Taranaki-King Country, 
# !50 Taupo, 51 Tauranga, 52 Te Atatu, 53 Tukituki, 54 Waikato, 55 Waimakariri, 
# !56 Wairarapa, 57 Waitakere, 58 Waitaki, 59 Wellington Central,60 West Coast
# !-Tasman, 61 Whanganui, 62 Whangarei, 63 Wigram, 64 Hauraki-Waikato, 
# !65 Ikaroa-Rawhiti, 66 Tamaki Makaurau, 67 Te Tai Hauauru, 68 Te Tai Tokerau, 
# !69 Te Tai Tonga, 70 Waiariki)
# RFun1T35  !RelFund01.T35
# RFun2T35  !RelFund02.T35
# RFun3T35  !RelFund03r.T35
# RQue1T35  !RelQuest01.T35
# RQue2T35  !RelQuest02.T35
# RQue3T35  !RelQuest03.T35
# RInt1T35  !RelIntrin01.T35
# RInt2T35  !RelIntrin02.T35
# RInt3T35  !RelIntrin03.T35
# RExP1T35  !RelExtrnPers01.T35
# RExP2T35  !RelExtrnPers02.T35
# RExP3T35  !RelExtrnPers03.T35
# RExS1T35  !RelExtrnSoc01.T35
# RExS2T35  !RelExtrnSoc02.T35
# RExS3T35  !RelExtrnSoc03.T35
# RLoc1T35  !RelGodLOC01.T35
# RLoc2T35  !RelGodLOC02.T35
# RLoc3T35  !RelGodLOC03.T35
# RNar1T35  !RelCollectNarc01.T35
# RNar2T35  !RelCollectNarc02.T35
# RNar3T35  !RelCollectNarc03.T35
# PSop1T35  !PolSopSCORED01.T35 (1 correct, 0 error)
# PSop2T35  !PolSopSCORED02.T35 (1 correct, 0 error)
# PSop3T35  !PolSopSCORED03.T35 (1 correct, 0 error)
# PSop4T35  !PolSopSCORED04.T35 (1 correct, 0 error)
# PSop5T35  !PolSopSCORED05.T35 (1 correct, 0 error)
# PSop6T35  !PolSopSCORED06.T35 (1 correct, 0 error)
# PSop7T35  !PolSopSCORED07.T35 (1 correct, 0 error)
# PSop8T35  !PolSop08.LieDetect.T35 (1 yes, 0 no)
# PrLo1T35  !LOCpers01.T35
# PrLo2T35  !LOCpers02.T35
# PrLo3T35  !LOCpers03.T35
# PoLo1T35  !LOCpol01.T35
# PoLo2T35  !LOCpol02.T35
# PoLo3T35  !LOCpol03.T35
# HeLo1T35  !LOChealth01.T35
# HeLo2T35  !LOChealth02.T35
# HeLo3T35  !LOChealth03.T35
# EtNa1T35  !EthCollectNarc01.T35
# EtNa2T35  !EthCollectNarc02.T35
# EtNa3T35  !EthCollectNarc03.T35
# IWCm1T35  !IWAH.q1.community.T35
# IWCm2T35  !IWAH.q2.community.T35
# IWCm3T35  !IWAH.q3.community.T35
# IWCm4T35  !IWAH.q4.community.T35
# IWCm5T35  !IWAH.q5.community.T35
# IWCm6T35  !IWAH.q6.community.T35
# IWCm7T35  !IWAH.q7.community.T35
# IWCm8T35  !IWAH.q8.community.T35
# IWCm9T35  !IWAH.q9.community.T35
# IWNZ1T35  !IWAH.q1.NZers.T35
# IWNZ2T35  !IWAH.q2.NZers.T35
# IWNZ3T35  !IWAH.q3.NZers.T35
# IWNZ4T35  !IWAH.q4.NZers.T35
# IWNZ5T35  !IWAH.q5.NZers.T35
# IWNZ6T35  !IWAH.q6.NZers.T35
# IWNZ7T35  !IWAH.q7.NZers.T35
# IWNZ8T35  !IWAH.q8.NZers.T35
# IWNZ9T35  !IWAH.q9.NZers.T35
# IWAH1T35  !IWAH.q1.humanity.T35
# IWAH2T35  !IWAH.q2.humanity.T35
# IWAH3T35  !IWAH.q3.humanity.T35
# IWAH4T35  !IWAH.q4.humanity.T35
# IWAH5T35  !IWAH.q5.humanity.T35
# IWAH6T35  !IWAH.q6.humanity.T35
# IWAH7T35  !IWAH.q7.humanity.T35
# IWAH8T35  !IWAH.q8.humanity.T35
# IWAH9T35  !IWAH.q9.humanity.T35
# MFHa1T35  !MFQ.harm01.T35
# MFHa2T35  !MFQ.harm02.T35
# MFHa3T35  !MFQ.harm03.T35
# MFHa4T35  !MFQ.harm04.T35
# MFHa5T35  !MFQ.harm05.T35
# MFHa6T35  !MFQ.harm06.T35
# MFFa1T35  !MFQ.fairness01.T35
# MFFa2T35  !MFQ.fairness02.T35
# MFFa3T35  !MFQ.fairness03.T35
# MFFa4T35  !MFQ.fairness04.T35
# MFFa5T35  !MFQ.fairness05.T35
# MFFa6T35  !MFQ.fairness06.T35
# MFIn1T35  !MFQ.ingroup01.T35
# MFIn2T35  !MFQ.ingroup02.T35
# MFIn3T35  !MFQ.ingroup03.T35
# MFIn4T35  !MFQ.ingroup04.T35
# MFIn5T35  !MFQ.ingroup05.T35
# MFIn6T35  !MFQ.ingroup06.T35
# MFAu1T35  !MFQ.authority01.T35
# MFAu2T35  !MFQ.authority02.T35
# MFAu3T35  !MFQ.authority03.T35
# MFAu4T35  !MFQ.authority04.T35
# MFAu5T35  !MFQ.authority05.T35
# MFAu6T35  !MFQ.authority06.T35
# MFPu1T35  !MFQ.purity01.T35
# MFPu2T35  !MFQ.purity02.T35
# MFPu3T35  !MFQ.purity03.T35
# MFPu4T35  !MFQ.purity04.T35
# MFPu5T35  !MFQ.purity05.T35
# MFPu6T35  !MFQ.purity06.T35
# 
# !===MEASURES COMPLETED ONLY BY MAORI PARTICIPANTS===
# 
# MIwiT35   !Maori.IwiKNOW.T35 (1 yes, 0 no)
# MHapuT35  !Maori.HapuKNOW.T35 (1 yes, 0 no)
# MmareT35  !Maori.MaraeKNOW.T35 (1 yes, 0 no)
# MmalcT35  !Maori.MaraeLOCATION.T35 (1 yes, 0 no)
# MspkT35   !Maori.Speak.T35
# MtvT35    !Maori.WatchMaoriTV.T35 (1 yes, 0 no)
# MtvhrT35  !Maori.WatchMaoriTVHours.T35
# MLbkT35   !Maori.learn_books.T35 (1 yes, 0 no)
# MLmaT35   !Maori.learn_marae.T35 (1 yes, 0 no)
# MLgpT35   !Maori.learn_grandparents.T35 (1 yes, 0 no)
# MLfrdT35  !Maori.learn_friends.T35 (1 yes, 0 no)
# MLnetT35  !Maori.learn_net.T35 (1 yes, 0 no)
# MLkapT35  !Maori.learn_kapa.T35 (1 yes, 0 no)
# MLkohT35  !Maori.learn_kohanga.T35 (1 yes, 0 no)
# MLofmT35  !Maori.learn_othfam.T35 (1 yes, 0 no)
# MLtvT35   !Maori.learn_tv.T35 (1 yes, 0 no)
# MLprtT35  !Maori.learn_parents.T35 (1 yes, 0 no)
# MGME1T35  !Maori.GME01.T35
# MGME2T35  !Maori.GME02.T35
# MGME3T35  !Maori.GME03.T35
# MGME4T35  !Maori.GME04r.T35
# MGME5T35  !Maori.GME05r.T35
# MGME6T35  !Maori.GME06.T35
# MGME7T35  !Maori.GME07r.T35
# MGME8T35  !Maori.GME08r.T35
# MCEA1T35  !Maori.CEAIE01r.T35
# MCEA2T35  !Maori.CEAIE02r.T35
# MCEA3T35  !Maori.CEAIE03r.T35
# MCEA4T35  !Maori.CEAIE04.T35
# MCEA5T35  !Maori.CEAIE05.T35
# MCEA6T35  !Maori.CEAIE06.T35
# MCEA7T35  !Maori.CEAIE07.T35
# MCEA8T35  !Maori.CEAIE08r.T35
# MISC1T35  !Maori.ISC01.T35
# MISC2T35  !Maori.ISC02.T35
# MISC3T35  !Maori.ISC03.T35
# MISC4T35  !Maori.ISC04.T35
# MISC5T35  !Maori.ISC05.T35
# MISC6T35  !Maori.ISC06r.T35
# MISC7T35  !Maori.ISC07.T35
# MS1T35    !Maori.S01.T35
# MS2T35    !Maori.S02r.T35
# MS3T35    !Maori.S03.T35
# MS4T35    !Maori.S04.T35
# MS5T35    !Maori.S05.T35
# MS6T35    !Maori.S06r.T35
# MS7T35    !Maori.S07r.T35
# MS8T35    !Maori.S08.T35
# MSPC1T35  !Maori.SPC01r.T35
# MSPC2T35  !Maori.SPC02r.T35
# MSPC3T35  !Maori.SPC03r.T35
# MSPC4T35  !Maori.SPC04r.T35
# MSPC5T35  !Maori.SPC05.T35
# MSPC6T35  !Maori.SPC06r.T35
# MSPC7T35  !Maori.SPC07.T35
# MSPC8T35  !Maori.SPC08.T35
# MAB1T35   !Maori.AB01.T35
# MAB2T35   !Maori.AB02.T35
# MAB3T35   !Maori.AB03.T35
# MAB4T35   !Maori.AB04.T35
# MAB5T35   !Maori.AB05.T35
# MAB6T35   !Maori.AB06.T35
# MAB7T35   !Maori.AB07r.T35
# MAB8T35   !Maori.AB08r.T35
# MVEM1T35  !Maori.VEM01.T35
# MVEM2T35  !Maori.VEM02.T35
# MVEM3T35  !Maori.VEM03r.T35
# MVEM4T35  !Maori.VEM04r.T35
# MVEM5T35  !Maori.VEM05.T35
# MVEM6T35  !Maori.VEM06r.T35
# MVEM7T35  !Maori.VEM07r.T35
# 
# !===MEASURES COMPLETED ONLY BY PACIFIC PARTICIPANTS===
# 
# PWBF1T35  !Pacific.WBF01.T35
# PWBF2T35  !Pacific.WBF02.T35
# PWBF3T35  !Pacific.WBF03.T35
# PWBF4T35  !Pacific.WBF04.T35
# PWBF5T35  !Pacific.WBF05.T35
# PWBF6T35  !Pacific.WBF06.T35
# PWBF7T35  !Pacific.WBF07.T35
# PWBS1T35  !Pacific.WBS01.T35
# PWBS2T35  !Pacific.WBS02.T35
# PWBS3T35  !Pacific.WBS03.T35
# PWBS4T35  !Pacific.WBS04.T35
# PWBS5T35  !Pacific.WBS05.T35
# PWBS6T35  !Pacific.WBS06.T35
# PWBS7T35  !Pacific.WBS07.T35
# PPCB1T35  !Pacific.PCB01.T35
# PPCB2T35  !Pacific.PCB02.T35
# PPCB3T35  !Pacific.PCB03.T35
# PPCB4T35  !Pacific.PCB04.T35
# PPCB5T35  !Pacific.PCB05.T35
# PPCB6T35  !Pacific.PCB06r.T35
# PGME1T35  !Pacific.GME01.T35
# PGME2T35  !Pacific.GME02.T35
# PGME3T35  !Pacific.GME03.T35
# PGME4T35  !Pacific.GME04.T35
# PGME5T35  !Pacific.GME05.T35
# PREL1T35  !Pacific.Rel01.T35
# PREL2T35  !Pacific.Rel02.T35
# PREL3T35  !Pacific.Rel03r.T35
# PREL4T35  !Pacific.Rel04.T35
# PREL5T35  !Pacific.Rel05.T35
# PREL6T35  !Pacific.Rel06.T35
# PCE1T35   !Pacific.CultEff01.T35
# PCE2T35   !Pacific.CultEff02.T35
# PCE3T35   !Pacific.CultEff03.T35
# PCE4T35;  !Pacific.CultEff04r.T35
# 
# !CATEGORICAL ARE
# !NOMINAL ARE
# 
# USEVARIABLE ARE
# 
# MFHa1T35 MFHa2T35 MFHa3T35 MFHa4T35 MFHa5T35 MFHa6T35 
# MFFa1T35 MFFa2T35 MFFa3T35 MFFa4T35 MFFa5T35 MFFa6T35 
# MFIn1T35 MFIn2T35 MFIn3T35 MFIn4T35 MFIn5T35 MFIn6T35 
# MFAu1T35 MFAu2T35 MFAu3T35 MFAu4T35 MFAu5T35 MFAu6T35 
# MFPu1T35 MFPu2T35 MFPu3T35 MFPu4T35 MFPu5T35 MFPu6T35 
# RQue1T35 RQue2T35 RQue3T35 
# RInt1T35 RInt2T35 RInt3T35 
# RExP1T35 RExP2T35 RExP3T35 
# RExS1T35 RExS2T35 RExS3T35;
# 
# USEOBSERVATIONS ARE (ReligT35 eq 1);
# 
# DEFINE: 
# !=========================================================================================
# !Define statements for creating scale mean scores. These create the same mean scores 
# !(and names) as in the NZAVS T3.5 2012 SPSS Base File. As in the SPSS file, contrait items
# !have already been recoded. (These statements also provide a quick reference guide for 
# !estimating the same constructs as latent variables).
# 
# MFHAT35 = MEAN(MFHa1T35 MFHa2T35 MFHa3T35 
# MFHa4T35 MFHa5T35 MFHa6T35);                    !T35.MFQ.HARM
# MFFAT35 = MEAN(MFFa1T35 MFFa2T35 MFFa3T35
# MFFa4T35 MFFa5T35 MFFa6T35);                    !T35.MFQ.FAIRNESS
# MFINT35 = MEAN(MFIn1T35 MFIn2T35 MFIn3T35
# MFIn4T35 MFIn5T35 MFIn6T35);                    !T35.MFQ.INGROUP
# MFAUT35 = MEAN(MFAu1T35 MFAu2T35 MFAu3T35
# MFAu4T35 MFAu5T35 MFAu6T35);                    !T35.MFQ.AUTHORITY
# MFPUT35 = MEAN(MFPu1T35 MFPu2T35 MFPu3T35
# MFPu4T35 MFPu5T35 MFPu6T35);                    !T35.MFQ.PURITY
# RFUNT35 = MEAN(RFun1T35 RFun2T35 RFun3T35);                    !T35.REL.FUND
# RQUET35 = MEAN(RQue1T35 RQue2T35 RQue3T35);                    !T35.REL.QUEST
# RINTT35 = MEAN(RInt1T35 RInt2T35 RInt3T35);                    !T35.REL.INTRIN
# REXPT35 = MEAN(RExP1T35 RExP2T35 RExP3T35);                    !T35.REL.EXTR_PERS
# REXST35 = MEAN(RExS1T35 RExS2T35 RExS3T35);                    !T35.REL.EXTR_SOC
# RNARCT35 = MEAN(RNar1T35 RNar2T35 RNar3T35);                   !T35.REL.COLLNARC
# RGLOCT35 = MEAN(RLoc1T35 RLoc2T35 RLoc3T35);                   !T35.REL.GODLOC
# IWAHT35 = MEAN(IWAH1T35 IWAH2T35 IWAH3T35 IWAH4T35 IWAH5T35 
# IWAH6T35 IWAH7T35 IWAH8T35 IWAH9T35);           !T35.IWAH.Humanity
# IWNZT35 = MEAN(IWNZ1T35 IWNZ2T35 IWNZ3T35 IWNZ4T35 IWNZ5T35
# IWNZ6T35 IWNZ7T35 IWNZ8T35 IWNZ9T35);           !T35.IWAH.NewZealand
# IWCMT35 = MEAN(IWCm1T35 IWCm2T35 IWCm3T35 IWCm4T35 IWCm5T35 
# IWCm6T35 IWCm7T35 IWCm8T35 IWCm9T35);           !T35.IWAH.Community
# PRLOT35 = MEAN(PrLo1T35 PrLo2T35 PrLo3T35);                    !T35.LOC.Personal
# POLOT35 = MEAN(PoLo1T35 PoLo2T35 PoLo3T35);                    !T35.LOC.Political
# HELOT35 = MEAN(HeLo1T35 HeLo2T35 HeLo3T35);                    !T35.LOC.Health
# ETNAT35 = MEAN(EtNa1T35 EtNa2T35 EtNa3T35);                    !T35.EthnicNarcissism
# PSOPT35 = MEAN(PSop1T35 PSop2T35 PSop3T35 PSop4T35 
# PSop5T35 PSop6T35 PSop7T35);                    !T35.PolSop
# MGMET35 = MEAN(MGME1T35 MGME2T35 MGME3T35 MGME4T35 
# MGME5T35 MGME6T35 MGME7T35 MGME8T35);           !T35.Maori.GME
# MCEAT35 = MEAN(MCEA1T35 MCEA2T35 MCEA3T35 MCEA4T35   
# MCEA5T35 MCEA6T35 MCEA7T35 MCEA8T35);           !T35.Maori.CEAIE
# MISCT35 = MEAN(MISC1T35 MISC2T35 MISC3T35 MISC4T35
# MISC5T35 MISC6T35 MISC7T35);                    !T35.Maori.ISC
# MST35 = MEAN(MS1T35 MS2T35 MS3T35 MS4T35 
# MS5T35 MS6T35 MS7T35 MS8T35);                     !T35.Maori.S
# MSPCT35 = MEAN(MSPC1T35 MSPC2T35 MSPC3T35 MSPC4T35
# MSPC5T35 MSPC6T35 MSPC7T35 MSPC8T35);           !T35.Maori.SPC
# MABT35 = MEAN(MAB1T35 MAB2T35 MAB3T35 MAB4T35
# MAB5T35 MAB6T35 MAB7T35 MAB8T35);                !T35.Maori.AB
# MVEMT35 = MEAN(MVEM1T35 MVEM2T35 MVEM3T35 MVEM4T35 
# MVEM5T35 MVEM6T35 MVEM7T35);                    !T35.Maori.VEM
# PWBFT35 = MEAN(PWBF1T35 PWBF2T35 PWBF3T35 PWBF4T35 
# PWBF5T35 PWBF6T35 PWBF7T35);                    !T35.Pacific.WBF
# PWBST35 = MEAN(PWBS1T35 PWBS2T35 PWBS3T35 PWBS4T35 
# PWBS5T35 PWBS6T35 PWBS7T35);                    !T35.Pacific.WBS
# PPCBT35 = MEAN(PPCB1T35 PPCB2T35 PPCB3T35 PPCB4T35 
# PPCB5T35 PPCB6T35);                             !T35.Pacific.PCB
# PGMET35 = MEAN(PGME1T35 PGME2T35 PGME3T35 PGME4T35 
# PGME5T35);                                      !T35.Pacific.GME
# PRELT35 = MEAN(PREL1T35 PREL2T35 PREL3T35 PREL4T35 
# PREL5T35 PREL6T35);                             !T35.Pacific.REL
# PCET35 = MEAN(PCE1T35 PCE2T35 PCE3T35 PCE4T35);                !T35.Pacific.CE
# 
# MODEL: 
# 



# 
# URBAN CODING
# 
#                     IF (Terr06T1 EQ 1) THEN UrbanT1 = 0;  !1 Far North District
#                   IF (Terr06T1 EQ 2) THEN UrbanT1 = 0;  !2 Whangarei District
#                   IF (Terr06T1 EQ 3) THEN UrbanT1 = 0;  !3 Kaipara District
#                   IF (Terr06T1 EQ 4) THEN UrbanT1 = 0;  !4 Rodney District
#                   IF (Terr06T1 EQ 5) THEN UrbanT1 = 1;  !5 North Shore City
#                   IF (Terr06T1 EQ 6) THEN UrbanT1 = 1;  !6 Waitakere City
#                   IF (Terr06T1 EQ 7) THEN UrbanT1 = 1;  !7 Auckland City
#                   IF (Terr06T1 EQ 8) THEN UrbanT1 = 1;  !8 Manukau City
#                   IF (Terr06T1 EQ 9) THEN UrbanT1 = 0;  !9 Papakura District
#                   IF (Terr06T1 EQ 10) THEN UrbanT1 = 0; !10 Franklin District
#                   IF (Terr06T1 EQ 11) THEN UrbanT1 = 0; !!11 Thames-Coramandel District
#                   IF (Terr06T1 EQ 12) THEN UrbanT1 = 0; !12 Hauraki District
#                   IF (Terr06T1 EQ 13) THEN UrbanT1 = 0; !13 Waikato District
#                   IF (Terr06T1 EQ 15) THEN UrbanT1 = 0; !15 Matamata-Piako District
#                   IF (Terr06T1 EQ 16) THEN UrbanT1 = 1; !16 Hamilton City
#                   IF (Terr06T1 EQ 17) THEN UrbanT1 = 0; !17 Waipa District
#                   IF (Terr06T1 EQ 18) THEN UrbanT1 = 0; !18 Otorohanga District
#                   IF (Terr06T1 EQ 19) THEN UrbanT1 = 0; !19 South Waikato District
#                   IF (Terr06T1 EQ 20) THEN UrbanT1 = 0; !20 Waitomo District
#                   IF (Terr06T1 EQ 21) THEN UrbanT1 = 0; !21 Taupo District
#                   IF (Terr06T1 EQ 22) THEN UrbanT1 = 0; !22 Western Bay of Plenty District
#                   IF (Terr06T1 EQ 23) THEN UrbanT1 = 1; !23 Tauranga City
#                   IF (Terr06T1 EQ 24) THEN UrbanT1 = 0; !24 Rotorua District
#                   IF (Terr06T1 EQ 25) THEN UrbanT1 = 0; !25 Whakatane District
#                   IF (Terr06T1 EQ 26) THEN UrbanT1 = 0; !26 Kawerau District
#                   IF (Terr06T1 EQ 27) THEN UrbanT1 = 0; !27 Opotiki District
#                   IF (Terr06T1 EQ 28) THEN UrbanT1 = 0; !28 Gisborne District
#                   IF (Terr06T1 EQ 29) THEN UrbanT1 = 0; !29 Wairoa District
#                   IF (Terr06T1 EQ 30) THEN UrbanT1 = 0; !30 Hastings District
#                   IF (Terr06T1 EQ 31) THEN UrbanT1 = 1; !31 Napier City
#                   IF (Terr06T1 EQ 32) THEN UrbanT1 = 0; !32 Central Hawke's Bay District
#                   IF (Terr06T1 EQ 33) THEN UrbanT1 = 0; !33 New Plymouth District
#                   IF (Terr06T1 EQ 34) THEN UrbanT1 = 0; !34 Stratford District
#                   IF (Terr06T1 EQ 35) THEN UrbanT1 = 0; !35 South Taranaki District
#                   IF (Terr06T1 EQ 36) THEN UrbanT1 = 0; !36 Ruapehu District
#                   IF (Terr06T1 EQ 37) THEN UrbanT1 = 0; !37 Wanganui District
#                   IF (Terr06T1 EQ 38) THEN UrbanT1 = 0; !38 Rangitikei District
#                   IF (Terr06T1 EQ 39) THEN UrbanT1 = 0; !39 Manawatu District
#                   IF (Terr06T1 EQ 40) THEN UrbanT1 = 1; !40 Palmerston North City
#                   IF (Terr06T1 EQ 41) THEN UrbanT1 = 0; !41 Tararu District
#                   IF (Terr06T1 EQ 42) THEN UrbanT1 = 0; !42 Horowhenua District
#                   IF (Terr06T1 EQ 43) THEN UrbanT1 = 0; !43 Kapiti Coast District
#                   IF (Terr06T1 EQ 44) THEN UrbanT1 = 1; !44 Porirua City
#                   IF (Terr06T1 EQ 45) THEN UrbanT1 = 1; !45 Upper Hutt City
#                   IF (Terr06T1 EQ 46) THEN UrbanT1 = 1; !46 Lower Hutt City
#                   IF (Terr06T1 EQ 47) THEN UrbanT1 = 1; !47 Wellington City
#                   IF (Terr06T1 EQ 48) THEN UrbanT1 = 0; !48 Masterton District
#                   IF (Terr06T1 EQ 49) THEN UrbanT1 = 0; !49 Carterton District
#                   IF (Terr06T1 EQ 50) THEN UrbanT1 = 0; !50 South Wairarapa District
#                   IF (Terr06T1 EQ 51) THEN UrbanT1 = 0; !51 Tasman District
#                   IF (Terr06T1 EQ 52) THEN UrbanT1 = 1; !52 Nelson City
#                   IF (Terr06T1 EQ 53) THEN UrbanT1 = 0; !53 Marlborough District
#                   IF (Terr06T1 EQ 54) THEN UrbanT1 = 0; !54 Kaikoura District
#                   IF (Terr06T1 EQ 55) THEN UrbanT1 = 0; !55 Buller District
#                   IF (Terr06T1 EQ 56) THEN UrbanT1 = 0; !56 Grey District
#                   IF (Terr06T1 EQ 57) THEN UrbanT1 = 0; !57 Westland District
#                   IF (Terr06T1 EQ 58) THEN UrbanT1 = 0; !58 Hurunui District
#                   IF (Terr06T1 EQ 59) THEN UrbanT1 = 0; !59 Waimakariri District
#                   IF (Terr06T1 EQ 60) THEN UrbanT1 = 1; !60 Christchurch City
#                   IF (Terr06T1 EQ 62) THEN UrbanT1 = 0; !62 Selwyn District
#                   IF (Terr06T1 EQ 63) THEN UrbanT1 = 0; !63 Ashburton District
#                   IF (Terr06T1 EQ 64) THEN UrbanT1 = 0; !64 Timaru District
#                   IF (Terr06T1 EQ 65) THEN UrbanT1 = 0; !65 Mackenzie District
#                   IF (Terr06T1 EQ 66) THEN UrbanT1 = 0; !66 Waimate District
#                   IF (Terr06T1 EQ 67) THEN UrbanT1 = 0; !67 Chatham Islands Territory
#                   IF (Terr06T1 EQ 68) THEN UrbanT1 = 0; !68 Waitaki District
#                   IF (Terr06T1 EQ 69) THEN UrbanT1 = 0; !69 Central Otago District
#                   IF (Terr06T1 EQ 70) THEN UrbanT1 = 0; !70 Queenstown-Lakes District
#                   IF (Terr06T1 EQ 71) THEN UrbanT1 = 1; !71 Dunedin City
#                   IF (Terr06T1 EQ 72) THEN UrbanT1 = 0; !72 Clutha District
#                   IF (Terr06T1 EQ 73) THEN UrbanT1 = 0; !73 Southland District
#                   IF (Terr06T1 EQ 74) THEN UrbanT1 = 0; !74 Gore District
#                   IF (Terr06T1 EQ 75) THEN UrbanT1 = 0; !75 Invercargill District
# )

