---
title: "Ch12_2e_ChapterExamples"
output: pdf_document
---

#### Preparation
```{r setup, include=TRUE, message = FALSE, warning = FALSE}
require(knitr)
require(haven)
#library(probitmfx) #install.packages("probitmfx")
library(AER)
library(apsrtable)
library(tidyverse)

opts_chunk$set(echo = TRUE)
options(digits = 3)
opts_knit$set(root.dir = 
"C:/Users/baileyma/Dropbox/Bailey_RealStats_RealEconometrics_2e_Supplements_March2021/ChapterExamples_R")

# Remove objects from the previous session
rm(list = ls(all = TRUE))   
```
```{r tidy = FALSE}
dog 		= read_dta("Ch12_ChapterExample_DogPolitics.dta")	
DogProb 	= glm(obamavote~ 	dog + as.numeric(libcon), 	family = binomial(link = "probit"), na.action=na.exclude,  data=dog)
summary(DogProb)
#probitmfx(obamavote~ 	dog + as.numeric(libcon), data=dog, atmean= FALSE)
```

```{r tidy = FALSE}
## Continuous variable, observed-value approach, discrete change
## ProbFit	= fitted(DogProb)		# Missing data can mess up fitted values unless use na.exclude option c(min(ProbFit, na.rm=T), max(ProbFit, na.rm=T))
P1b		= pnorm(DogProb$coef[1] + DogProb$coef[2]*dog$dog + DogProb$coef[3]*as.numeric(dog$libcon))
P2		= pnorm(DogProb$coef[1] + DogProb$coef[2]*dog$dog + DogProb$coef[3]*(as.numeric(dog$libcon)+1))
mean(P2[is.na(dog$obamavote)==0]-P1b[is.na(dog$obamavote)==0], na.rm=T)

	# marginal effects
	Marg		= DogProb$coef[3]* dnorm(DogProb$coef[1] + DogProb$coef[2]*dog$dog + DogProb$coef[3]*as.numeric(dog$libcon))
	mean(Marg[is.na(dog$obamavote)==0], na.rm=T)		
	Marg[1:10]

## Dummy variable, observed-value approach, discrete change
P0		= pnorm(DogProb$coef[1] + DogProb$coef[2]*0 + DogProb$coef[3]*as.numeric(dog$libcon))
P1		= pnorm(DogProb$coef[1] + DogProb$coef[2]*1 + DogProb$coef[3]*as.numeric(dog$libcon))
#mean(P1[is.na(ProbFit)==0]-P0[is.na(ProbFit)==0], na.rm=T)


# Logit
DogLog 	= glm(obamavote~ 	dog + as.numeric(libcon), 	
              family = binomial(link = "logit"),  
              na.action= na.exclude, 
              data=dog)
summary(DogLog)
P1 = fitted(DogLog)
P1b = exp(DogLog$coef[1] + DogLog$coef[2]*dog$dog + DogLog$coef[3]*as.numeric(dog$libcon)) / (1+ exp(DogLog$coef[1] + DogLog$coef[2]*dog$dog + DogLog$coef[3]*as.numeric(dog$libcon)))
P2  = exp(DogLog$coef[1] + DogLog$coef[2]*dog$dog + DogLog$coef[3]*(as.numeric(dog$libcon)+1)) / (1+ exp(DogLog$coef[1] + DogLog$coef[2]*dog$dog + DogLog$coef[3]*(as.numeric(dog$libcon)+1)))
mean(P2-P1, na.rm=T)		#mean(P1b-P1, na.rm=T)
c(min(P1, na.rm=T), max(P1, na.rm=T))

```

```{r tidy = FALSE}
## Dummy variable, observed-value approach, discrete change
P0		= exp(DogLog$coef[1] + DogLog$coef[2]*0 + DogLog$coef[3]*as.numeric(dog$libcon)) / (1+ exp(DogLog$coef[1] + DogLog$coef[2]*0 + DogLog$coef[3]*as.numeric(dog$libcon)))
P1		= exp(DogLog$coef[1] + DogLog$coef[2]*1 + DogLog$coef[3]*as.numeric(dog$libcon)) / (1+ exp(DogLog$coef[1] + DogLog$coef[2]*1 + DogLog$coef[3]*as.numeric(dog$libcon)))

#LPM
DogLPM 	= lm(obamavote~ 	dog + as.numeric(libcon), 	na.action= na.exclude, data=dog)
summary(DogLPM)
max(fitted(DogLPM), na.rm=T)
min(fitted(DogLPM), na.rm=T)
	# Standardized
	DogLPMStd 	= lm(obamavote~ 	dog + scale(as.numeric(libcon)), 	na.action= na.exclude, data=dog)
	summary(DogLPMStd)

sqrt(var(as.numeric(dog$libcon), na.rm=T))
coefficients(DogLPM)[3]* sqrt(var(as.numeric(dog$libcon), na.rm=T))

table(dog$libcon)
table(as.numeric(dog$libcon))
# Fitted values
coefficients(DogLPM)[1] + 1*coefficients(DogLPM)[2] + 7*coefficients(DogLPM)[3]
coefficients(DogLPM)[1] + 0*coefficients(DogLPM)[2] + 1*coefficients(DogLPM)[3]

# Displays a table in Latex form
apsrtable(DogLPM, DogProb, DogLog, model.names= c("LPM", "Probit", "Logit"),
	coef.names=c("Constant", "Dog owner", "Ideology"), digits=3) 
```

```{r tidy = FALSE}
LibConJitter 	= as.numeric(dog$libcon) + rnorm(length(dog$libcon))/3
LibConSpan 		= seq(1, 7, by = 0.1)
ProbFitLineDog0 	= pnorm(coef(DogProb)[1] + coef(DogProb)[2]*0 + coef(DogProb)[3]*LibConSpan )
ProbFitLineDog1 	= pnorm(coef(DogProb)[1] + coef(DogProb)[2]*1 + coef(DogProb)[3]*LibConSpan )
LogFitLineDog0 	= 1/(1+exp(-(coef(DogLog)[1] + coef(DogLog)[2]*0 + coef(DogLog)[3]*LibConSpan )))
LogFitLineDog1 	= 1/(1+exp(-(coef(DogLog)[1] + coef(DogLog)[2]*1 + coef(DogLog)[3]*LibConSpan )))
LPMFitLineDog0 	= coef(DogLPM)[1] + coef(DogLPM)[2]*0 + coef(DogLPM)[3]*LibConSpan
LPMFitLineDog1 	= coef(DogLPM)[1] + coef(DogLPM)[2]*1 + coef(DogLPM)[3]*LibConSpan
```
```{r tidy = FALSE}
par(mfrow=c(1, 3), mar=c(2.0, 2, 1.3,0.5), oma=c(2, 2.5, 3., 0.3))	 # mar(south, west, north, east)
plot(as.numeric(dog$libcon), fitted(DogLPM),  ylim=c(-0.05, 1.05), xlim=c(1, 7), type = "n", xlab="", xaxt='n',  yaxt='n', ylab="", cex=0.75, lwd =0.5, col="blue")
	mtext("Probability\nvote\nfor\nObama", 	side = 3, line = -1.2, at = 0, cex = 0.8)
	mtext("Linear probability model (LPM)", 	side = 3, line = 0.2, cex = 0.75)
	axis(1, at=seq(1, 7, by=1),	labels=seq(1, 7 , by=1), cex.axis=0.85, padj=-1.)
	axis(2, at = seq(0.0, 1.0 , by=0.25), labels =  seq(0.0, 1.0 , by=0.25), 
	     tick = T, las = 1, cex.axis = .85, mgp = c(2,.7,0))
	mtext("Conservative ideology", 		side = 1, line = 1.5, cex = 0.8)
	lines(LibConSpan, LPMFitLineDog0, type="l", lwd=2.0, col="darkblue") 
	lines(LibConSpan, LPMFitLineDog1, type="l", lwd=2.0, col="cornflowerblue", lty = 2) 
	text(6.1, 0.9, "non-dog owner", col="darkblue", cex=0.75)
	text(5.85, 0.85, "dog owner", col="cornflowerblue", cex=0.75)
	lines(c(4, 5.), c(0.90, 0.90) ,type="l", lwd=1.0, col="darkblue") 
	lines(c(4, 5.), c(0.85, 0.85) ,type="l", lwd=1.0, col="cornflowerblue", lty = 2) 

plot(as.numeric(dog$libcon), fitted(DogLPM),  ylim=c(-0.05, 1.05), xlim=c(1, 7), type = "n", xlab="", xaxt='n',  yaxt='n', ylab="", cex=0.75, lwd =0.5, col="blue")
	axis(1, at=seq(1, 7, by=1),	labels=seq(1, 7 , by=1), cex.axis=0.85, padj=-1.)
	axis(2, at = seq(0.0, 1.0 , by=0.25), labels =  seq(0.0, 1.0 , by=0.25), 
	     tick = T, las = 1, cex.axis = .85, mgp = c(2,.7,0))
	mtext("Conservative ideology", 		side = 1, line = 1.5, cex = 0.8)
	mtext("Probit", 					side = 3, line = 0.2, cex = 0.75)
	lines(LibConSpan, ProbFitLineDog0,type="l", lwd=2.0, col="darkblue") 
	lines(LibConSpan, ProbFitLineDog1,type="l", lwd=2.0, col="cornflowerblue", lty = 2) 
	text(6.1, 0.9, "non-dog owner", col="darkblue", cex=0.75)
	text(5.85, 0.85, "dog owner", col="cornflowerblue", cex=0.75)
	lines(c(4, 5.), c(0.90, 0.90) ,type="l", lwd=1.0, col="darkblue") 
	lines(c(4, 5.), c(0.85, 0.85) ,type="l", lwd=1.0, col="cornflowerblue", lty = 2) 

plot(as.numeric(dog$libcon), fitted(DogLPM),  ylim=c(-0.05, 1.05), xlim=c(1, 7), type = "n", xlab="", xaxt='n',  yaxt='n', ylab="", cex=0.75, lwd =0.5, col="blue")
	axis(1, at=seq(1, 7, by=1),	labels=seq(1, 7 , by=1), cex.axis=0.85, padj=-1.)
	axis(2, at = seq(0.0, 1.0 , by=0.25), labels =  seq(0.0, 1.0 , by=0.25), 
	     tick = T, las = 1, cex.axis = .85, mgp = c(2,.7,0))
	mtext("Conservative ideology", 		side = 1, line = 1.5, cex = 0.8)
	mtext("Logit", 					side = 3, line = 0.2, cex = 0.75)
	lines(LibConSpan, LogFitLineDog0 ,type="l", lwd=2.0, col="darkblue") 
	lines(LibConSpan, LogFitLineDog1 ,type="l", lwd=2.0, col="cornflowerblue", lty = 2) 
	text(6.1, 0.9, "non-dog owner", col="darkblue", cex=0.75)
	text(5.85, 0.85, "dog owner", col="cornflowerblue", cex=0.75)
	lines(c(4, 5.), c(0.90, 0.90) ,type="l", lwd=1.0, col="darkblue") 
	lines(c(4, 5.), c(0.85, 0.85) ,type="l", lwd=1.0, col="cornflowerblue", lty = 2) 
```

## Example for computing corner
```{r tidy = FALSE}
	DogProbitExample 	= glm(obamavote~ 	as.numeric(libcon), 	family = binomial(link = "probit"),  na.action= na.exclude, data=dog)
	P1Example 		= fitted(DogProbitExample)
	Xsequence = seq(1, 7, 0.05)
	plot(Xsequence, pnorm(DogProbitExample$coef[1] + DogProbitExample$coef[2]*Xsequence), type = "l")
```

## Dogs and cats example
```{r tidy = FALSE}

## Estimate unrestricted model
DogCatUR 	= glm(obamavote~ 	dog + cat + as.numeric(libcon), 	family = binomial(link = "probit"), na.action=na.exclude,  data=dog)
summary(DogCatUR)

## Estimate restricted model
dog$DogCat	= dog$dog + dog$cat
DogCatR 	= glm(obamavote~ 	DogCat + as.numeric(libcon), 	family = binomial(link = "probit"), na.action=na.exclude,  data=dog)
summary(DogCatR)

## Calculate LR test statistic
logLik(DogCatUR)[1]
logLik(DogCatR)[1]
LR 	= 2*(logLik(DogCatUR)[1] - logLik(DogCatR)[1])		# LR	= 2*(-6726.7899 +6733.9411)
14.3024

## Calculate p value and critical value for LR test statistic
1-pchisq(LR, 1)		# p value
qchisq(0.95, 1)		# Critical value

## Wald test
library(car)
linearHypothesis(DogCatUR,c("dog = cat"),test="F")

# library(lmtest)	#waldtest(DogCatUR, DogCatR)

apsrtable(DogCatUR , DogCatR, model.names= c("UR", "R"),
	coef.names=c("Constant", "Dog owner", "Cat owner", "Ideology", "Dog owner + Cat owner" ), digits=3) # Displays a table in Latex form
	# APSRTABLE gets likelihood wrong
```

## More dog politics
```{r tidy = FALSE}
# marginal effects
	Marg		= DogProb$coef[3]* dnorm(DogProb$coef[1] + DogProb$coef[2]*dog$dog + DogProb$coef[3]*as.numeric(dog$libcon))
	mean(Marg[is.na(dog$obamavote)==0], na.rm=T)		
	Marg[1:10]

## Dummy variable, observed-value approach, discrete change
P0		= pnorm(DogProb$coef[1] + DogProb$coef[2]*0 + DogProb$coef[3]*as.numeric(dog$libcon))
P1		= pnorm(DogProb$coef[1] + DogProb$coef[2]*1 + DogProb$coef[3]*as.numeric(dog$libcon))

```

# Case Study from Real Econometrics

```{r tidy = FALSE}
library(foreign)
#ketchup	= read.dta("C:\\Users\\baileyma\\Dropbox\\Textbook_drop\\TextbookData\\Ketchup\\Ketchup_Stata12.dta")	
#save(ketchup, file = "Ch12_ChapterExample_Ketchup.RData")
load(file = "Ch12_ChapterExample_Ketchup.RData")

# Subset
ketchupExam = ketchup %>% 
  filter(BoughtKetchup ==1) %>%
  select(StoreBrand, Pdiff, "display" = dp4, "featured" = ft4)
save(ketchupExam, file = "PracticeExam_Ketchup.RData")

# LPM for core model
KetchLPM 	= lm(StoreBrand~ 	Pdiff + DisplayBrand + dp4 + FeatureBrand + ft4 +
                 hhold_income + hhold_members, 	na.action=na.exclude,
               data=ketchup[ketchup$BoughtKetchup ==1,])

# Estimate heteroscedastiticy-consistent standard errors for LPM
library(lmtest)
library(sandwich) ## 
lmtest::coeftest(KetchLPM, sandwich::vcovHC(KetchLPM, type = "HC1"))
summary(KetchLPM)

# Missing data can mess up fitted values unless use na.exclude option
LPMFit	= fitted(KetchLPM)		

# Min and max of fitted values
c(min(LPMFit), max(LPMFit))

# Probit
KetchProbit 	= glm(StoreBrand~ 	Pdiff + DisplayBrand + dp4 + FeatureBrand + ft4 +
                     hhold_income + hhold_members, 	
                   family = binomial(link = "probit"), na.action=na.exclude,
                   data=ketchup[ketchup$BoughtKetchup ==1,])
summary(KetchProbit)

#log likelihood in the probit model
logLik(KetchProbit)

# Missing data can mess up fitted values unless use na.exclude option
ProbFit	= fitted(KetchProbit)
c(min(ProbFit), max(ProbFit))

KetchLogit 	= glm(StoreBrand~ 	Pdiff + DisplayBrand + dp4 + FeatureBrand + ft4 +
                    hhold_income + hhold_members, 	
                  family = binomial(link = "logit"), na.action=na.exclude,
                  data=ketchup[ketchup$BoughtKetchup ==1,])
summary(KetchLogit)
# Missing data can mess up fitted values unless use na.exclude option
LogitFit	= fitted(KetchLogit )
c(min(LogitFit), max(LogitFit))
logLik(KetchLogit )	# Note apsrtable function provides wrong likelihood (this function produces same value as Stata)

## Continuous variable, observed-value approach, discrete change
PdiffPlus	= ketchup$Pdiff+ 1		# 0.16 is one std dev
P1b		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members )

P2		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*PdiffPlus+ KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members )
mean(P2-P1b, na.rm=T)

# marginal effects of Pdiff
Marg		= KetchProbit$coef[2]* dnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
mean(Marg, na.rm=T)		
Marg[1:10]

P2		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff+ KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + 
			KetchProbit$coef[7]*(ketchup$hhold_income + 1 ) + KetchProbit$coef[8]*ketchup$hhold_members )
mean(P2-P1b, na.rm=T)

P2		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff+ KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + 
			KetchProbit$coef[7]*ketchup$hhold_income + KetchProbit$coef[8]*(ketchup$hhold_members + 1))
mean(P2-P1b, na.rm=T)

## Dummy variable, observed-value approach, discrete change
P0dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*0  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
P1dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*1  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
mean(P1dumProbit-P0dumProbit, na.rm=T)


P0dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*0 + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
P1dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*1+ KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
mean(P1dumProbit-P0dumProbit, na.rm=T)


P0dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*0 + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
P1dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand+ KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*1 + KetchProbit$coef[6]*ketchup$ft4 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
mean(P1dumProbit-P0dumProbit, na.rm=T)


P0dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand + KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*0 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
P1dumProbit		= pnorm(KetchProbit$coef[1] + KetchProbit$coef[2]*ketchup$Pdiff + KetchProbit$coef[3]*ketchup$DisplayBrand+ KetchProbit$coef[4]*ketchup$dp4  + KetchProbit$coef[5]*ketchup$FeatureBrand + KetchProbit$coef[6]*1 + KetchProbit$coef[7]*ketchup$hhold_income + 
			KetchProbit$coef[8]*ketchup$hhold_members)
mean(P1dumProbit-P0dumProbit, na.rm=T)




# Logit
P1Logit 	= fitted(KetchLogit)
P1bLogit 	= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
P2Logit 	= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*PdiffPlus + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*PdiffPlus + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P2Logit-P1bLogit, na.rm=T)		#mean(P1b-P1, na.rm=T)	# mean(P2Logit, na.rm=T)
c(min(P1, na.rm=T), max(P1, na.rm=T))


P2Logit 	= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff+ KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*(ketchup$hhold_income + 1) +
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff+ KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*(ketchup$hhold_income +1) + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P2Logit-P1bLogit, na.rm=T)

P2Logit 	= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff+ KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*(ketchup$hhold_members+1)) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff+ KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*(ketchup$hhold_members+1)))
mean(P2Logit-P1bLogit, na.rm=T)


## Dummy variable, observed-value approach, discrete change
P0dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*0    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*0    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
P1dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*1    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*1    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P1dumLogit-P0dumLogit, na.rm=T)

P0dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*0+ KetchLogit$coef[4]*ketchup$dp4   + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*0 + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
P1dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*1 + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*1 + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P1dumLogit-P0dumLogit, na.rm=T)

P0dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4   + KetchLogit$coef[5]*0 + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*0 + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
P1dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*1 + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*1 + KetchLogit$coef[6]*ketchup$ft4 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P1dumLogit-P0dumLogit, na.rm=T)


P0dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4   + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*0 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*0 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
P1dumLogit		= exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*1 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members) / (1+ exp(KetchLogit$coef[1] + KetchLogit$coef[2]*ketchup$Pdiff + KetchLogit$coef[3]*ketchup$DisplayBrand + KetchLogit$coef[4]*ketchup$dp4    + KetchLogit$coef[5]*ketchup$FeatureBrand + KetchLogit$coef[6]*1 + KetchLogit$coef[7]*ketchup$hhold_income + 
			KetchLogit$coef[8]*ketchup$hhold_members))
mean(P1dumLogit-P0dumLogit, na.rm=T)

library(apsrtable)
apsrtable(KetchLPM, KetchProbit, KetchLogit, model.names= c("LPM", "Probit", "Logit"),
	coef.names=c("Constant", "Price difference", "Brand name displayed", "Store brand displayed", "Brand name featured", "Store brand featured", "Household income", "Household size"), digits=3) # Displays a table in Latex form

## Other analyses
# Estimate household income as a factor (instead of psuedo-continuous as above)
KetchLPMfactor 	= lm(StoreBrand~ 	Pdiff + DisplayBrand + dp4 + FeatureBrand + ft4 +
                       factor(hhold_income) + hhold_members, 	na.action=na.exclude,
                     data=ketchup[ketchup$BoughtKetchup ==1,])
summary(KetchLPMfactor)
```