# NOT RUN {
# generate data
set.seed(123)
n = 6; t = 3
id = rep(1:n, each = t)
rand.int = rep(rnorm(n, sd = 0.7), each = t)
group = rep(c(0,1), each = n*t/2)
time = rep(0:(t-1), n)
offset = rnorm(n*t, sd = 0.3)
beta = c(3, 0.3, 0.1)
X = model.matrix(~group + time)
mu = exp(X %*% beta + rand.int + offset)
y = rep(NA, n*t)
library(tweeDEseq)
for (i in 1:(n*t)) y[i] = rPT(1, mu = mu[i], D = 2, a = 0, max = 1000)
data.long = data.frame(y, group, time, id, offset)
rm(list = setdiff(ls(), 'data.long'))
# 1) Quick example (5 quadrature points, hessian and SEs not computed)
# estimate the model
fit1 = nbmixed(fixef.formula = y ~ group + time, id = data.long$id,
offset = data.long$offset, data = data.long, npoints = 5,
freq.updates = 200, hessian = FALSE, trace = TRUE)
# print summary:
summary(fit1, wald = FALSE)
# }
# NOT RUN {
# 2) Full computation, including hessian evaluation and using more quadrature points
# estimate the model
fit2 = nbmixed(fixef.formula = y ~ group + time, id = data.long$id,
offset = data.long$offset, data = data.long, npoints = 10,
freq.updates = 200, hessian = TRUE, trace = TRUE)
# print and get summary:
results = summary(fit2, wald = TRUE)
ls(results)
# view table with estimates of regression coefficients, standard errors and Wald test:
results$coefficients
# }
Run the code above in your browser using DataCamp Workspace