###################################################
### chunk number 1: mix1
###################################################
#line 308 "MixtureModels.Rnw"
options(useFancyQuotes=F)
library(flexmix)
library(foreign)
dat <- na.omit(read.dta("~/Desktop/ICPSR_Slides/Collinearity/mixture_data.dta"))
dat$gdppc10k <- dat$gdppc/10000
model <- FLXMRglmfix(family = "gaussian", nested=list(k=c(1,1), 
	formula = c(~lgates + aclp + bnr + xrcomp + parcomp,
		 ~log_checks + polconiii + xconst + laworder + subfed)), 
	fixed = ~ gdppc10k + logpop + cwar + iwar)
out <- stepFlexmix(rep1 ~ 1 |ccode, k=2, model=model, data=dat, nrep=20)
 
 
###################################################
### chunk number 2: groupprob
###################################################
#line 330 "MixtureModels.Rnw"
probs <- out@posterior$scaled
apply(probs, 2, mean)
summary(out)
 
 
###################################################
### chunk number 3: histprob
###################################################
#line 344 "MixtureModels.Rnw"
pdf("hsitprob.pdf", heigh=6, width=6)
hist(probs[,1], xlab="Posterior Probabilities of Elections Group", main="")
invisible(dev.off())
 
 
###################################################
### chunk number 4: id
###################################################
#line 358 "MixtureModels.Rnw"
source("mixtureTools.R")
id.list <- IdentifyList(out@posterior$scaled, dat, dat$ccode, cluster=T, alpha=.05)
id.list$caseID.1
id.list$caseID.2
 
 
###################################################
### chunk number 5: refit
###################################################
#line 374 "MixtureModels.Rnw"
out.refit <- refit(out)
summary(out.refit)
 
 
###################################################
### chunk number 6: map
###################################################
#line 387 "MixtureModels.Rnw"
library(maptools)
world <- readShapePoly("~/Desktop/ICPSR_Slides/Lecture 3/WorldCountries.shp")
world$group <- NA
world$group[which(world$ccode %in% id.list$caseID.1)] <- 1
world$group[which(world$ccode %in% id.list$caseID.2)] <- 2
world$group <- factor(world$group, levels=1:2, 
	labels=c("Elections", "Checks"))
cols <- c("red", "blue")
print(spplot(world["group"], col.regions=cols))
 
 
###################################################
### chunk number 7: concom
###################################################
#line 414 "MixtureModels.Rnw"
out.a <- stepFlexmix(rep1 ~ 1 | ccode , k=2, model=model, data=dat, nrep=20, 
		concomitant = FLXPmultinom( ~ parlresp_c))
out.a.refit <- refit(out.a)
 
 
###################################################
### chunk number 8: sums
###################################################
#line 426 "MixtureModels.Rnw"
summary(out.a.refit)
out.a.refit@concomitant
 
 
###################################################
### chunk number 9: dftest
###################################################
#line 441 "MixtureModels.Rnw"
mod1 <- lm(rep1 ~lgates + aclp + bnr + xrcomp + parcomp + 
	gdppc10k + logpop + cwar + iwar, data=dat, y=T)
mod2 <- lm(rep1 ~log_checks + polconiii + xconst + laworder + 
	subfed + gdppc10k + logpop + cwar + iwar, data=dat, y=T)
smod1 <- summary(mod1)
smod2 <- summary(mod2)
il1 <- log(dnorm(mod1$y, mod1$fitted, smod1$sigma))
il2 <- log(dnorm(mod2$y, mod2$fitted, smod2$sigma))
cuts <- qbinom(c(.05,.95), size=length(mod1$y), prob=.5)
d <- sum(il1 > il2)
1-pbinom(d, length(il1), .5)