a1 <- .5
a2 <- .5
b1 <- .5
b2 <- .5
b3 <- .5
x1 <- rnorm(1000)
x2 <- rnorm(1000)
x3 <- rnorm(1000)
denom <- 1 - a1*a2
y1star <- b1*x1/denom + a1*b2*x2/denom + (b3 + a1*b3)*x3/denom + rnorm(1000)
y2star <- b2*x2/denom + a2*b1*x1/denom + (b3 + a2*b3)*x3/denom + rnorm(1000)
y1 <- (y1star>0)
y2 <- (y2star>0)
junk <- data.frame(y1=y1,y2=y2,x1=x1,x2=x2,x3=x3)
sprob <- function(data, ind) {
Data <- data[ind,]
rf1 <- glm(y1 ~ x1 + x2 + x3, family=binomial(link=probit), data=Data)
rf2 <- glm(y2 ~ x1 + x2 + x3, family=binomial(link=probit), data=Data)
y1hat <- predict(rf1, type="link")
y2hat <- predict(rf2, type="link")
sf1 <- glm(y1 ~ y2hat + x1 + x3, family=binomial(link=probit), data=Data)
sf2 <- glm(y2 ~ y1hat + x2 + x3, family=binomial(link=probit), data=Data)
cbind(sf1$coefficients,sf2$coefficients)
}
res.boot <- boot(junk, sprob, 1000)