###################################################
### chunk number 1: trans
###################################################
#line 111 "Effects_package.Rnw"
library(car)
library(effects)
data(Ornstein)
mod <- lm(interlocks ~ log(assets) + 
	sector + nation, data=Ornstein)
eff <- effect("log(assets)", mod, default.levels=100)
pdf("asset_eff.pdf", height=6, width=6)
print(
plot(eff)	
)
invisible(dev.off())


###################################################
### chunk number 2: eff2
###################################################
#line 141 "Effects_package.Rnw"
eff2 <- effect("log(assets)", mod, 
	default.levels=100, typical = "median")


###################################################
### chunk number 3: nmm
###################################################
#line 152 "Effects_package.Rnw"
colnames(eff2$model.matrix)


###################################################
### chunk number 4: eff3
###################################################
#line 164 "Effects_package.Rnw"
eff3 <- effect("log(assets)", mod, 
	default.levels=100, 
	given.value= c("sectorBNK" = 0, 
	"sectorCON" = 1, "sectorFIN" = 0,
	"sectorHLD" = 0, "sectorMAN" = 0, 
	"sectorMER" = 0, "sectorMIN" = 0, 
	"sectorTRN" = 0, "sectorWOD" = 0,
	"nationOTH" = 0, "nationUK" = 0, 
	"nationUS" = 0))


###################################################
### chunk number 5: op1
###################################################
#line 184 "Effects_package.Rnw"
mod <- lm(prestige ~ income + 
	education + poly(women, 2), 
	data=Prestige)
eff <- effect("poly(women, 2)", mod,
	default.levels=100)
pdf("womeff.pdf", height=6, width=6)
plot(eff)
invisible(dev.off())


###################################################
### chunk number 6: splines
###################################################
#line 209 "Effects_package.Rnw"
library(splines)
mod <- lm(prestige ~ bs(income, df=6) + 
	education + women, data=Prestige)
eff <- effect("bs(income, df=6)", mod, 
	default.levels=100)
pdf("inceff.pdf", height=6, width=6)
plot(eff)
invisible(dev.off())


###################################################
### chunk number 7: inteff1
###################################################
#line 233 "Effects_package.Rnw"
data(SLID)
mod <- lm(log(wages) ~ education + age + 
	sex*language, data=SLID)
eff <- effect("sex*language", mod)

pdf("sl1.pdf", height=6, width=6)
print(plot(eff))
invisible(dev.off())
pdf("sl2.pdf", height=6, width=6)
print(plot(eff, x.var="sex", as.table=T))
invisible(dev.off())


###################################################
### chunk number 8: both
###################################################
#line 262 "Effects_package.Rnw"
mod <- lm(log(wages) ~ education + language + 
	sex*age, data=SLID)
eff <- effect("sex*age", mod)

pdf("sl3.pdf", height=6, width=6)
print(plot(eff, as.table=T))
invisible(dev.off())
pdf("sl4.pdf", height=6, width=6)
print(plot(eff, x.var="sex", as.table=T))
invisible(dev.off())


###################################################
### chunk number 9: contin
###################################################
#line 291 "Effects_package.Rnw"
mod <- lm(log(wages) ~ sex + language + 
	education*age, data=SLID)
eff <- effect("education*age", mod)

pdf("sl5.pdf", height=6, width=6)
print(plot(eff, as.table=T))
invisible(dev.off())
pdf("sl6.pdf", height=6, width=6)
print(plot(eff, x.var="age", as.table=T))
invisible(dev.off())


###################################################
### chunk number 10: bin
###################################################
#line 322 "Effects_package.Rnw"
library(foreign)
nes <- read.dta("anes1992.dta")
mod <- glm(votedem ~ age + pid + 
	black + retnat, data=nes, family=binomial)
eff <- effect("pid", mod, 
	given.values = c(
		"pid" = 4, 
		"black" = 0, 
		"retnatsame" = 0, 
		"retnatworse" = 1))

pdf("pid1.pdf", height=6, width=6)
print(plot(eff))
invisible(dev.off())
pdf("pid2.pdf", height=6, width=6)
print(plot(eff, rescale.axis=F))
invisible(dev.off())


###################################################
### chunk number 11: ologit_ex
###################################################
#line 401 "Effects_package.Rnw"
options(useFancyQuotes=F)
library(foreign)
fh <- read.dta("fh2006_ss.dta")
library(MASS)
mod <- polr(freedomstat ~ log(pop) + gdppc10k + civ, data=fh)
summary(mod)


###################################################
### chunk number 12: pvals
###################################################
#line 418 "Effects_package.Rnw"
smod <- summary(mod)

p.val <- 2*pt(abs(smod$coef[,3]), 
  mod$df.residual, lower.tail=F) 

round(cbind(smod$coef, p.val), 4)


###################################################
### chunk number 13: geff1
###################################################
#line 434 "Effects_package.Rnw"
gdp.eff <- effect("gdppc10k", mod, 
  default.levels=100)


###################################################
### chunk number 14: geff_graph
###################################################
#line 445 "Effects_package.Rnw"
print(
plot(gdp.eff, xlab = "GDP/capita (in $10000)", 
  ylab = "Pr(Freedom Status)")
)


###################################################
### chunk number 15: geff_stacked
###################################################
#line 460 "Effects_package.Rnw"
print(
plot(gdp.eff, xlab="GDP/capita (in $10000 US)", 
  ylab = "Pr(Freedom Status)", style="stacked", 
  colors=c("gray75", "gray50", "gray25"), 
  main = "")
)


###################################################
### chunk number 16: olex
###################################################
#line 501 "Effects_package.Rnw"
plot.dat <- data.frame(
  probs = c(gdp.eff$prob), 
  status = rep(levels(fh$freedomstat), each=100), 
  lower = c(gdp.eff$lower.prob), 
  upper = c(gdp.eff$upper.prob), 
  gdppc10k = c(gdp.eff$x))


###################################################
### chunk number 17: olg1
###################################################
#line 517 "Effects_package.Rnw"
trellis.par.set(strip.background = list(col=NA))
print(
xyplot(probs ~ gdppc10k | status, data=plot.dat, 
  xlab = "GDP/capita (in $10000 US)",
  ylab = "Pr(Freedom Status)", 
  layout=c(1,3), as.table=T,
  ylim=c(-.05, 1.05),
  panel = function(x,y,subscripts){
    panel.lines(x,y,lty=1, col="black")
    panel.lines(x,plot.dat$lower[subscripts], 
      lty=2, col="black")
    panel.lines(x,plot.dat$upper[subscripts], 
      lty=2, col="black")},
  key=list(space="top", 
    lines=list(col=c("black", "black"), 
      lty=c(1,2)), 
    text=list(c("Fitted Probabilities", 
      "95% CI"))))
)


###################################################
### chunk number 18: olg1a
###################################################
#line 553 "Effects_package.Rnw"
plot(range(plot.dat$gdppc10k), 
  range(c(plot.dat$lower, plot.dat$upper)), 
  type = "n", xlab= "GDP/capita (in $10000 US)", 
  ylab = "Pr(Freedom Status)")

sp <- split(plot.dat, plot.dat$status)

lines(sp[[1]]$gdppc10k, sp[[1]]$probs, lty=1)
lines(sp[[2]]$gdppc10k, sp[[2]]$probs, lty=2)
lines(sp[[3]]$gdppc10k, sp[[3]]$probs, lty=3)
legend("topleft", c("Free", "Partly Free", "Not Free"),
  lty=c(1,3,2), inset=.01)


###################################################
### chunk number 19: civeff
###################################################
#line 582 "Effects_package.Rnw"
civ.eff <- effect(
  "civ", mod)
print(
plot(civ.eff)
)


###################################################
### chunk number 20: olciv
###################################################
#line 598 "Effects_package.Rnw"
plot.dat <- data.frame(
  probs = c(civ.eff$prob), 
  lower = c(civ.eff$lower.prob),
  upper = c(civ.eff$upper.prob), 
  status = rep(c("Not Free", 
    "Partly Free", "Free"), each=6), 
  civ = as.factor(rep(civ.eff$x[[1]], 3)))


###################################################
### chunk number 21: olg2
###################################################
#line 613 "Effects_package.Rnw"
trellis.par.set(strip.background = 
  list(col=NA))
print(
xyplot(probs ~ civ | status, 
  data=plot.dat, as.table=T, 
  ylim=c(-.05, 1.05), xlab = "", 
  ylab = "Pr(Freedom Status)", 
  layout = c(1,3), 
  panel = function(x,y,subscripts){
    panel.points(x,y,pch=16, 
      col="black", cex=1.25)
    panel.arrows(
      x, plot.dat$lower[subscripts], 
      x, plot.dat$upper[subscripts],
      col="black", angle=90, code=3)})
)


###################################################
### chunk number 22: olg2a
###################################################
#line 646 "Effects_package.Rnw"
trellis.par.set(strip.background=
  list(col=NA))
print(
xyplot(probs ~ status | civ, 
  data=plot.dat, as.table=T, 
  ylim=c(-.05, 1.05), xlab = "", 
  ylab = "Pr(Freedom Status)", 
  layout = c(2,3), 
  panel = function(x,y,subscripts){
    panel.points(x,y,pch=16, 
      col="black", cex=1.25)
    panel.arrows(
      x, plot.dat$lower[subscripts], 
      x, plot.dat$upper[subscripts],
      col="black", angle=90, code=3)})
)


###################################################
### chunk number 23: frdat
###################################################
#line 721 "Effects_package.Rnw"
library(foreign)
frdat1 <- read.dta("fra94ds_rintro.dta")
frdat <- read.dta("fra94ds_rintro.dta", convert.factors=F)
frdat$vote <- frdat1$vote
frdat$vote <- as.factor(as.character(frdat$vote))
frdat$retnat <- factor(frdat$retnat, 
  levels=1:3, labels=c("better", "same", "worse"))
rm(frdat1)

library(nnet)
multi.mod <- multinom(vote ~ lrself + retnat + age + male + 
    edulevel + notrelig + urban + demsat + occupa2 + occupa3 + 
    eusup, data=frdat)


###################################################
### chunk number 24: multisum
###################################################
#line 744 "Effects_package.Rnw"
summary(multi.mod)


###################################################
### chunk number 25: multip
###################################################
#line 756 "Effects_package.Rnw"
smulti <- summary(multi.mod)
multi.t <- smulti$coefficients/smulti$standard.errors
multi.p <- pt(abs(multi.t), 
  nrow(multi.mod$fitted.values)-multi.mod$edf, 
  lower.tail=F)
b <- round(smulti$coefficients, 3)
b[which(multi.p > .05, arr.ind=T)] <- ""


###################################################
### chunk number 26: printb
###################################################
#line 772 "Effects_package.Rnw"
noquote(b)


###################################################
### chunk number 27: lrs_eff
###################################################
#line 784 "Effects_package.Rnw"
lrs.eff <- effect("lrself", 
  multi.mod, default.levels=25)


###################################################
### chunk number 28: lrsg1
###################################################
#line 797 "Effects_package.Rnw"
plot(
plot(lrs.eff, rug=FALSE)
)


###################################################
### chunk number 29: lrs_stack
###################################################
#line 815 "Effects_package.Rnw"
plot(
plot(lrs.eff, 
  rug=FALSE, 
  style="stacked", 
  colors = 
    c("gray0", "gray25", 
    "gray50", "gray75", 
    "gray100"),
  main = "")
)


###################################################
### chunk number 30: lrsout
###################################################
#line 839 "Effects_package.Rnw"
plot.dat <- data.frame(
  prob = c(lrs.eff$prob), 
  lower = c(lrs.eff$lower.prob),
  upper = c(lrs.eff$upper.prob), 
  lrs = c(lrs.eff$x[[1]]), 
  party = rep(lrs.eff$y.levels, each=25))


###################################################
### chunk number 31: lrsg2
###################################################
#line 852 "Effects_package.Rnw"
trellis.par.set(
  strip.background = list(col=NA))
print(
xyplot(prob ~ lrs | party, 
  data = plot.dat, layout=c(1,5), 
  xlab = "Left-Right Self Placement",
  ylab = "Pr(Vote)", 
  panel = function(x,y,subscripts){
    panel.lines(x,y,col="black")
    panel.lines(x, plot.dat$lower[subscripts], 
      col="black", lty=2)
    panel.lines(x, plot.dat$upper[subscripts], 
      col="black", lty=2)}, 
  key=list(space="top", 
    lines=list(lty=c(1,2), 
      col=c("black", "black")),
    text=list(c("Fitted Probabilities", 
      "95% CI"))))
)


###################################################
### chunk number 32: lrsg3
###################################################
#line 889 "Effects_package.Rnw"
print(
xyplot(prob ~ lrs, 
  groups=party, 
  type="l", 
  data = plot.dat, 
  xlab = 
  "Left-Right Self Placement", 
  ylab = "Pr(Vote)",
  lty=c(1:5), 
  col="black", 
  lwd=2,
  key=list(space="top", 
    lines=list(
      lty=1:5, 
      col="black"), 
    text=list(
      lrs.eff$y.levels)))
)


###################################################
### chunk number 33: ret_eff
###################################################
#line 920 "Effects_package.Rnw"
ret.eff <- effect("retnat", multi.mod)
print(
plot(ret.eff, as.table=T, 
  layout=c(2,3),
  ylim = c(-.05, 1.05), 
  xlab = "Economic Perceptions", 
  ylab = "Pr(Vote)")
)


###################################################
### chunk number 34: retdat
###################################################
#line 938 "Effects_package.Rnw"
plot.dat <- data.frame(
  prob = c(ret.eff$prob), 
  lower = c(ret.eff$lower.prob),
  upper = c(ret.eff$upper.prob), 
  ret = factor(c(ret.eff$x[[1]]), 
    levels=1:3,
    labels=c("better", "same", "worse")),
  party = rep(ret.eff$y.levels, each=3))


###################################################
### chunk number 35: retg2
###################################################
#line 954 "Effects_package.Rnw"
trellis.par.set(
  strip.background = list(col=NA))
print(
xyplot(prob ~ ret | party, 
  data=plot.dat, 
  as.table=T, 
  xlab = "", ylab = "Pr(Vote)", 
  ylim= c(-.05, 1.05), 
  panel = function(x,y,subscripts){
    panel.arrows(
      x, plot.dat$lower[subscripts], 
      x, plot.dat$upper[subscripts], 
      code = 3, angle=90, length=.05)
    panel.points(x,y,
      pch=16, col="black")})
)


###################################################
### chunk number 36: retg3
###################################################
#line 987 "Effects_package.Rnw"
trellis.par.set(
  strip.background = list(col=NA))
print(
xyplot(prob ~ party | ret, 
  data=plot.dat, 
  as.table=T, 
  layout = c(3,1),
  scales=list(x=list(rot=45)),
  xlab = "", ylab = "Pr(Vote)", 
  ylim= c(-.05, 1.05), 
  panel = function(x,y,subscripts){
    panel.arrows(
      x, plot.dat$lower[subscripts], 
      x, plot.dat$upper[subscripts], 
      code = 3, angle=90, length=.05)
    panel.points(x,y,
      pch=16, col="black")})
)