Comparar texto

Encuentra la diferencia entre dos archivos de texto

Editor live

Ocultar sin cambios

Sin ajuste de línea

Vista

Nivel de detalle

Resaltado de sintaxis

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 hace 9 mesesEl diff nunca expira

13 eliminaciones

Líneas
Total
Eliminado

Caracteres
Total
Eliminado

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

164 líneas

21 adiciones

Líneas
Total
Añadido

Caracteres
Total
Añadido

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

168 líneas

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

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

# 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)

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 <- 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 %*% fit$var_hat %*% grad)

# ---- 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 <- grad(

func = ev_fn,

x = fit$beta_hat,

X = X_c[i, ])

se_ev_hat[i] <- sqrt(grad %*% fit$var_hat %*% grad)

}

# 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 <- grad(

func = fd_fn,

x = fit$beta_hat,

hi = X_hi,

lo = X_lo)

se_fd_hat <- sqrt(grad %*% fit$var_hat %*% grad)

# 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

Diferencias guardadas

Texto original

Abrir archivo

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

# 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 <- 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 %*% fit$var_hat %*% grad)

# ---- 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 <- grad(
    func = ev_fn, 
    x = fit$beta_hat, 
    X = X_c[i, ]) 
  se_ev_hat[i] <- sqrt(grad %*% fit$var_hat %*% grad)
}

# 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 <- grad(
  func = fd_fn, 
  x = fit$beta_hat, 
  hi = X_hi,
  lo = X_lo)  
se_fd_hat <- sqrt(grad %*% fit$var_hat %*% grad)

# 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() ----

# 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_poisson(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 <- grad(
  func = fd_fn, 
  x = fit$beta_hat, 
  hi = X_hi,
  lo = X_lo)  
se_fd_hat <- sqrt(grad %*% fit$var_hat %*% grad)

# 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