Comparar texto

Encontre a diferença entre dois arquivos de texto

Editor live

Recolher inalteradas

Sem quebra de linha

Layout

Nível de detalhe

Realce de sintaxe

Diffchecker Desktop A maneira mais segura de usar o Diffchecker. Obtenha o aplicativo Diffchecker Desktop: seus diffs nunca saem do seu computador!Obter Desktop

Poisson versus NB

Criado há 9 mesesO diff nunca expira

13 remoções

Linhas
Total
Removido

Caracteres
Total
Removido

Para continuar usando este recurso, atualize para Diffchecker Pro Ver preços

164 linhas

21 adições

Linhas
Total
Adicionado

Caracteres
Total
Adicionado

Para continuar usando este recurso, atualize para Diffchecker Pro Ver preços

168 linhas

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

Diferenças salvas

Texto original

Abrir arquivo

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

Abrir arquivo

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