Pieces of R code for conference #10

######################
### conference #10 ###
######################

library("arm")

############
### 10.2 ###
############

Bypass <- read.table("bypass.data.txt", header=T)
str(Bypass)

### b
logit1 <- glm(new ~ I(age<80)
, family=binomial(link="logit"), data=Bypass)

display(logit1)

logit2 <- glm(new ~ I(age<80) + severity
, family=binomial(link="logit"), data=Bypass)

display(logit2, digits=3)

logit3 <- glm(new ~ I(age<80) + severity
, family=binomial(link="logit"), data=Bypass, subset=(age>=70 & age <=90))

display(logit3, digits=3)

logit4 <- glm(new ~ I(age<80) + severity
, family=binomial(link="logit"), data=Bypass, subset=(age>=75 & age <=85))

display(logit4, digits=3)

invlogit(coef(logit4)[1] + coef(logit4)[2] + coef(logit4)[3]*mean(Bypass$severity[Bypass$age>=75 & Bypass$age<80]))
invlogit(coef(logit4)[1] + coef(logit4)[3]*mean(Bypass$severity[Bypass$age>=80 & Bypass$age<=85]))

invlogit(coef(logit4)[1] + coef(logit4)[2] + coef(logit4)[3]*mean(Bypass$severity[Bypass$age>=75 & Bypass$age<80]))*(1-invlogit(coef(logit4)[1]>=75 & Bypass$age<80])))

invlogit(coef(logit4)[1] + coef(logit4)[3]*mean(Bypass$severity[Bypass$age>=80 & Bypass$age<=85]))*(1-invlogit(coef(logit4)[1] + coef(logit4)[3]*mean(Bypass$severity[Bypass$age>=80 & Bypass$age<=85])))

lm1 <- lm(stay ~ severity, data=Bypass)
display(lm1)

lm2 <- lm(stay ~ severity + age, data=Bypass)
display(lm2)

lm3 <- lm(stay ~ severity + age + I(age^2), data=Bypass)
display(lm3)

plot(stay ~ severity, data=Bypass, subset=new==0, ylim=c(15,45))
points(stay ~ severity, pch=3, data=Bypass, subset=new==1)

plot(stay ~ age, data=Bypass, subset=new==0, ylim=c(15,45), xlim=c(50,90))
points(stay ~ age, pch=3, data=Bypass, subset=new==1)

plot(severity ~ age, data=Bypass, subset=(new==1), xlim=range(Bypass$age), ylim=range(Bypass$severity))
points(severity ~ age, data=Bypass, subset=(new==0), pch=3)


### BONUS c ###

# without severity
Bypass$age80 <- Bypass$age - 80

TE1 <- lm(stay ~ new + age80 + age80:new, data=Bypass)
display(TE1)


# using the age threshold as a local IV
library("sem")

RDD.tsls <- tsls(stay ~ new + age80 , ~ I(age<80) + age80:I(age<80) + age80 , data=Bypass)
summary(RDD.tsls)

RDD.tsls2 <- tsls(stay ~ new + age80 + age80:new, ~ I(age<80) + age80:I(age<80) + age80 , data=Bypass)
summary(RDD.tsls2)

RDD.tsls3 <- tsls(stay ~ new + age80 , ~ I(age<80) + age80 , data=Bypass)
summary(RDD.tsls3a)


# example of a corresponding Wald estimate
RDD.itt <- lm(stay ~ I(age<80) + age80 + age80:I(age<80), data=Bypass)
display(RDD.itt)

RDD.comply <- lm(new ~ I(age<80) + age80 + age80:I(age<80), data=Bypass)
display(RDD.comply)

RDD.Wald <- coef(RDD.itt)[2]/coef(RDD.comply)[2]
RDD.Wald


# using severity
TE <- lm(stay ~ new + severity + age, data=Bypass)
display(TE)


############
### 10.3 ###
############

### a

n <- 1000
ni <- 500

h0 <- rnorm(n, mean=69.1, sd=2.9)
h1 <- h0 + rnorm(n, mean=2, sd=.5)

milk0 <- rbinom(n, size=1, prob=.3)
milk1 <- ifelse(milk0==1, 1, rbinom(n, size=1, prob=.5))

never <- milk0==0 & milk1==0
always <- milk0==1 & milk1==1
complier <- milk0==0 & milk1==1

sum(never)
sum(always)
sum(complier)

voucher <- rep(c(1,0), c(ni, n-ni))

milk <- ifelse(voucher, milk1, milk0)
h <- ifelse(milk, h1, h0)


ITT <- lm(h ~ voucher)
display(ITT)

comply <- lm(milk ~ voucher)
display(comply)

LATE.Wald <- coef(ITT)[2]/coef(comply)[2]
LATE.Wald

library("sem")

LATE.tsls <- tsls(h ~ milk, ~ voucher )
summary(LATE.tsls)


### b
n.b <- 1000

never.b <- rep(c(1, 0, 0), c(100, 100, 800))
always.b <- rep(c(0, 1, 0), c(100, 100, 800))
complier.b <- rep(c(0, 0, 1), c(100, 100, 800))

h0.b <- rnorm(n.b, mean=69.1, sd=2.9)
h1.b <- ifelse(complier.b, h0.b + rnorm(n.b, mean=2, sd=.5), h0.b + rnorm(n.b, mean=1, sd=.5))

voucher.b <- rbinom(n.b, 1, .5)

milk.b <- ifelse(never.b, 0, ifelse(always.b, 1, ifelse(voucher.b, 1, 0)))
h.b <- ifelse(milk.b, h1.b, h0.b)


ITT.b <- lm(h.b ~ voucher.b)
display(ITT.b)

comply.b <- lm(milk.b ~ voucher.b)
display(comply.b)

LATE.Wald.b <- coef(ITT.b)[2]/coef(comply.b)[2]
LATE.Wald.b

LATE.tsls.b <- tsls(h.b ~ milk.b, ~ voucher.b )
summary(LATE.tsls.b)




### c
n.c <- 1000

never.c <- rep(c(1, 0, 0), c(400, 400, 200))
always.c <- rep(c(0, 1, 0), c(400, 400, 200))
complier.c <- rep(c(0, 0, 1), c(400, 400, 200))

h0.c <- rnorm(n.c, mean=69.1, sd=2.9)
h1.c <- ifelse(complier.c, h0.c + rnorm(n.c, mean=2, sd=.5), h0.c + rnorm(n.c, mean=1, sd=.5))

voucher.c <- rbinom(n.c, 1, .5)

milk.c <- ifelse(never.c, 0, ifelse(always.c, 1, ifelse(voucher.c, 1, 0)))
h.c <- ifelse(milk.c, h1.c, h0.c)


ITT.c <- lm(h.c ~ voucher.c)
display(ITT.c)

comply.c <- lm(milk.c ~ voucher.c)
display(comply.c)

LATE.Wald.c <- coef(ITT.c)[2]/coef(comply.c)[2]
LATE.Wald.c

LATE.tsls.c <- tsls(h.c ~ milk.c, ~ voucher.c )
summary(LATE.tsls.c)

First 5 problems of practice final exam for Sciences Po course, for mardi, 15 dec, 2009

You may use the textbook and one feuille of notes (both sides). No other references. You can use a pocket calculator but no computer.

Also please note that all these examples are fake; they are not based on real data.

Below are the first five problems; I'll send more soon.

Problem 1.

1000 people are surveyed and asked how many political activities they have done in the previous year. 300 people say “none,” 300 say “one,” 200 say “two,” 100 say “three,” and 100 say “4 or more.”

Give an estimate and standard error of the average number of political activities done by people in the population. (Hint: to compute the standard error of the mean, you need to estimate the standard distribution of the responses. If you draw a histogram of the responses, that should help.)

Problem 2.

A survey is done asking people if they support government-supplied child care. 60% of the men in the survey answer yes (with a standard error of 5%), and 70% of the women answer yes (with a standard error of 5%).

Give an estimate and standard error for the difference in support between men and women.

Assuming the data were collected using a simple random sample, how many people responded to the survey?

How large would the sample size have to be for the difference in support between men and women to have a standard error of 1%?

Problem 3.

A survey is done asking people if they support affirmative action for ethnic minorities. 30% of whites say yes (with a standard error of 3%), 80% of blacks say yes (with a standard error of 5%), and70% of others say yes (with a standard error of 10%).

The population of interest is 70% white, 20% black, and 10% other. Assuming the survey is a simple random sample, give an estimate and standard error of the percentage of people in the population who support affirmative action for ethnic minorities.

Problem 4.

1100 people are surveyed: 500 say they voted for John McCain and 600 say they voted for Barack Obama. Of the McCain voters, 40% support the health-care bill that is being debated in Congress. Of the Obama voters, 70% support the bill. (People who did not vote in the past election are excluded from this study.)

Suppose you take the data (y=1 for supporting the bill and 0 for opposing it) and run a linear regression, predicting y from a variable called “obama,” which is 1 for people who say they voted for Obama and 0 for people who say they voted for McCain.

The regression output looks like this:

lm(formula = y ~ obama)
coef.est coef.se
(Intercept) 0.__ 0.__
obama 0.__ 0.__
---
n = ____, k = _
residual sd = 0.__, R-Squared = 0.__

Fill in the blanks.

Problem 5.

Students in a laboratory experiment are given the choice of which of two video games to play, a storybook adventure or a shoot-em-up game. A logistic regression is fit predicting the students’ choice (y=0 if they choose the adventure or 1 if they choose the shooting game), given their age in years. (The students’ ages in the experiment range from 18 to 25 years.)

The logistic regression output looks like this:

glm(formula = y ~ age, family = binomial(link = "logit"))
coef.est coef.se
(Intercept) 5.64 1.61
age -0.30 0.08
---
n = 230, k = 2
residual deviance = 263.5, null deviance = 281.0 (difference = 17.4)

From this model, what is the probability of choosing the shooting game if you are 20 years old?

Graph the fitted regression line (from x ranging from 18 to 25) , and put in a bunch of dots representing values of x and y that are consistent with the fitted line.

Research article to be discussed mardi 15 dec

Regression models for decision making: a cost-benefit analysis of incentives in telephone surveys. Journal of Business and Economic Statistics (2003). (Andrew Gelman, Matt Stevens, and Valerie Chan)

some R code for conference #10

######################
### conference #10 ###
######################

library("arm")

###########
### 9.8 ###
###########


n.C <- 100
n.T <- 100

x.C <- rnorm(n.C, mean=10, sd=4)
x.T <- rnorm(n.T, mean=10, sd=4)

e.C <- rnorm(n.C, mean=0, sd=1)
e.T <- rnorm(n.C, mean=0, sd=1)

y.C <- 1 + 2*x.C + e.C

y.Ta <- 1 + 2*x.T + e.T
y.Tb <- 1 + 2*x.T + e.T + 10
y.Tc <- 1 + 2*x.T + e.T + .5*x.T

plot(x.C, y.C, pch=1)
abline(lm(y.C ~ x.C))

points(x.T, y.Ta, pch=3)
abline(lm(y.Ta ~ x.T))


plot(x.C, y.C, pch=1)
abline(lm(y.C ~ x.C))

points(x.T, y.Tb, pch=3)
abline(lm(y.Tb ~ x.T))

plot(x.C, y.C, pch=1)
abline(lm(y.C ~ x.C))

points(x.T, y.Tc, pch=3)
abline(lm(y.Tc ~ x.T))


###########
### 9.9 ###
###########

### a

n.C <- 100
n.T <- 100

x.C <- rnorm(n.C, mean=10, sd=4)
x.T1 <- rnorm(n.T, mean=20, sd=4)

e.C <- rnorm(n.C, mean=0, sd=1)
e.T <- rnorm(n.T, mean=0, sd=1)

y.C <- 1 + 2*x.C + e.C
y.T1 <- 1 + 2*x.T1 + e.T + 10

y.a <- c(y.C, y.T1)
x.a <- c(x.C, x.T1)
T.a <- rep(c(0,1), c(n.C, n.T))


fit1 <- lm(y.a ~ T.a)
display(fit1)

fit2 <- lm(y.a ~ T.a + x.a)
display(fit2)

plot(x.C, y.C, pch=1, xlim=c(0,40), ylim=c(0,70))
abline(lm(y.C ~ x.C))

points(x.T1, y.T1, pch=3)
abline(lm(y.T1 ~ x.T1))

abline(h=mean(y.C))
abline(h=mean(y.T1))


### b

y.C2 <- 1 + 2*x.C + .2*I(x.C^2) + e.C
y.T2 <- 1 + 2*x.T1 + .2*I(x.T1^2) + e.T + 10

y.b <- c(y.C2, y.T2)
x.b <- c(x.C, x.T1)
T.b <- rep(c(0,1), c(n.C, n.T))


plot(x.C, y.C2, pch=1, xlim=c(0,40), ylim=c(0,200))
points(x.T1, y.T2, pch=3)

abline(lm(y.C2 ~ x.C))
abline(lm(y.T2 ~ x.T1))

fit3 <- lm(y.b ~ T.b + x.b)
display(fit3)

fit3b <- lm(y.b ~ T.b + x.b + T.b*x.b)
display(fit3b)

fit4 <- lm(y.b ~ T.b + x.b + I(x.b^2))
display(fit4)

Homework assignment due mardi, 8 dec

1. Exercise 10.2: Regression discontinuity analysis: suppose you are trying to evaluate the effect of a new procedure for coronary bypass surgery that is supposed to help with the post-operative healing process. The new procedure is risky, however, and is rarely performed in patients who are over 80 years old. Data from this (hypothetical) example are displayed in Figure 10.10.

(a) Does this seem like an appropriate setting in which to implement a regression discontinuity analysis?

(b) The folder bypass contains data for this example: stay is the length of hospital stay after surgery, age is the age of the patient, and new is the indicator variable indicating that the new surgical procedure was used. Pre-operative disease severity (severity) was unobserved by the researchers but we have access to it for illustrative purposes. Can you find any evidence using these data that the regression discontinuity design is inappropriate?

2. Exercise 10.3: Instrumental variables: come up with a hypothetical example in which it would be appropriate to estimate treatment effects using an instrumental variables strategy. For simplicity, stick to an example with a binary instrument and binary treatment variable.

(a) Simulate data for this imaginary example assuming all the assumptions are met. Estimate the local average treatment effect (LATE) for the data by dividing the intent-to-treat effect by the percentage of compliers. Show that two-stage least squares yields the same point estimate.

(b) Now simulate data in which the exclusion restriction is not met (so, for instance, those whose treatment level is left unaffected by the instrument have a treatment effect of half the magnitude of the compliers) but the instrument is strong (say, 80% of the population are compliers) and see how far off your estimate is.

(c) Finally simulate data in which the exclusion restriction is violated in the same way, however, the instrument is weak (only 20% of the population are compliers) and see how far off your estimate is.

Research article to be discussed mardi 8 dec 2009

Can losing lead to winning?, by Jonah Berger and Devin Pope.

Also see discussion here.

Homework assignment due mardi, 1 dec

1. Exercise 9.7: Gain-score models: in the discussion of gain-score models in Section 9.3, we noted that if we include the pre-treatment measure of the outcome in a gain score model, the coefficient on the treatment indicator will be the same as if we had just run a standard regression of the outcome on the treatment indicator and the pre-treatment measure. Show why this is true.

2. Exercise 9.8: Assume that linear regression is appropriate for the regression of an outcome, y, on treatment indicator, T , and a single confounding covariate, x. Sketch hypothetical data (plotting y versus x, with treated and control units indicated by circles and dots, respectively) and regression lines (for treatment and control group) that represent each of the following situations:

(a) No treatment effect,

(b) Constant treatment effect,

(c) Treatment effect increasing with x.

3. Exercise 9.9: Consider a study with an outcome, y, a treatment indicator, T , and a single confounding covariate, x. Draw a scatterplot of treatment and control observations that demonstrates each of the following:

(a) A scenario where the difference in means estimate would not capture the true treatment effect but a regression of y on x and T would yield the correct estimate.

(b) A scenario where a linear regression would yield the wrong estimate but a nonlinear regression would yield the correct estimate.

some R code for conference #8

#####################
### Conference #8 ###
#####################

library("arm")

###########
### 9.4 ###
###########

categoryZ <- 1:8
nbpersZ <- c(300, 300, 500, 500, 200, 200, 200, 200)
xZ <- rep(c(0,1),4)
xZ
TZ <- rep(rep(c(0,1),c(2,2)),2)
T
y0Z <- rep(c(4,10), c(4,4))
y1Z<- rep(c(6,12), c(4,4))

base <- data.frame(category = rep(categoryZ, nbpersZ),
x = rep(xZ, nbpersZ),
T = rep(TZ, nbpersZ),
y0 = rep(y0Z, nbpersZ),
y1 = rep(y1Z, nbpersZ))
str(base)
base$y <- ifelse(base$T == 0, base$y0, base$y1)
base$z <- ifelse(base$y0 == 4, 0, 1)

### a

# ATE = 2

# estimations

fit1 <- lm(y ~ T, data=base)
display(fit1)

mean(base$y[base$T==1]) - mean(base$y[base$T==0])

fit2 <- lm(y ~ T + x, data=base)
display(fit2)


fit3 <- lm(y ~ T + x + T:x, data=base)
display(fit3)

m0 <- mean(base$y[base$T==1 & base$x==0]) - mean(base$y[base$T==0 & base$x==0])
m1 <- mean(base$y[base$T==1 & base$x==1]) - mean(base$y[base$T==0 & base$x==1])
m0
m1


####
fit4 <- lm(y ~ T + z + T:z, data=base)
display(fit4)

R codes for book

Go to my webpage, then to the link for my book with Jennifer, then click on Data and software, then Directories with data, then you can download the zip files with data, code for the book, and code for the book.

Recently in class we've been discussing the examples from chapters 7 and 8.

Homework assignment due mardi, 24 nov

1. Exercise 9.1: Suppose you are interested in the effect of the presence of vending machines in schools on childhood obesity. What randomized experiment would you want to
do (in a perfect world) to evaluate this question?

2. Exercise 9.4 (b,c,d): See statement of Exercise 9.4 from last week's homework assignment. Then:

(b) Is it plausible to believe that these data came from a randomized experiment? Defend your answer.

(c) Another population quantity is the mean of y for those who received the treatment minus the mean of y for those who did not. What is the relation between this quantity and the average treatment effect?

(d) For these data, is it plausible to believe that treatment assignment is ignorable given sex? Defend your answer.

3. Exercise 9.5: For the hypothetical study in the previous exercise, figure out the estimate and the standard error of the coefficient of T in a regression of y on T and x.

4. Exercise 9.6: You are consulting for a researcher who has performed a randomized trial where the treatment was a series of 26 weekly therapy sessions, the control was no therapy, and the outcome was self-report of emotional state one year later. However, most people in the treatment group did not attend every therapy session. In fact there was a good deal of variation in the number of therapy sessions actually attended. The researcher is concerned that her results represent “watered down” estimates because of this variation and suggests adding in another predictor to the model: number of therapy sessions attended. What would you advise her?

some R code for conference #7

#####################
### conference #7 ###
#####################

library("arm")
library("car")

###########
### 2.3 ###
###########

n.sum <- 20
n.sim <- 1000

simCLT <- replicate(n.sim, sum(runif(n.sum)))
hist(simCLT, probability=TRUE, xlim=c(n.sum/2-3*sqrt(n.sum/12), n.sum/2+3*sqrt(n.sum/12)))
curve(dnorm(x, mean=n.sum/2, sd=sqrt(n.sum/12)), add=T)

# use number of bins according to car:
hist(simCLT, probability=TRUE, xlim=c(n.sum/2-3*sqrt(n.sum/12), n.sum/2+3*sqrt(n.sum/12)), nclass=n.bins(simCLT))
curve(dnorm(x, mean=n.sum/2, sd=sqrt(n.sum/12)), add=T)


###########
### 8.1 ###
###########

### y = 3 + 0.1*x1 + 0.5*x2 + 5*error, error ~ T(mean=0, scale=5, df=4)

### a

n.obs <- 100

x1 <- 1:n.obs
x2 <- rbinom(n.obs, size=1, prob=0.5)
error <- rt(n.obs, df=4)
true.coef <- c(3, 0.1, 0.5, 5)
y <- cbind(rep(1, n.obs), x1, x2, error) %*% true.coef
y[1:10]
plot(x1, y)

fit1 <- lm(y ~ x1 + x2)
display(fit1)
hat.coef <- coef(fit1)
hat.coef.se <- se.coef(fit1)
hat.coef.se <- sqrt(diag(vcov(fit1)))
cover.68 <- abs(hat.coef - true.coef[1:3]) < hat.coef.se
cover.68

### b

sim.T <- function(c.int, c1, c2, scale, nobs){
x1 <- 1:nobs
x2 <- rbinom(nobs, size=1, prob=0.5)
error <- rt(nobs, df=4)
true.coef <- c(c.int, c1, c2, scale)
y <- cbind(rep(1, nobs), x1, x2, error) %*% true.coef
fit1 <- lm(y ~ x1 + x2)
hat.coef <- coef(fit1)
hat.coef.se <- sqrt(diag(vcov(fit1)))
abs(hat.coef - true.coef[1:3]) < hat.coef.se
}

res.sim <- replicate(1000, sim.T(3, 0.1, 0.5, 5, 100))

#dim(res.sim)
#is.array(res.sim)

apply(res.sim, 1, mean)


### c
library("hett")

fit2 <- tlm(y ~ x1 + x2, start=list(dof=4), estDof=FALSE)
summary(fit2)
fit2
hat.coef <- coef(fit2$loc.fit)
hat.coef.se <- sqrt(diag(summary.tlm(fit2)$loc.summary[["cov.unscaled"]]))
abs(hat.coef - true.coef[1:3]) < hat.coef.se

sim.T2 <- function(c.int, c1, c2, scale, nobs){
x1 <- 1:nobs
x2 <- rbinom(nobs, size=1, prob=0.5)
error <- rt(nobs, df=4)
true.coef <- c(c.int, c1, c2, scale)
y <- cbind(rep(1, nobs), x1, x2, error) %*% true.coef
fit2 <- tlm(y ~ x1 + x2, start=list(dof=4), estDof=FALSE)
hat.coef <- coef(fit2$loc.fit)
hat.coef.se <- sqrt(diag(summary.tlm(fit2)$loc.summary[["cov.unscaled"]]))
abs(hat.coef - true.coef[1:3]) < hat.coef.se
}

res.sim2 <- replicate(1000, sim.T2(3, 0.1, 0.5, 5, 100))

#dim(res.sim)
#is.array(res.sim)

apply(res.sim2, 1, mean)


### compare standard error of the parameter estimates

sim.T.se <- function(c.int, c1, c2, scale, nobs){
x1 <- 1:nobs
x2 <- rbinom(nobs, size=1, prob=0.5)
error <- rt(nobs, df=4)
true.coef <- c(c.int, c1, c2, scale)
y <- cbind(rep(1, nobs), x1, x2, error) %*% true.coef
fit1 <- lm(y ~ x1 + x2)
se.coef(fit1)
}
sim.T2.se <- function(c.int, c1, c2, scale, nobs){
x1 <- 1:nobs
x2 <- rbinom(nobs, size=1, prob=0.5)
error <- rt(nobs, df=4)
true.coef <- c(c.int, c1, c2, scale)
y <- cbind(rep(1, nobs), x1, x2, error) %*% true.coef
fit2 <- tlm(y ~ x1 + x2, start=list(dof=4), estDof=FALSE)
sqrt(diag(summary.tlm(fit2)$loc.summary[["cov.unscaled"]]))
}

sim.se1 <- replicate(1000, sim.T.se(3, 0.1, 0.5, 5, 100))
sim.se2 <- replicate(1000, sim.T2.se(3, 0.1, 0.5, 5, 100))
apply(sim.se1, 1, mean)
apply(sim.se2, 1, mean)

Homework assignment due mardi, 17 nov

1. Exercise 9.3: Suppose you are a consultant for a researcher who is interested in investigating the effects of teacher quality on student test scores. Use the strategy of mapping this question to a randomized experiment to help define the question more clearly. Write a memo to the researcher asking for needed clarifications to this study proposal.

2. Exercise 9.4 (a): The table below describes a hypothetical experiment on 2400 persons. Each row of the table specifies a category of person, as defined by his or her pre-treatment predictor x, treatment indicator T, and potential outcomes y0, y1. (For simplicity, we assume unrealistically that all the people in this experiment fit into these eight categories.)

Category #persons.in.category x T y0 y1
1 300 0 0 4 6
2 300 1 0 4 6
3 500 0 1 4 6
4 500 1 1 4 6
5 200 0 0 10 12
6 200 1 0 10 12
7 200 0 1 10 12
8 200 1 1 10 12

In making the table we are assuming omniscience, so that we know both y0 and y1 for all observations. But the (nonomniscient) investigator would only observe x, T, and yT for each unit. (For example, a person in category 1 would have x=0, T=0, y=4, and a person in category 3 would have x=0, T=1, y=6.)

(a) What is the average treatment effect in this population of 2400 persons?

Voter turnout in the United States

As requested, here is a link to the research on the political attitudes of nonvoters and the potential effects of increased voter turnout.

Research article to be discussed mardi, 17 nov

Follow this link and read the linked article in Behavioral and Brain Sciences.

Also this.

Research article to be discussed mardi 10 nov

The Accident Externality from Driving, by Aaron Edlin and Pinar Mandic.

We [Edlin and Karaca-Mandic] estimate auto accident externalities (more specifically insurance externalities) using panel data on state-average insurance premiums and loss costs. Externalities appear to be substantial in traffic-dense states: in California, for example, we find that the increase in traffic density from a typical additional driver increases total statewide insurance costs of other drivers by $1,725-$3,239 per year, depending on the model. High-traffic density states have large economically and statistically significant externalities in all specifications we check. In contrast, the accident externality per driver in low-traffic states appears quite small. On balance, accident externalities are so large that a correcting Pigouvian tax could raise $66 billion annually in California alone, more than all existing California state taxes during our study period, and over $220 billion per year nationally.

Homework assignment, due mardi, 10 nov

1. Exercise 2.3: Demonstration of the Central Limit Theorem: let x = x1 + · · · + x20, the sum of 20 independent Uniform(0,1) random variables. In R, create 1000 simulations of x and plot their histogram. On the histogram, overlay a graph of the normal density function. Comment on any differences between the histogram and the curve.

2. Exercise 8.1: Fitting the wrong model: suppose you have 100 data points that arose from the following model: y = 3+0.1x1+0.5x2+error, with errors having a t distribution with mean 0, scale 5, and 4 degrees of freedom. We shall explore the implications of fitting a standard linear regression to these data.

(a) Simulate data from this model. For simplicity, suppose the values of x1 are simply the integers from 1 to 100, and that the values of x2 are random and equally likely to be 0 or 1. (In R, you can define x.1 <- 1:100, simulate x.2 using rbinom(), then create the linear predictor, and finally simulate the random errors in y using the rt() function.) Fit a linear regression (with normal errors) to these data and see if the 68% confidence intervals for the regression coefficients (for each, the estimates ±1 standard error) cover the true values.

(b) Put the above step in a loop and repeat 1000 times. Calculate the confidence
coverage for the 68% intervals for each of the three coefficients in the model.

(c) Repeat this simulation, but instead fit the model using t errors (see Exercise 6.6).

Some R code for conference #6

#####################
### conference #6 ###
#####################

library("arm")

###########
### 7.1 ###
###########

### a

f.shots <- function(p.success=.6, hits.stop=2){
shots <- 0
missed <- 0
success <- 0
while (missed < hits.stop){
if (rbinom(1,1,p.success)==1) {
shots <- shots+1
success <- success + 1
missed <- 0
} else {
shots <- shots+1
missed <- missed + 1
}
}
c(shots=shots,success=success)
}
f.shots()

### b
sim.shots <- rep(NA,1000)
mytime <- proc.time()
for (i in 1:1000){
sim.shots[i] <- f.shots()["shots"]
}
proc.time() - mytime

mean(sim.shots)
sd(sim.shots)

# slightly faster with "replicate"
mytime <- proc.time()
sim.shots2 <- replicate(1000, f.shots()["shots"])
proc.time() - mytime

mean(sim.shots2)
sd(sim.shots2)


# theoretical results:

is the probability of success





sim.shots.success <- data.frame(shots=rep(NA,1000), success=rep(NA,1000))
for (i in 1:1000){
sim.shots.success[i,] <- f.shots()
}
str(sim.shots.success)
sim.shots.success[1:10,]
sim.shots.success$prop.success <- sim.shots.success$succes/sim.shots.success$shots

mean(sim.shots.success$shots)
mean(sim.shots.success$success)
sd(sim.shots.success$shots)
sd(sim.shots.success$success)

### c
plot(sim.shots.success$shots, sim.shots.success$prop.success)

#############
### BONUS ###
#############

### if what we want is the total number of non missed shots
f.hits <- function(p.success=.6, hits.stop=2){
hits <- 0
missed <- 0
while (missed < hits.stop){
if (rbinom(1,1,p.success)==1) {
hits <- hits+1
missed <- 0
} else {
missed <- missed + 1
}
}
hits
}
f.hits()

sim.hits <- rep(NA,10000)
mytime <- proc.time()
for (i in 1:10000){
sim.hits[i] <- f.hits()
}
proc.time() - mytime

mean(sim.hits)
sd(sim.hits)

## mat10000 <- matrix(c(.6,2), nrow=2, ncol=10000)
## mytime <- proc.time()
## sim.hits2 <- apply(mat10000, 2, f.hits)
## proc.time() - mytime

## mean(sim.hits2)
## sd(sim.hits2)

# theoretical results in that case:
# compute the probability to take n hits by induction
# p(0)=.4*.4
# a sequence leading to n hits starts either by 1 or by 0
# the set of sequences leading to n hits and starting with a 0 is the same as the set of sequences leading to n hits and starting with a 1, but adding a 0 at the beginning of the sequence
# n -> n+1
# the set of sequences leading to n+1 hits and starting with a 1 is the same as the set of sequences leading to n hits and starting with a 1, but adding a 1 at the beginning plus the set of sequences starting with a 1 but adding 10 at the beginning






# let H be the random variable equal to the number of hits











###########
### 7.3 ###
###########

# let's assume all random variables are independent, continuous and normally distributed

f.savings <- function(m.unit=5, sd.unit=4, m.market=40000, s.market=10000){
market.size <- round(rnorm(1, mean=m.market, sd=s.market),0)
if (market.size<0) {0 <- market.size}
savings <- rnorm(market.size, mean=m.unit, sd=sd.unit)
sum(savings)
}
f.savings()

total.savings <- rep(NA,1000)
for (i in 1:1000){
total.savings[i] <- f.savings()
}
mean(total.savings)
sd(total.savings)

# the order of magnitude should be about s.market*m.unit=10000*5

# another story:
# let's assume that there are 40000 firms
# each firm will independently enter the market with probability 1/2501
# each firm which enters the market will buy 2501 widgets
# otherwise it will buy 0 widgets

f.savings2 <- function(m.unit=5, sd.unit=4){
market.size <- 2501*rbinom(1,40000,1/2501)
savings <- rnorm(market.size, mean=m.unit, sd=sd.unit)
sum(savings)
}
f.savings2()
total.savings2 <- rep(NA,1000)
for (i in 1:1000){
total.savings2[i] <- f.savings2()
}
mean(total.savings2)
sd(total.savings2)

###########
### 7.8 ###
###########

unidim.sim <- function(estimate, sd.estimate, df, n){
X <- rchisq(n, df=df)
rnorm(n, mean=estimate, sd=sd.estimate*sqrt(df/X))
}
unidim.sim(estimate=600, sd.estimate=400, df=50, n=5)

plot(unidim.sim(estimate=600, sd.estimate=400, df=50, n=1000), unidim.sim(estimate=3, sd.estimate=1, df=100, n=1000))

cost.eff <- rep(NA,1000)
for (i in 1:1000){
cost.eff[i] <- mean(unidim.sim(estimate=600, sd.estimate=400, df=50, n=1000))/mean(unidim.sim(estimate=3, sd.estimate=1, df=100, n=1000))
}

hist(cost.eff)
quantile(cost.eff,probs=c(.025, .25, .75, .975))

cost.eff2 <- rep(NA,1000)
for (i in 1:1000){
cost.eff2[i] <- mean(unidim.sim(estimate=600, sd.estimate=400, df=50, n=1000))/mean(unidim.sim(estimate=3, sd.estimate=2, df=100, n=1000))
}

hist(cost.eff2)
quantile(cost.eff2,probs=c(.025, .25, .75, .975))

Hints for conference #6

# You can complete the following pieces of code

###########
### 7.1 ###
###########

### a

f.shots <- function(p.success=.6, hits.stop=2){
shots <- 0
missed <- 0
success <- 0
while ( ### condition to complete ### ){

### actions to complete here ###

}

# result:
c(shots=shots,success=success)
}
f.shots()

### b
# initialize a vector of size 1000 with missing values
sim.shots <- rep(NA,1000)
# fill in the missing values...

mean(sim.shots)
sd(sim.shots)

# do the same with both shots and successes:
sim.shots.success <- data.frame(shots=rep(NA,1000), success=rep(NA,1000))


mean(sim.shots.success$shots)
mean(sim.shots.success$success)
sd(sim.shots.success$shots)
sd(sim.shots.success$success)

### c
# plot the result


###########
### 7.3 ###
###########

# let's assume all random variables are independent, continuous and normally distributed

f.savings <- function(m.unit=5, sd.unit=4, m.market=40000, s.market=10000){
# calculate market size:
market.size <-

# calculate vector of savings and sum it

f.savings()

# initialize
total.savings <- rep(NA,1000)
# fill in

mean(total.savings)
sd(total.savings)

# another story:
# let's assume that there are 40000 firms
# each firm will independently enter the market with probability 1/2501
# each firm which enters the market will buy 2501 widgets
# otherwise it will buy 0 widgets

# do the same

###########
### 7.8 ###
###########

# create function according to the procedure p.143
# use rchisq and rnorm
unidim.sim <- function(estimate, sd.estimate, df, n){
}

# try it
unidim.sim(estimate=600, sd.estimate=400, df=50, n=5)

# plot the result with 1000 points

# 1000 simulations:
cost.eff <- rep(NA,1000)
# fill in

# histogram
hist(cost.eff)

# use function quantile for the confidence

# try again with different values

some R code for conference #5

### conf #5

library("arm")

###########
### 5.5 ###
###########

### a

# plot of invlogit
curve(invlogit(x), from=-10, to=10)

# plot of invlogit(-24 +.4*x) with different ranges of x
curve(invlogit(-24 +.4*x), from=0, to=100)
curve(invlogit(-24 +.4*x), from=45, to=75)
curve(invlogit(-24 +.4*x), from=30, to=90)

midterm <- rnorm(50, mean=60, sd=15)
midterm <- ifelse(midterm<0, 0, ifelse(midterm>100, 100, midterm))

jitter.binary <- function(a, jitt=.05){
ifelse(a==0, runif(length(a), 0, jitt), runif(length(a), 1-jitt, 1))
}

# generate consistent final results using the model and midterm results
# do not just use the predicted value...
final.bad <- ifelse(invlogit(-24 + .4*midterm)>.5, 1, 0)
final.good <- rbinom(length(midterm), 1, invlogit(-24 + .4*midterm))

# little test to check the syntax above
# test <- c(rep(.1,100),rep(.9,100))
# rbinom(length(test), 1, test)

plot(midterm,jitter.binary(final.good), cex=.25)
curve(invlogit(-24 + .4*x), add=T)
abline(v=60)

plot(midterm,jitter.binary(final.bad), cex=.25)
curve(invlogit(-24 + .4*x), add=T)
abline(v=60)

### b








# test the transformation
logit1 <- glm(final.good ~ midterm, family=binomial(link="logit"))
display(logit1)

midterm.cs <- (midterm-mean(midterm))/sd(midterm)
logit2 <- glm(final.good ~ midterm.cs, family=binomial(link="logit"))
display(logit2)

### c
newpred <- rnorm(length(midterm), 0, 1)
logit3 <- glm(final.good ~ midterm + newpred, family=binomial(link="logit"))
display(logit3)

###########
### 5.6 ###
###########

curve(dlogis(x), from=-6, to=6)

curve(dlogis(x, location=4.5), from=-2, to=11)
abline(h=0)
#abline(v=0)
polygon(c(0,seq(0,12, length=100),12), c(0,dlogis(seq(0,12, length=100), location=4.5),0), col="grey")
polygon(c(0,seq(0,12, length=10),12), c(0,dlogis(seq(0,12, length=10), location=4.5),0), col="grey")

Homework assignment, due mardi, 27 oct

1. Exercise 7.1: Discrete probability simulation: suppose that a basketball player has a 60% chance of making a shot, and he keeps taking shots until he misses two in a row. Also assume his shots are independent (so that each shot has 60% probability of success, no matter what happened before).

(a) Write an R function to simulate this process.

(b) Put the R function in a loop to simulate the process 1000 times. Use the simulation to estimate the mean, standard deviation, and distribution of the total number of shots that the player will take.

(c) Using your simulations, make a scatterplot of the number of shots the player will take and the proportion of shots that are successes.

2. Exercise 7.3: Propagation of uncertainty: we use a highly idealized setting to illustrate the use of simulations in combining uncertainties. Suppose a company changes its technology for widget production, and a study estimates the cost savings at $5 per unit, but with a standard error of $4. Furthermore, a forecast estimates the size of the market (that is, the number of widgets that will be sold) at 40,000, with a standard error of 10,000. Assuming these two sources of uncertainty are independent, use simulation to estimate the total amount of money saved by the new product (that is, savings per unit, multiplied by size of the market).

3. Exercise 7.8: Inference for the ratio of parameters: a (hypothetical) study compares the costs and effectiveness of two different medical treatments.

• In the first part of the study, the difference in costs between treatments A and B is estimated at $600 per patient, with a standard error of $400, based on a regression with 50 degrees of freedom.

• In the second part of the study, the difference in effectiveness is estimated at 3.0 (on some relevant measure), with a standard error of 1.0, based on a regression with 100 degrees of freedom.

• For simplicity, assume that the data from the two parts of the study were collected independently. Inference is desired for the incremental cost-effectiveness ratio: the difference between the average costs of the two treatments, divided by the difference between their average effectiveness. (This problem is discussed further by Heitjan, Moskowitz, and Whang, 1999.)

(a) Create 1000 simulation draws of the cost difference and the effectiveness difference, and make a scatterplot of these draws.

(b) Use simulation to come up with an estimate, 50% interval, and 95% interval for the incremental cost-effectiveness ratio.

(c) Repeat this problem, changing the standard error on the difference in effectiveness to 2.0.