Comparar texto

Encuentra la diferencia entre dos archivos de texto

Visualización

Vista

Editor liveOcultar espacios en blancoOcultar sin cambiosSin ajuste de línea

Nivel de detalle

Sintaxis

Procesamiento

Diffchecker Desktop La forma más segura de usar Diffchecker. ¡Obtén la app de Diffchecker Desktop: tus diffs nunca salen de tu computadora!Obtener Desktop

Poisson versus NB

Creado el mes pasadoEl diff nunca expira

12 eliminaciones

Líneas
Total
Eliminado

Caracteres
Total
Eliminado

Para continuar usando esta función, actualice a Diffchecker Pro Ver precios

165 líneas

22 adiciones

Líneas
Total
Añadido

Caracteres
Total
Añadido

Para continuar usando esta función, actualice a Diffchecker Pro Ver precios

170 líneas

# ---- fit poisson model with optim() ----

# ---- fit negative binomial model with optim() ----

# packages

library(tidyverse) # for glimpse() and ggplot()

# data; see ?crdata::holland2015

holland <- crdata::holland2015 |>

filter(city == "santiago")

# formula corresponds to model 1 for each city in holland (2015) table 2

f <- operations ~ lower + vendors + budget + population

# ---- create a function to fit the model ----

# log-likelihood function

# log-likelihood function (NB2: Var(y|X)=mu + mu^2/theta)

poisson_ll <- function(beta, y, X) {

nb_ll <- function(par, y, X) {

beta <- par[1:ncol(X)]

theta <- par[ncol(X) + 1]

linpred <- X%*%beta # perhaps denoted eta

lambda <- exp(linpred)

mu <- exp(linpred)

ll <- sum(dpois(y, lambda = lambda, log = TRUE))

ll <- sum(dnbinom(y, mu = mu, size = theta, log = TRUE))

return(ll)

}

# function to fit model

est_poisson <- function(f, data) {

est_nb <- function(f, data) {

# make X and y

mf <- model.frame(f, data = data)

X <- model.matrix(f, data = mf)

y <- model.response(mf)

# create starting values

par_start <- rep(0, ncol(X))

par_start <- c(rep(0, ncol(X)), 1) # append theta start (>0)

# run optim()

est <- optim(par_start,

fn = poisson_ll,

fn = nb_ll,

y = y,

X = X,

hessian = TRUE, # for SEs!

control = list(fnscale = -1),

method = "BFGS")

# check convergence; print warning if not

if (est$convergence != 0) print("Model did not converge!")

# create list of objects to return

res <- list(beta_hat = est$par,

k <- ncol(X)

var_hat = solve(-est$hessian))

res <- list(beta_hat = est$par[1:k],

var_hat = solve(-est$hessian)[1:k, 1:k],

theta_hat = est$par[k + 1],

var_theta_hat = solve(-est$hessian)[k + 1, k + 1])

# return the list

return(res)

}

# fit model

fit <- est_poisson(f, data = holland)

fit <- est_nb(f, data = holland)

print(fit, digits = 2) # print estimates w/ reasonable digits

# ---- compute the expected value given X_c ----

# create chosen values for X

# note 1: naming columns helps a bit later

# note 2: can also do with f, model.matrix(..., newdata = ...)

X_c <- cbind(

"constant" = 1, # intercept

"lower" = median(holland$lower),

"vendors" = median(holland$vendors),

"budget" = median(holland$budget),

"population" = median(holland$population)

)

# function to compute qi

ev_fn <- function(beta, X) {

exp(X%*%beta)

}

# invariance property

ev_hat <- ev_fn(fit$beta_hat, X_c)

# delta method

library(numDeriv) # for grad()

grad_ev <- grad(

func = ev_fn, # what function are we taking the derivative of?

x = fit$beta_hat, # what variable(s) are we taking the derivative w.r.t.?

X = X_c) # what other values are needed?

se_ev_hat <- sqrt(grad_ev %*% fit$var_hat %*% grad_ev)

# ---- compute the ev given X_c (w/ range of values) ----

# create chosen values for X

X_c <- cbind(

"constant" = 1, # intercept

"lower" = seq(min(holland$lower), max(holland$lower), by = 1),

"vendors" = median(holland$vendors),

"budget" = median(holland$budget),

"population" = median(holland$population)

)

# containers for estimated quantities of interest and ses

ev_hat <- numeric(nrow(X_c))

se_ev_hat <- numeric(nrow(X_c))

# loop over each row of X_c and compute qi and se

for (i in 1:nrow(X_c)) { # for the ith row of X...

Diferencias guardadas

Texto original

Abrir archivo

# ---- fit poisson model with optim() ----

# packages
library(tidyverse)  # for glimpse() and ggplot()

# data; see ?crdata::holland2015
holland <- crdata::holland2015 |>
  filter(city == "santiago")

# formula corresponds to model 1 for each city in holland (2015) table 2
f <- operations ~ lower + vendors + budget + population

# ---- create a function to fit the model ----

# log-likelihood function
poisson_ll <- function(beta, y, X) {
  linpred <- X%*%beta  # perhaps denoted eta
  lambda <- exp(linpred) 
  ll <- sum(dpois(y, lambda = lambda, log = TRUE))
  return(ll)
}

# function to fit model
est_poisson <- function(f, data) {
  
  # make X and y
  mf <- model.frame(f, data = data)
  X <- model.matrix(f, data = mf)
  y <- model.response(mf)
  
  # create starting values
  par_start <- rep(0, ncol(X))
  
  # run optim()
  est <- optim(par_start, 
               fn = poisson_ll, 
               y = y,
               X = X,
               hessian = TRUE, # for SEs!
               control = list(fnscale = -1),
               method = "BFGS") 
  
  # check convergence; print warning if not
  if (est$convergence != 0) print("Model did not converge!")
  
  # create list of objects to return
  res <- list(beta_hat = est$par,
              var_hat = solve(-est$hessian))
  
  # return the list
  return(res)
}

# fit model
fit <- est_poisson(f, data = holland)
print(fit, digits = 2)  # print estimates w/ reasonable digits

# ---- compute the expected value given X_c ----

# create chosen values for X
# note 1: naming columns helps a bit later
# note 2: can also do with f, model.matrix(..., newdata = ...)
X_c <- cbind(
  "constant" = 1, # intercept
  "lower"      = median(holland$lower), 
  "vendors"    = median(holland$vendors),
  "budget"     = median(holland$budget),
  "population" = median(holland$population)
)

# function to compute qi
ev_fn <- function(beta, X) {
  exp(X%*%beta)
}

# invariance property
ev_hat <- ev_fn(fit$beta_hat, X_c)

# delta method
library(numDeriv)  # for grad()
grad_ev <- grad(
  func = ev_fn,     # what function are we taking the derivative of?
  x = fit$beta_hat, # what variable(s) are we taking the derivative w.r.t.?
  X = X_c)          # what other values are needed?
se_ev_hat <- sqrt(grad_ev %*% fit$var_hat %*% grad_ev)

# ---- compute the ev given X_c (w/ range of values) ----

# create chosen values for X
X_c <- cbind(
  "constant" = 1, # intercept
  "lower"      = seq(min(holland$lower), max(holland$lower), by = 1), 
  "vendors"    = median(holland$vendors),
  "budget"     = median(holland$budget),
  "population" = median(holland$population)
)

# containers for estimated quantities of interest and ses
ev_hat <- numeric(nrow(X_c))
se_ev_hat <- numeric(nrow(X_c))

# loop over each row of X_c and compute qi and se
for (i in 1:nrow(X_c)) {   # for the ith row of X...
  # invariance property
  ev_hat[i] <- ev_fn(fit$beta_hat, X_c[i, ])
  # delta method
  grad_ev <- grad(
    func = ev_fn, 
    x = fit$beta_hat, 
    X = X_c[i, ]) 
  se_ev_hat[i] <- sqrt(grad_ev %*% fit$var_hat %*% grad_ev)
}

# put X_c, qi estimates, and se estimates in data frame
qi <- cbind(X_c, ev_hat, se_ev_hat) |>
  data.frame() |>
  glimpse()

# plot
ggplot(qi, aes(x = lower, y = ev_hat, 
               ymin = ev_hat - 1.64*se_ev_hat, 
               ymax = ev_hat + 1.64*se_ev_hat)) + 
  geom_ribbon() + 
  geom_line()

# ---- compute first difference ----

# make X_lo
X_lo <- cbind(
  "constant" = 1, # intercept
  "lower"      = quantile(holland$lower, probs = 0.25), 
  "vendors"    = median(holland$vendors),
  "budget"     = median(holland$budget),
  "population" = median(holland$population)
)

# make X_hi by modifying the relevant value of X_lo
X_hi <- X_lo
X_hi[, "lower"] <- quantile(holland$lower, probs = 0.75)

# function to compute first difference
fd_fn <- function(beta, hi, lo) {
  exp(hi%*%beta) - exp(lo%*%beta)
}

# invariance property
fd_hat <- fd_fn(fit$beta_hat, X_hi, X_lo)

# delta method
grad_fd <- grad(
  func = fd_fn, 
  x = fit$beta_hat, 
  hi = X_hi,
  lo = X_lo)  
se_fd_hat <- sqrt(grad_fd %*% fit$var_hat %*% grad_fd)

# estimated fd
fd_hat

# estimated se
se_fd_hat

# 90% ci
fd_hat - 1.64*se_fd_hat  # lower
fd_hat + 1.64*se_fd_hat  # upper

Texto modificado

Abrir archivo

# ---- fit negative binomial model with optim() ----

# packages
library(tidyverse)  # for glimpse() and ggplot()

# data; see ?crdata::holland2015
holland <- crdata::holland2015 |>
  filter(city == "santiago")

# formula corresponds to model 1 for each city in holland (2015) table 2
f <- operations ~ lower + vendors + budget + population

# ---- create a function to fit the model ----

# log-likelihood function (NB2: Var(y|X)=mu + mu^2/theta)
nb_ll <- function(par, y, X) {
  beta  <- par[1:ncol(X)]
  theta <- par[ncol(X) + 1]
  linpred <- X%*%beta  # perhaps denoted eta
  mu <- exp(linpred) 
  ll <- sum(dnbinom(y, mu = mu, size = theta, log = TRUE))
  return(ll)
}

# function to fit model
est_nb <- function(f, data) {
  
  # make X and y
  mf <- model.frame(f, data = data)
  X <- model.matrix(f, data = mf)
  y <- model.response(mf)
  
  # create starting values
  par_start <- c(rep(0, ncol(X)), 1)  # append theta start (>0)
  
  # run optim()
  est <- optim(par_start, 
               fn = nb_ll, 
               y = y,
               X = X,
               hessian = TRUE, # for SEs!
               control = list(fnscale = -1),
               method = "BFGS") 
  
  # check convergence; print warning if not
  if (est$convergence != 0) print("Model did not converge!")
  
  # create list of objects to return
  k <- ncol(X)
  res <- list(beta_hat = est$par[1:k],
              var_hat = solve(-est$hessian)[1:k, 1:k], 
              theta_hat = est$par[k + 1], 
              var_theta_hat = solve(-est$hessian)[k + 1, k + 1])
  
  # return the list
  return(res)
}

# fit model
fit <- est_nb(f, data = holland)
print(fit, digits = 2)  # print estimates w/ reasonable digits

# ---- compute the expected value given X_c ----

# function to compute qi
ev_fn <- function(beta, X) {
  exp(X%*%beta)
}

# invariance property
ev_hat <- ev_fn(fit$beta_hat, X_c)

# ---- compute the ev given X_c (w/ range of values) ----

# containers for estimated quantities of interest and ses
ev_hat <- numeric(nrow(X_c))
se_ev_hat <- numeric(nrow(X_c))

# put X_c, qi estimates, and se estimates in data frame
qi <- cbind(X_c, ev_hat, se_ev_hat) |>
  data.frame() |>
  glimpse()

# plot
ggplot(qi, aes(x = lower, y = ev_hat, 
               ymin = ev_hat - 1.64*se_ev_hat, 
               ymax = ev_hat + 1.64*se_ev_hat)) + 
  geom_ribbon() + 
  geom_line()

# ---- compute first difference ----

# make X_hi by modifying the relevant value of X_lo
X_hi <- X_lo
X_hi[, "lower"] <- quantile(holland$lower, probs = 0.75)

# function to compute first difference
fd_fn <- function(beta, hi, lo) {
  exp(hi%*%beta) - exp(lo%*%beta)
}

# invariance property
fd_hat <- fd_fn(fit$beta_hat, X_hi, X_lo)

# delta method
grad_fd <- grad(
  func = fd_fn, 
  x = fit$beta_hat, 
  hi = X_hi,
  lo = X_lo)  
se_fd_hat <- sqrt(grad_fd %*% fit$var_hat %*% grad_fd)

# estimated fd
fd_hat

# estimated se
se_fd_hat

# 90% ci
fd_hat - 1.64*se_fd_hat  # lower
fd_hat + 1.64*se_fd_hat  # upper