#define pvw function, which performs predictive value weighting for the particular setup we will
#use to simulate data

pvw <- function(data, indices, sens, spec) {

localData <- data[indices,]

zyc_mod <- glm(z~y+c+y:c, family="binomial", data=localData)
localData$piycstar <- predict(zyc_mod,type=c("response"))

localData$se <- rep(sens,n)
localData$se[localData$piycstar>=localData$se] <- localData$piycstar[localData$piycstar>=localData$se]+0.001

localData$sp <- rep(spec,n)
localData$sp[localData$piycstar<=(1-localData$sp)] <- 1-localData$piycstar[localData$piycstar<=(1-localData$sp)]+0.001

localData$term1 <- (localData$se-1)*localData$piycstar*(localData$se*(localData$piycstar-1))^(-1)
localData$term2 <- (localData$sp-1)*(localData$piycstar-1)*(localData$sp*localData$piycstar)^(-1)
localData$det <- 1/(localData$term1*localData$term2-1)
localData$ppv <- localData$det*(localData$term2-1)
localData$npv <- localData$det*(localData$term1-1)

expandedData <- rbind(localData,localData)
expandedData$x <- 0
expandedData$x[(n+1):(2*n)] <- 1

expandedData$wgt <- 0
expandedData$wgt[expandedData$x==1 & expandedData$z==1] <- expandedData$ppv[expandedData$x==1 & expandedData$z==1]
expandedData$wgt[expandedData$x==0 & expandedData$z==1] <- 1-expandedData$ppv[expandedData$x==0 & expandedData$z==1]
expandedData$wgt[expandedData$x==1 & expandedData$z==0] <- 1-expandedData$npv[expandedData$x==1 & expandedData$z==0]
expandedData$wgt[expandedData$x==0 & expandedData$z==0] <- expandedData$npv[expandedData$x==0 & expandedData$z==0]

wgtMod <- glm(y~x+c, family="quasibinomial",weights=wgt,data=expandedData)
return(coef(wgtMod))
}


#perform simulation study of pvw

sims <- 1000
simEsts <- array(0, dim=c(sims,3))

for (i in 1:sims) {
print(i)
n <- 1000
x <- 1*(runif(n)<0.5)
c <- 1*(x==1)*(runif(n)<0.7) + 1*(x==0)*(runif(n)<0.3)
pr <- exp(x+c)/(1+exp(x+c))
y <- 1*(runif(n)<pr)
summary(glm(y~x+c, family="binomial"))

z <- 1*(x==1)*(runif(n)<0.9) + 1*(x==0)*(runif(n)<0.2)
summary(glm(y~z+c, family="binomial"))

originalData <- data.frame(y,z,c)

pvw(originalData,1:n,0.9,0.8)

simEsts[i,] <- pvw(originalData,1:n)

}

#finally we look to see how the mean of the estimates compare to their true values
colMeans(simEsts)