# This is a time-series dataset of Germany cities.

dat <- read_csv("germany_data.csv") %>%
  # Let's select our key variables of interest
  select(id, # Unique ID Code
         muni, # Name of German District Name 
         solar, # No. of Solar Farms built between 2008 and 2011
         pvout, # Sunlight (Average Photovoltaic Output)
         area, # Area in square kilometers
         pop, # Population in 2000
         land_price, # Cost of Land in 2000 in euros per square kilometer
         income, # Income per taxpayer in 2001 in thousands of euros
         unemployment, # Unemployment rate in 2001
         crime_rate, # Crime Rate in 2007
         # (a common proxy for bonding social capital
         # when survey measures are unavailable)
         voter_turnout, # Voter Turnout in 2005
         # (proxy for bridging social capital)
         winning_party, # % Votes for Ruling Paer (CDU-CSU) in 2005
         green) %>% # % Votes for Green Party in 2005
    # Let's zoom into cities with at least 1 solar farm,
  # Since the other cities are probably not very comparable to the rest
  filter(solar > 0) 

# Make an Negative Binomial model using our dataset
z <- dat %>%

  # Let's model the association between solar farm adoption 
  # and voter turnout in that community
  zelig(formula = solar ~ voter_turnout +
       # While controlling for other proxies for other types of social capital
                 crime_rate + winning_party +
         # environmentalism and socioeconomics
         green + income +  unemployment +
         # and technical conditions
         pvout + area + pop + land_price, model = "negbin", cite = FALSE) 

# Let's check out the effect of voter turnout
z

# Let's get the 20th, 40th, 60th, and 80th percentiles
range <- quantile(dat$voter_turnout, probs = c(0.2,0.4,0.6,0.8), na.rm = TRUE)

# And let's generate simulated effects
z_sim <- z %>%
  # We can set the specified range of values here (5 evenly spaced points between the 20th and 80th percentiles),
  # and then all other covariates will be held at their means
  setx(voter_turnout = range) %>%
  # and running the simulations
  sim()

# We can view them here.
# We get 1000 simulated values for each point in our range.
# for both predicted values (adjusts for estimation uncertainty)
# and also for expected values (adjust for both estimation and fundamental uncertainty)
# You'll notice that predicted values have wider confidence intervals, while expected values are narrower.
# You should probably used predicted values for things like election forecasting, 
# because they give you a wider margin of error,
# while expected values are great for explaining past outcomes.
z_sim

z_sim %>%
  # Let's get these in a dataset
  zelig_qi_to_df() %>%
  # and shrink them down to just the median values and 95% confidence intervals
  # Let's get them for just expected values
  qi_slimmer(qi_type = "ev", ci = 0.95) %>%
  # And we can map them here!
  ggplot(mapping = aes(x = voter_turnout, y = qi_ci_median, 
                       ymin = qi_ci_min, ymax = qi_ci_max)) +
  geom_ribbon(alpha = 0.5, fill = "firebrick") +
  geom_line() +
  theme_minimal(base_size = 14) +
  theme(plot.caption = element_text(hjust = 0.5)) +
  labs(x = "(%) Voter Turnout",
       y = "Expected Solar Farms\n(in an average city)",
       caption = "Based on 1000 simulations made with the sim() function in Zelig.")

# First, let's extract the model from its zelig format, 
# so we can work more easily with it
m <- from_zelig_model(z)

# Get observed alpha and beta coefficients from your model, which we will call "gamma"
gamma <- m$coefficients

# Get variance covariance matrix from your model, called "sigma"
sigma <- vcov(m)


# Grab some average traits from our data, a vector we will call "x"
# (This is nearly the same code that we used to make fakedata before, just from a new object)
x <- m$model %>% 
  summarize(
    # Set Voter turnout at a specific value (20th percentile)
    voter_turnout = range[1],
    # But leave everything else at their means
    crime_rate = mean(crime_rate, na.rm = TRUE),
    winning_party = mean(winning_party, na.rm = TRUE), 
    green = mean(green, na.rm = TRUE), 
    income = mean(income, na.rm = TRUE),
    unemployment = mean(unemployment, na.rm = TRUE),
    pvout = mean(pvout, na.rm = TRUE),
    area = mean(area, na.rm = TRUE),
    pop = mean(pop, na.rm = TRUE),
    land_price = mean(land_price, na.rm = TRUE)) %>%
  # Let's convert it into a matrix too.
  as.matrix()

# Let's grab our alpha coefficient, the intercept
alpha <- gamma[1] # grab just the first, which is the intercept

# Let's grab our beta coefficients
beta <- gamma[-1] # grab all but the first (which was the intercept)

# Let's get theta (sometimes written as eta for negative binomial models)
# Theta is the total of X*Beta + Alpha
theta <- sum(beta * x + alpha)

# Now, we will multiply theta by the link function to get the outcome. 
# In the case of ordinary least squares models,
# the link function is non-existent; 
# it is equal to multiplying theta by 1.
# However, many other models have link functions, 
# which transform the total product of our coefficients and data into an outcome.
# Negative binomial models have a very simple log-link function 
# (it looks like: log(y) = X * Beta + Alpha).
# So all we have to do is exponentiate it.
y_hat <- exp(theta)

# Let's take a look
y_hat

# Generate 1000 simulations,
# drawing from a multivariate normal distribution
gamma_hat <- MASS::mvrnorm(n = 1000, mu = gamma, Sigma = sigma)
# Each of these variables are related, based on the variance-covariance matrix,
# even though they were randomly drawn.

# Let's check out the first two lines
gamma_hat %>% head(2)

# We have 11 columns, one for each coefficient, and 1000 rows, one per simulation.

# Get simulated version of alpha (meaning 1000 rows of the simulated Intercept coefficient)
alpha_hat <- gamma_hat[,1] # this little [,1] says, please keep just the first column

# Get simulated version of beta coefficients (meaning 1000 rows of the 10 beta coefficients)
beta_hat <- gamma_hat[,-1] # this little [,-1] says, please remove the first column

# Let's multiply each simulated beta coefficient by its corresponding set x value,
# and then add them up as "eta_hat" 
# (This will give us 5 columns, because we had 5 rows of fakedata in x to simulate)
sim_out <- exp(beta_hat %*% t(x) + alpha_hat)

# These are our 1000 predicted values
pv <- MASS::rnegbin(n = 1000, mu = sim_out, theta = m$theta)

# Let's repeat our zelig simulation again
z_sim <- z %>%
  setx(voter_turnout = range) %>%
  sim()
# Let's bind both together
bind_rows(
  # Grabbing our self-made predicted values,
  data.frame(pv = pv,
           type = "Do-it-yourself"), 
  # And comparing them against zelig's predicted values
  # (As you can see from the code, I had to go digging into the model object to find them)
  data.frame(pv = z_sim$sim.out$range[[1]]$pv %>% unlist(),
           type = "Zelig")) %>%
  # And mapping each to different parts of a visual
  ggplot(mapping = aes(x = type, y = pv, fill = type)) +
  # Displaying the points, jittered
  geom_jitter(alpha = 0.5, color = "grey") +
  # Using violin plots to map their distributions.
  geom_violin(alpha = 0.75) +
  theme_minimal(base_size = 14) +
  labs(x = "Type of Simulation", y = "Predicted # of Solar Farms Adopted\ngiven low voter turnout (20th percentile)") +
  guides(fill = "none")

#If you re-run this code many times, you'll find they have almost identical shapes.

# Get simulated coefficients
gamma_hat <- MASS::mvrnorm(n = 1000, mu = gamma, Sigma = sigma)
alpha_hat <- gamma_hat[,1]
beta_hat <- gamma_hat[,-1]
# Now let's take this portion of the model and add the simulated intercept,
# and then exponentiate both, to deal with the log-link function in the model
sim_out <- exp(beta_hat %*% t(x) + alpha_hat)

# Let's write a quick function to produce 1 expected value averaged over 100 simulations from a negative binomial distribution,
fundamental = function(i){
  # Take this simulated outcome from a multivariate normal distribution,
  # and now make a negative binomial distribution with a mean of that value,
  # and grab 1000 random simulations.
  # This approximates the effect of fundamental uncertainty on your value
  rnbinom(n = 1000, size = 1, mu = sim_out[i]) %>%
    # Now let's take the mean of each, averaging over fundamental uncertainty.
    mean() %>%
    return()
}

# These are our 1000 expected values
ev <- 1:1000 %>%
  map(~fundamental(.)) %>%
  unlist()

# Let's repeat our zelig simulation again
z_sim <- z %>%
  setx(voter_turnout = range[1]) %>%
  sim()


# Let's bind both together
bind_rows(
  # Grabbing our self-made expected values,
  data.frame(ev = ev,
           type = "Do-it-yourself"), 
  # And comparing them against zelig's expected values
  # (As you can see from the code, I had to go digging into the model object to find them)
  data.frame(ev = z_sim$sim.out[[1]]$ev %>% unlist(),
           type = "Zelig")) %>%
  # And mapping each to different parts of a visual
  ggplot(mapping = aes(x = type, y = ev, fill = type)) +
  # Displaying the points, jittered
  geom_jitter(alpha = 0.5, color = "grey") +
  # Using violin plots to map their distributions.
  geom_violin(alpha = 0.75) +
  theme_minimal(base_size = 14) +
  labs(x = "Type of Simulation", y = "Expected # of Solar Farms Adopted\ngiven low voter turnout (20th percentile)") +
  guides(fill = "none")

#If you re-run this code many times, you'll find they have almost identical shapes.