# 계량경제학 강의 R 코드

# 제16장 이항반응모형

# 16.1 선형확률모형

# 예제 16.1 여성의 경제활동 참여

#install.packages("sampleSelection")
data(Mroz87, package="sampleSelection")
names(Mroz87)

sum(Mroz87$lfp)

model <- lfp~nwifeinc+educ+exper+I(exper^2)+age+kids5+kids618
ols <- lm(model,data=Mroz87)
library(lmtest)
coeftest(ols)

library(sandwich)
coeftest(ols,vcov=vcovHC)

Mroz87$phat <- fitted(ols)
fgls <- lm(model,data=Mroz87,weights=1/(phat*(1-phat)))

summary(Mroz87$phat)

# 16.2 일반화된 선형모형

# 그림 16.1

curve(log(x/(1-x)),0,1,axes=F,lwd=3)
axis(1)
axis(2)
abline(h=0,v=c(0,1),lty=3)

# 예제 16.2 여성의 경제활동 참여에 관한 로짓과 프로빗 추정

model <- lfp~nwifeinc+educ+exper+I(exper^2)+age+kids5+kids618
logit <- glm(model,data=Mroz87,family=binomial(link="logit"))
coeftest(logit)
probit <- glm(model,data=Mroz87,family=binomial(link="probit"))
coeftest(probit)

# 16.3 맞춘값과 설명력

# 예제 16.3 로짓과 프로빗에서 확률 예측

Mroz87$phatl <- predict(logit,type="response")
Mroz87$phatp <- predict(probit,type="response")
Mroz87[1:10,c('lfp','phatl','phatp')]

plot(phatp~phatl,data=Mroz87,xlab="logit p-hat",ylab="probit p-hat")
abline(0,1,lty=2)
abline(v=0.5,h=0.5,lty=2)

Mroz87$yhatl <- as.numeric(Mroz87$phatl>.5)
with(Mroz87, mean(lfp==yhatl))
Mroz87$yhatp <- as.numeric(Mroz87$phatp>.5)
with(Mroz87, mean(lfp==yhatp))

# 예제 16.4 로짓과 프로빗에서 유사 R제곱

with(Mroz87, cor(lfp,yhatl)^2)
with(Mroz87, cor(lfp,yhatp)^2)

# 16.4 로짓과 프로빗의 관계

cbind(ols$coef, logit$coef, probit$coef)

logit$coef/probit$coef

curve(qlogis(pnorm(x)),-2,2)
abline(0,1.6,lty=2)
abline(0,1.7,lty=3)
abline(0,1.8,lty=4)

# 16.6 계수의 해석과 한계효과

# 예제 16.5 로짓과 프로빗에 의한 평균에서의 한계효과와 평균 한계효과

model <- lfp~nwifeinc+educ+exper+I(exper^2)+age+kids5+kids618
logit <- glm(model,data=Mroz87,family=binomial(link="logit"))
b1logit <- logit$coef['nwifeinc']
b1logit
Mroz87$xb <- predict(logit,type='link')
mea1 <- b1logit*dlogis(mean(Mroz87$xb))
mea1
ame1 <- b1logit*mean(dlogis(Mroz87$xb))
ame1

probit <- glm(model,data=Mroz87,family=binomial(link="probit"))
b1probit <- probit$coef['nwifeinc']
b1probit
Mroz87$xb <- predict(probit,type='link')
mea1 <- b1probit*dnorm(mean(Mroz87$xb))
mea1
ame1 <- b1probit*mean(dnorm(Mroz87$xb))
ame1

lm(model,data=Mroz87)

# 16.7 복잡한 경우의 한계효과

# 예제 16.6 나이 증가의 효과

Mroz1 <- Mroz87
Mroz1$age <- Mroz1$age+1
Mroz87$phat <- predict(logit,type='response')
Mroz1$phat <- predict(logit,newdata=Mroz1,type='response')
me <- Mroz1$phat-Mroz87$phat
ame <- mean(me)
ame
logit$coef['age']*mean(dlogis(predict(logit)))

# 예제 16.7 범주형 설명변수의 효과

Mroz87$haskids5 <- as.numeric(Mroz87$kids5>0)
model <- lfp~nwifeinc+educ+exper+I(exper^2)+age+haskids5+kids618
logit <- glm(model,data=Mroz87,family=binomial(link="logit"))
coeftest(logit)

Mroz0 <- Mroz87
Mroz0$haskids5 <- 0
Mroz0$phat <- predict(logit,Mroz0,type='response')
mean(Mroz0$phat)
Mroz1 <- Mroz87
Mroz1$haskids5 <- 1
Mroz1$phat <- predict(logit,Mroz1,type='response')
mean(Mroz1$phat)

# 예제 16.8 제곱항이 있을 때의 평균 한계효과

model <- lfp~nwifeinc+educ+exper+I(exper^2)+age+kids5+kids618
logit <- glm(model,data=Mroz87,family=binomial(link="logit"))
h <- 0.01
Mroz2 <- Mroz1 <- Mroz87
Mroz1$exper <- Mroz1$exper-h
Mroz2$exper <- Mroz2$exper+h
Mroz1$phat <- predict(logit,Mroz1,type='response')
Mroz2$phat <- predict(logit,Mroz2,type='response')
mean((Mroz2$phat-Mroz1$phat)/(2*h))

Mroz87$expersq <- Mroz87$exper^2
model <- lfp~nwifeinc+educ+exper+expersq+age+kids5+kids618
logit <- glm(model,data=Mroz87,family=binomial(link="logit"))
h <- 0.01
Mroz2 <- Mroz1 <- Mroz87
Mroz1$exper <- Mroz1$exper-h
Mroz2$exper <- Mroz2$exper+h
Mroz1$phat <- predict(logit,Mroz1,type='response')
Mroz2$phat <- predict(logit,Mroz2,type='response')
mean((Mroz2$phat-Mroz1$phat)/(2*h))