Nothing
# Multivariate negative binomial model with robust estimation of regression coefficients.
# This function fits a multivariate negative binomial model by Maximum Likelihood and calculates robust standard errors
# of the regression coefficients.
multinbmod<-
function (formula, data, id, offset, start.coef = NULL,
start.phi = NULL,control=list())
# initial settings for optimisation
{ inv <- match.call()
if (missing(data))
data <- environment(formula)
if (!is.list(control))
stop("control must be a list")
mf <- match.call(expand.dots = FALSE)
mf$start.coef <- mf$start.phi<- mf$control<-NULL
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, environment(formula))
mt <- attr(mf, "terms")
xvars <- as.character(attr(mt, "variables"))[-1]
if ((yvar <- attr(mt, "response")) > 0)
xvars <- xvars[-yvar]
xlev <- if (length(xvars) > 0) {
xlev <- lapply(mf[xvars], levels)
xlev[!sapply(xlev, is.null)]
}
X <- if (!is.empty.model(mt))
model.matrix(mt, mf, contrasts)
p <- NCOL(X)
Y <- model.response(mf, "numeric")
id <- mf$"(id)"
no.id <- (missing(id) || is.null(id))
if (no.id) {
warning("id variable is missing!")
return(NULL)
}
offset <- model.offset(mf)
if(is.null(offset)){
offset <- rep(1,length(Y))}
if (!is.null(offset)&& length(offset) != NROW(Y))
stop(paste("Number of offsets is", length(offset), ", should equal",
NROW(Y), "(number of observations)"))
if (is.null(start.coef)) {
start.coef <- numeric(p)
start.coef[1] <- log(mean(Y + 0.5))
}
else {
if (length(start.coef) != p)
stop("beta.start has wrong length")
}
if (is.null(start.phi)) {
start.phi <- 0.5
}
else {
if (length(start.phi) != 1)
stop("start.sigma has wrong length")
if (start.phi<=0)
stop("start.phi must be positive")
}
fit<-multinb.fit(y=Y, x=X, offset, id=id, start.par=c(start.coef,start.phi),control)
res<-list()
res$converged <- fit$"converged?"
res$coefficients<-fit$"estimated regression coefficients"
res$model.coef.se<- fit$"se from model"
res$robust.coef.se<- fit$"robust se"
res$robust.t.values<- fit$"robust t-values"
res$mle.phi <- fit$"estimated phi"
res$phi.se<- fit$"se(phi)"
res$minus2.loglik<-fit$"-2 x loglikelihood"
res$iterations<- fit$iterations
res$call<- inv
class(res) <- "multinbmod"
return(res)
}
#Summary of a multinbmod object
summary.multinbmod<-function(object,...)
{
TAB<-cbind(Estimate=object$coefficients,
ModelSE = object$model.coef.se,
RobustSE = object$robust.coef.se,
Robust.t= object$robust.t.values)
res<- list( call=object$call,converged=object$converged,coefficients=TAB,MLE_of_phi=object$mle.phi,SE_of_phi=object$phi.se,
minus2.loglik=object$minus2.loglik,iterations=object$iterations)
class(res) <- "summary.multinbmod"
res
}
#Multivariate negative binomial model with robust estimation of regression coefficients
#This function is called by "multinbmod", but it can also be called directly
multinb.fit<-
function(y, x, offset=1, id, start.par,control=list())
{
if (!is.list(control))
stop("control must be a list")
np<-NCOL(x)+1
if (missing(start.par)) {
start.par <- numeric(np)
start.par[1] <- log(mean(y + 0.5))
start.par[np]<- 0.5
}
else{
if (length(start.par) != np)
stop("start.par has wrong length")
if (start.par[np]<=0)
stop("starting value of phi must be positive")}
no.id <- (missing(id) || is.null(id))
if (no.id) {
warning("id variable is missing!")
return(NULL)
}
if (offset!=1 && length(offset) != NROW(y))
stop(paste("Number of offsets is", length(offset), ", should equal",
NROW(y), "(number of observations)"))
tapplys.fun<-function(x,y){
res<-tapply(x,y,sum)
return(res[order(unique(y))])}
ydot<-tapplys.fun(y, id)
mnbloglik <- function(param, y, x, id,ydot)
{
p <- length(param)
beta <- param[1:(p - 1)]
a <- param[p]
mu <- offset*exp(x %*% beta)
mudot <- tapplys.fun(mu, id)
minlog<- -1 * sum(lgamma(a^(-1) + ydot) - lgamma(a^(-1)) - (a^(-1) +
ydot) * log(1 + mudot * a) + ydot * log(a) + ydot *
log(mudot)) - sum(tapplys.fun(y * log(mu), id) - log(
mudot) * ydot)
return(minlog)
}
grr <- function(param, y, x,id,ydot)
{
p <- length(param)
beta <- param[1:(p - 1)]
a <- param[p]
mu <- offset*exp(x %*% beta)
mudot <- tapplys.fun(mu, id)
derivbeta <- rep(0, p)
for(i in 1:(p - 1)) {
derivbeta[i] <- -1*sum(tapplys.fun(y*x[,i],id) - tapplys.fun(x[,i]*mu,id)*(1+ydot*a)/(1+mudot*a))
}
theta<- 1/a
derivbeta[p] <- (1/a^2)*sum(digamma(theta+ydot)-log(1+mudot/theta)+
(theta+ydot)*mudot/(theta*(theta+mudot))-ydot/theta-digamma(theta))
return(derivbeta)
}
hes<-function(param,y,x,id,ydot)
{
p<- length(param)
beta <- param[1:(p - 1)]
a <- param[p]
mu <- offset*exp(x %*% beta)
mudot <- tapplys.fun(mu, id)
hesmat<-matrix(numeric(p^2),p,p)
for (j in (1:(p-1))){
for (i in (1:(p-1))){
hesmat[i,j]<- -1*sum(-1*tapplys.fun(mu*x[,i]*x[,j],id)*(1+a*ydot)/(1+a*mudot) + tapplys.fun(mu*x[,i],id)*tapplys.fun(mu*x[,j],id)*a*
(1+a*ydot)/(1+a*mudot)^2)
}}
for (i in (1:(p-1))){
hesmat[i,p]<-hesmat[p,i]<- sum(tapplys.fun(mu*x[,i],id)*(ydot-mudot)/(1+mudot*a)^2)
}
hesmat[p,p]<- -1*sum(-ydot/a^2 + trigamma(a^(-1)+ydot)/a^4+2*digamma(a^(-1)+ydot)/a^3-trigamma(a^(-1))/a^4-
2*digamma(a^(-1))/a^3 + 2*mudot/(a^2*(1+mudot*a))-2*log(1+mudot*a)/a^3+(a^(-1)+ydot)*mudot^2/(1+mudot*a)^2)
return(hesmat)
}
res<- nlminb(start.par, mnbloglik,grr,lower = c(rep( - Inf, (np - 1)), 0.01),
y = y, x = x, id = id,ydot=ydot,control=control)
betamle <- res$par[1:(np - 1)]
names(betamle)<-names(as.data.frame(x))
amle <- as.numeric(res$par[np])
muest <- offset*exp(x %*% betamle)
mudotest <- tapplys.fun(muest, id)
svec <- (y - muest)
Nid <- length(unique(id))
var0<- var1 <- matrix(numeric((np - 1)^2), (np - 1))
for(i in 1:Nid) {
eks <- x[(pick <- id == unique(id)[i]), , drop = F]
eks <- as.matrix(eks)
mus <- muest[pick]
ders <- eks * mus
if(length(mus) > 1) {
visinv <- diag(1/mus) - (amle/(1 + amle * sum(mus))) * rep(1, length(mus)) %*% t(rep(1, length(
mus)))
var0 <- var0 + t(ders) %*% visinv %*% ders
vsan<- svec[pick] %*% t(svec[pick])
var1 <- var1 + t(ders) %*% visinv %*% vsan %*% visinv %*% ders
}
else {
visinv <- 1/(mus * (1 + amle * mus))
var0 <- var0 + visinv * (t(ders) %*% ders)
var1 <- var1 + (t(ders) %*% ders) * (visinv^2 * (svec[
pick])^2)
}
}
modcov<-solve(var0)
robcovbetamle <- solve(var0) %*% var1 %*% solve(var0)
robsebetamle <- sqrt(diag(robcovbetamle))
names(robsebetamle)<-names(betamle)
tvalues <- betamle/robsebetamle
modsebemle<- sqrt(diag(modcov))
names(modsebemle)<-names(betamle)
amle<- res$par[np]
varesta<- 1/(-1*sum(-ydot/amle^2 + trigamma(amle^(-1)+ydot)/amle^4+2*digamma(amle^(-1)+ydot)/amle^3-trigamma(amle^(-1))/amle^4-
2*digamma(amle^(-1))/amle^3 + 2*mudotest/(amle^2*(1+mudotest*amle))-2*log(1+mudotest*amle)/amle^3+(amle^(-1)+ydot)*mudotest^2/(1+mudotest*amle)^2) )
myres <- list(betamle,modsebemle,robsebetamle, tvalues,modcov,robcovbetamle,
amle, sqrt(varesta), 2*(res$
objective+sum(lgamma(y+1))), ifelse(res$convergence==0,TRUE,FALSE),res$iterations)
names(myres) <- c("estimated regression coefficients","se from model","robust se", "robust t-values","covariance of beta estimates from model",
"robust covariance of beta estimates","estimated phi","se(phi)","-2 x loglikelihood","converged?",
"iterations")
myres
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.