Nothing
derivative_tD_n_W_inv_n <- function(r1, mu_n, eta_n, X_n, tD_n, W_n_inv, A_n_inv, A_n_inv_d, Working_corr, linkd1, linkd2, varfd1, cluster_condition, cluster_size)
{
if(cluster_condition==1)
Working_corr_inv <- solve(Working_corr)
else
Working_corr_inv <- solve(Working_corr[1:cluster_size, 1:cluster_size])
tmp <- (A_n_inv_d*varfd1(mu_n)*linkd1(eta_n)*X_n[,r1])%*%Working_corr_inv%*%A_n_inv
##r1: component of beta wrt derivative has been taken
t(linkd2(eta_n)*X_n[,r1]*X_n)%*%W_n_inv - .5*tD_n%*%( tmp+t(tmp) )
}
derivative_D_n <- function(r1, eta_n, X_n, linkd2)
{
##r1: component of beta wrt derivative has been taken
linkd2(eta_n)*X_n[,r1]*X_n
}
k_r1r2_r3_fun <- function(p, mu_n, eta_n, X_n, tD_n, Cov_n, tD_W_inv_n, W_n_inv, A_n_inv, A_n_inv_d, Working_corr, linkd1, linkd2, varfd1, cluster_condition, cluster_size)
{
k_r1r2_r3_vals <- lapply(1:p, function(PARS) Reduce("+", lapply(1:length(eta_n),
function(x) derivative_tD_n_W_inv_n(
r1=PARS,
mu_n=mu_n[[x]],
eta_n=eta_n[[x]],
X_n=X_n[[x]],
tD_n=tD_n[[x]],
W_n_inv=W_n_inv[[x]],
A_n_inv=A_n_inv[[x]],
A_n_inv_d=A_n_inv_d[[x]],
Working_corr=Working_corr,
linkd1=linkd1, linkd2=linkd2, varfd1=varfd1,
cluster_condition=cluster_condition,
cluster_size=cluster_size[x])%*%Cov_n[[x]]%*%t(tD_W_inv_n[[x]]))))
k_r1r2_r3_vals <- lapply(1:p, function(z) Reduce("rbind",lapply(k_r1r2_r3_vals,"[",z,)))
k_r1r2_r3_vals
}
k_r1r2r3_fun <- function(p, mu_n, eta_n, X_n, D_n, tD_n, Cov_n, tD_W_inv_n, W_n_inv, A_n_inv, A_n_inv_d, Working_corr, linkd1, linkd2, varfd1, cluster_condition, cluster_size=NULL)
{
k_r1r2r3_vals1 <- lapply(1:p, function(PARS) Reduce("+", lapply(1:length(eta_n),
function(x) -derivative_tD_n_W_inv_n(
r1=PARS,
mu_n=mu_n[[x]],
eta_n=eta_n[[x]],
X_n=X_n[[x]],
tD_n=tD_n[[x]],
W_n_inv=W_n_inv[[x]],
A_n_inv=A_n_inv[[x]],
A_n_inv_d=A_n_inv_d[[x]],
Working_corr=Working_corr,
linkd1=linkd1, linkd2=linkd2, varfd1=varfd1,
cluster_condition=cluster_condition,
cluster_size=cluster_size[x])%*%D_n[[x]])))
k_r1r2r3_vals2 <- lapply(1:p, function(PARS) Reduce("+", lapply(1:length(eta_n),
function(x) -tD_W_inv_n[[x]]%*%derivative_D_n(
r1=PARS,
eta_n=eta_n[[x]],
X_n=X_n[[x]], linkd2=linkd2))))
k_r1r2r3_vals1 <- lapply(1:p, function(z){ tmp <- Reduce("rbind",lapply(k_r1r2r3_vals1,"[",z,)); tmp+t(tmp)})
k_r1r2r3_vals2 <- lapply(1:p, function(z) Reduce("rbind",lapply(k_r1r2r3_vals2,"[",z,)))
k_r1r2r3_vals <- mapply("+", k_r1r2r3_vals1, k_r1r2r3_vals2, SIMPLIFY = FALSE)
k_r1r2r3_vals
}
##Link functions and their derivatives
##Supported functions are the same as in package gee (to date 4/14/17)
link_funcs <- function(link_name)
{
##link function
link <- switch(link_name,
Identity = function(x) x,
Logarithm = function(x) exp(x),
Logit = function(x) exp(x)/(1+exp(x)),
Reciprocal = function(x) 1/x,
Probit = function(x) pnorm(x),
Cloglog = function(x) 1-exp(-exp(x))
)
##first derivative of the link function
linkd1 <- switch(link_name,
Identity = function(x) rep(1, length(x)),
Logarithm = function(x) exp(x),
Logit = function(x) exp(x)/(1+exp(x))^2,
Reciprocal = function(x) -1/x^2,
Probit = function(x) dnorm(x),
Cloglog = function(x) exp(x-exp(x))
)
##second derivative of the link function
linkd2 <- switch(link_name,
Identity = function(x) rep(0, length(x)),
Logarithm = function(x) exp(x),
Logit = function(x) -exp(x)*(exp(x)-1)/(1+exp(x))^3,
Reciprocal = function(x) 2/x^3,
Probit = function(x) -dnorm(x)*x,
Cloglog = function(x) exp(x-exp(x))*(1-exp(x))
)
list(link=link, linkd1=linkd1, linkd2=linkd2)
}
##Variance functions and their derivatives
##Supported functions are the same as in package gee (to date 4/14/17)
var_funcs <- function(var_name)
{
##variance function
varf <- switch(var_name,
Gaussian = function(x) rep(1, length(x)),
Poisson = function(x) x,
Binomial = function(x) x*(1-x),
Gamma = function(x) x^2
)
##first derivative of the variance function
varfd1 <- switch(var_name,
Gaussian = function(x) rep(0, length(x)),
Poisson = function(x) rep(1, length(x)),
Binomial = function(x) 1-2*x,
Gamma = function(x) 2*x
)
list(varf=varf, varfd1=varfd1)
}
BCgee <- function(fit)
{
##check for class of fit
if(!any(class(fit)%in%"gee"))
stop("fit must be retrieved from function gee in package gee.")
links <- link_funcs(fit$model$link)
linkf <- links$link
linkd1 <- links$linkd1
linkd2 <- links$linkd2
vars <- var_funcs(fit$model$varfun)
varf <- vars$varf
varfd1 <- vars$varfd1
mnames <- c("", "formula", "data", "offset", "subset", "na.action", "id")
cnames <- names(fit$call)
cnames <- cnames[match(mnames, cnames, 0)]
mcall <- fit$call[cnames]
mcall[[1]] <- as.name("model.frame")
model.f <- eval(mcall, parent.frame())
Terms <- attributes(model.f)$terms
contrasts <- model.extract(model.f, contrasts)
X <- model.matrix(Terms, model.f, contrasts)
offset <- model.offset(model.f)
if(is.null(offset))
offset <- 0
y <- as.matrix(model.extract(model.f, "response"))
weights <- rep(1, length(y))
flag_weights <- 0
if(ncol(y)==2)
{
flag_weights <- 1
weights <- as.vector(y%*%c(1, 1))
y <- y[,1]
}
##linear predictor
eta_n <- tapply(fit$linear.predictors + offset, fit$id, function(x) x )
##number of clusters
n_cluster <- length(eta_n)
##fitted values
mu_n <- lapply(eta_n, "linkf" )
##
weights_n <- by(weights, fit$id, function(y) y )
##design matrix
X_n <- by(X, fit$id, function(y) as.matrix(y) )
##D matrix
D_n <- lapply(1:n_cluster, function(x) linkd1(eta_n[[x]])*X_n[[x]] )
tD_n <- lapply(D_n, "t")
#
if(!flag_weights)
{
A_n_inv <- lapply(mu_n, function(y) (fit$scale^(-1/2))*diag(varf(y)^(-1/2)))
A_n_inv_d <- lapply(mu_n, function(y) (fit$scale^(-1/2))*diag(varf(y)^(-3/2)))
}else{
lst_mu_wei <- mapply(function(x,y) list(mu=x,w=y), mu_n, weights_n, SIMPLIFY=FALSE)
A_n_inv <- lapply(lst_mu_wei, function(Z) (fit$scale^(-1/2))*diag(Z$w^(1/2)*varf(Z$mu)^(-1/2)))
A_n_inv_d <- lapply(lst_mu_wei, function(Z) (fit$scale^(-1/2))*diag(Z$w^(1/2)*varf(Z$mu)^(-3/2)))
}
Working_corr <- fit$working.correlation
##account for non-constant cluster size
cluster_size <- table(fit$id)
cluster_condition <- length(unique(cluster_size))
if(cluster_condition==1)
{
Working_corr_inv <- solve(Working_corr)
W_n_inv <- lapply(1:n_cluster, function(x) A_n_inv[[x]]%*%Working_corr_inv%*%A_n_inv[[x]])
}else{
W_n_inv <- lapply(1:n_cluster, function(x) A_n_inv[[x]]%*%solve(Working_corr[1:cluster_size[x],1:cluster_size[x]])%*%A_n_inv[[x]])
}
tD_W_inv_n <- lapply(1:n_cluster, function(x) tD_n[[x]]%*%W_n_inv[[x]])
residuals_n <- (fit$y-fit$fitted.values*weights)/weights
residuals_n <- tapply(residuals_n, fit$id, function(x) x )
Cov_n <- lapply(1:n_cluster, function(x) residuals_n[[x]]%*%t(residuals_n[[x]]))
###The sandwich (asymptotic variance of estimator of regression parameters)
Sandwich <- fit$robust.variance
#The bread in the sandwich
B <- fit$naive.variance
p <- length(coefficients(fit))
k_r1r2_r3_vals <- k_r1r2_r3_fun(p, mu_n, eta_n, X_n, tD_n, Cov_n, tD_W_inv_n, W_n_inv, A_n_inv, A_n_inv_d,
Working_corr, linkd1, linkd2, varfd1, cluster_condition, cluster_size)
k_r1r2r3_vals <- k_r1r2r3_fun(p, mu_n, eta_n, X_n, D_n, tD_n, Cov_n, tD_W_inv_n, W_n_inv, A_n_inv, A_n_inv_d,
Working_corr, linkd1, linkd2, varfd1, cluster_condition, cluster_size)
#first piece of the bias
b1 <- B%*%unlist( lapply(k_r1r2_r3_vals, function(x) sum(x*B)) )
#second piece of the bias
b2 <- .5*B%*%unlist( lapply(k_r1r2r3_vals, function(x) sum(x*Sandwich)) )
result <- fit
##bias estimates
result$bias <- as.numeric(b1+b2)
attr(result$bias,"names") <- attributes(coefficients(fit))$names
##bias corrected estimates
result$coefficients <- result$coefficients - result$bias
##updated linear predictors
result$linear.predictors <- as.numeric(X%*%result$coefficients)# + offset)
##updated fitted
result$fitted.values <- as.numeric(linkf(result$linear.predictors))
##updated residuals
result$residuals <- as.numeric(result$y-result$fitted.values*weights)/weights
##updated naive variance
eta_n_upd <- tapply(result$linear.predictors, result$id, function(x) x )
D_n_upd <- lapply(1:n_cluster, function(x) linkd1(eta_n_upd[[x]])*X_n[[x]] )
tD_n_upd <- lapply(D_n_upd, "t")
if(!flag_weights)
{
A_n_inv_upd <- by(result$fitted.values, fit$id, function(y) (fit$scale^(-.5))*diag(varf(y)^(-.5)))
}else{
mu_n_upd <- by(result$fitted.values, fit$id, function(y) y )
lst_mu_wei_upd <- mapply(function(x,y) list(mu=x,w=y), mu_n_upd, weights_n, SIMPLIFY=FALSE)
A_n_inv_upd <- lapply(lst_mu_wei_upd, function(Z) (fit$scale^(-1/2))*diag(Z$w^(1/2)*varf(Z$mu)^(-1/2)))
}
if(cluster_condition==1)
{
W_n_inv_upd <- lapply(1:n_cluster, function(x) A_n_inv_upd[[x]]%*%Working_corr_inv%*%A_n_inv_upd[[x]])
}else{
W_n_inv_upd <- lapply(1:n_cluster, function(x) A_n_inv_upd[[x]]%*%solve(Working_corr[1:cluster_size[x],1:cluster_size[x]])%*%A_n_inv_upd[[x]])
}
tD_W_inv_n_upd <- lapply(1:n_cluster, function(x) tD_n_upd[[x]]%*%W_n_inv_upd[[x]])
result$naive.variance <- solve(Reduce("+", lapply(1:n_cluster, function(x) tD_W_inv_n_upd[[x]]%*%D_n_upd[[x]])))
##updated robust variance
residuals_n_upd <- tapply(result$residuals, result$id, function(x) x )
Cov_n_upd <- lapply(1:n_cluster, function(x) residuals_n_upd[[x]]%*%t(residuals_n_upd[[x]]))
V_upd <- Reduce("+", lapply(1:n_cluster, function(x) tD_W_inv_n_upd[[x]]%*%Cov_n_upd[[x]]%*%t(tD_W_inv_n_upd[[x]])))
result$robust.variance <- result$naive.variance%*%V_upd%*%result$naive.variance
attr(result, "class") <- c("BCgee")
result
}
print.BCgee <- function(x, digits = NULL, quote = FALSE, prefix = "", ...)
{
if(is.null(digits)) digits <- options()$digits else options(digits =
digits)
# cat("\n BCgee: BIAS CORRECTED ESTIMATES OF REGRESSION PARAMATERS OF ")
# cat("\n GENERALIZED LINEAR MODELS FOR DEPENDENT DATA ")
cat("\n BCgee: Bias correction based on estimates retrieved from gee package")
cat("\n ", x$version, "\n")
cat("\nModel:\n")
cat(" Link: ", x$model$link, "\n")
cat(" Variance to Mean Relation:", x$model$varfun, "\n")
if(!is.null(x$model$M))
cat(" Correlation Structure: ", x$model$corstr, ", M =", x$
model$M, "\n")
else cat(" Correlation Structure: ", x$model$corstr, "\n")
cat("\nCall:\n")
dput(x$call) # cat("\nTerms:\n")
### ys <- matrix(rep(as.matrix(x$id, ncol = 1), 5), ncol = 5)
ys <- matrix(rep(matrix(x$id, ncol = 1), 5), ncol = 5)
ys[, 2] <- x$y
ys[, 3] <- x$linear.predictors
ys[, 4] <- x$fitted.values
ys[, 5] <- x$residuals
dimnames(ys) <- list(1:length(x$y), c("ID", "Y", "LP", "fitted",
"Residual")) # cat("\nFitted Values:\n")
cat("\nNumber of observations : ", x$nobs, "\n")
cat("\nMaximum cluster size : ", x$max.id, "\n")
nas <- x$nas
if(any(nas))
cat("\n\nCoefficients: (", sum(nas),
" not defined because of singularities)\n", sep = "")
else cat("\n\nBias-corrected coefficients:\n")
print(x$coefficients, digits = digits)
cat("\nEstimated Scale Parameter: ", format(round(x$scale, digits)))
cat("\nNumber of Iterations: ", x$iterations)
nc <- min(4, nrow(x$working.correlation))
cat("\n\nWorking Correlation\n")
print(x$working.correlation[1:nc, 1:nc], digits = digits)
cat("\n\nReturned Error Value:\n")
print(x$error)
invisible(x)
}
print.summary.BCgee <- function(x, digits = NULL, quote = FALSE, prefix = "", ... )
{
if(is.null(digits))
digits <- options()$digits
else options(digits = digits)
# cat("\n",x$title)
# cat("\n",x$version,"\n")
cat("\n BCgee: Bias correction based on estimates retrieved from gee package")
cat("\n ", x$version, "\n")
cat("\nModel:\n")
cat(" Link: ",x$model$link,"\n")
cat(" Variance to Mean Relation:",x$model$varfun,"\n")
if(!is.null(x$model$M))
cat(" Correlation Structure: ",x$model$corstr,", M =",x$model$M,"\n")
else cat(" Correlation Structure: ",x$model$corstr,"\n")
cat("\nCall:\n")
dput(x$call)
cat("\nSummary of Residuals:\n")
print(x$residual.summary, digits = digits)
nas <- x$nas
### if(any(nas))
if(!is.null(nas) && any(nas))
cat("\n\nCoefficients: (", sum(nas),
" not defined because of singularities)\n", sep = "")
else cat("\n\nCoefficients:\n")
print(x$coefficients, digits = digits)
cat("\nEstimated Scale Parameter: ", format(round(x$scale,digits)))
cat("\nNumber of Iterations: ", x$iterations)
cat("\n\nWorking Correlation\n")
print(x$working.correlation,digits=digits)
if(!is.null(x$naive.correlation)) {
cat("\n\nNaive Correlation of Estimates:\n")
print(x$naive.correlation)
}
if(!is.null(x$robust.correlation)) {
cat("\n\nRobust Correlation of Estimates:\n")
print(x$robust.correlation)
}
invisible(x)
}
summary.BCgee <- function(object, correlation = TRUE, ...)
{
coef <- object$coefficients
resid <- object$residuals
p <- object$rank
if(is.null(p))
p <- sum(!is.na(coef))
if(!p) {
warning("This model has zero rank --- no summary is provided")
return(object)
}
nas <- is.na(coef)
cnames <- names(coef[!nas])
coef <- matrix(rep(coef[!nas], 5), ncol = 5)
dimnames(coef) <- list(cnames, c("Estimate",
"Naive S.E.", "Naive z",
"Robust S.E.", "Robust z"))
rse <- sqrt(diag(object$robust.variance))
nse <- sqrt(diag(object$naive.variance))
coef[,2] <- nse
coef[,3] <- coef[,1]/coef[,2]
coef[,4] <- rse
coef[,5] <- coef[,1]/coef[,4]
summary <- list()
summary$call <- object$call
summary$version <- object$version
summary$nobs <- object$nobs
summary$residual.summary <- quantile(as.vector(object$residuals))
names(summary$residual.summary) <- c("Min", "1Q", "Median", "3Q", "Max")
summary$model<- object$model
summary$title <- object$title
summary$coefficients <- coef
summary$working.correlation <- object$working.correlation
summary$scale <- object$scale
summary$error <- paste("Error code was", object$error)
summary$working.correlation <- object$working.correlation
summary$iterations <- object$iterations
if ( correlation ) {
## rob.var <- object$robust.variance
## nai.var <- object$naive.variance
## summary$robust.correlation <- rob.var /
## outer(sqrt(diag(rob.var)),sqrt(diag(rob.var)))
## dimnames(summary$robust.correlation) <- list(object$xnames,object$xnames)
## summary$naive.correlation <- nai.var /
## outer(sqrt(diag(nai.var)),sqrt(diag(nai.var)))
## dimnames(summary$naive.correlation) <- list(object$xnames,object$xnames)
}
attr(summary,"class") <- "summary.BCgee"
summary
}
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.