Simulate fake data to assess model fit

simconf.r

#set the working directoy and read the foreign library
setwd("/.../")
library(foreign)

#read the data and remove missing values of the dependent variable
data <- read.dta("repdata.dta")
data <- data[data$onset != 4,]

#estimate the model
model <- glm(onset ~ warl + gdpenl + lpopl1 + lmtnest + ncontig + Oil + nwstate + instab + polity2l + ethfrac + relfrac, data = data, family = "binomial")
summary(model)

#create a matrix with the variables used in the model and remove rows with missing values
X <- as.matrix(cbind(1,model$data[,c("warl","gdpenl","lpopl1","lmtnest","ncontig","Oil","nwstate","instab","polity2l","ethfrac","relfrac")]))
X <- na.omit(X)

#get the beta coefficients and variance-covariance matrix from the model
beta <- coef(model)
covvar <- vcov(model)

#create a matrix to hold the fake data sets
n <- nrow(X)
out <- matrix(nr = n, nc = 1000)

#simulate a 1000 fake data sets
for (s in 1:1000){
  b <- mvrnorm(1,beta,covvar)
  xb <- X %*% b
  p <- 1/(1 + exp(-xb))
  y.fake <- rbinom(n,1,p)
  out[,s] <- y.fake
}

#for each fake data set calculate the proportion of 1s
pred <- apply(out,2,function(x){
  p <- sum(x == 1)/length(x)
  return(p)
})

#get the proportion from the observed data
p.true <- sum(data$onset == 1)/length(data$onset)

#put the simulated proportions in a data frame for ggplot2
pred <- data.frame(pred = pred)

#plot the density of the simulated proportions and mark the observed proportion
plot <- ggplot(pred, aes(x = pred)) + 
        geom_density() + 
        theme_bw() + 
        geom_segment(aes(x = p.true, xend = p.true, y = 0, yend = 100), colour = "red", arrow=arrow(length=unit(0.3,"cm"), ends = "first"))

plot