Description Usage Arguments Author(s) References Examples
This is the driver for the rank-based GEE fit. It is flly described in Section 8.6 of Kloke and McKean (2014).
1 |
y |
Response vector. |
xmat |
Design matrix. |
center |
Center (Block) indicator. |
scores |
Scores, default is Wilcoxon. |
geemod |
Link function, currently, only option is liear model ("LM"). |
structure |
Working covariance type. |
substructure |
Variance-component type. |
med |
TRUE for symmetric type scores (score function odd about (1/2)); else quantile is computed accordingly. |
varcovst |
Type of standard errors. |
maxstp |
Maximum number of iterations (default is 50). |
eps |
Convergence criteria (default is 0.00001). |
hbrs |
FALSE is default resulting in rank based fit. If TRUE a wteighted HBR fit is computed. |
delta |
Window option for scale estimator. Default is .8. Must be strictly less than 1. |
hparm |
Option for outlier scale correction. Default is 2. |
Jospeph W. McKean
Kloke and McKean (2014), Nonparametrics Using R, Boca Raton: Chapman-Hall.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (y, xmat, center, scores = wscores, geemod = "LM", structure = "CS",
substructure = "MM", med = TRUE, varcovst = "var2", maxstp = 50,
eps = 1e-05, hbrs = FALSE, delta = 0.8, hparm = 2)
{
int = TRUE
xmat <- centerx(xmat)
n <- length(y)
one <- rep(1, n)
p <- length(xmat[1, ])
if (int) {
if (hbrs == FALSE) {
fit0 <- rfit(y ~ xmat, scores = scores)
}
if (hbrs == TRUE) {
fit0 <- hbrfit(y ~ xmat)
}
beta0 <- fit0$coef[2:(p + 1)]
y <- y - fit0$coef[1]
alphat <- fit0$coef[1]
}
else {
if (hbrs == FALSE) {
fit0 <- rfit(y ~ xmat, scores = scores)
yhatint <- fitted.values(fit0)
}
if (hbrs == TRUE) {
fit0 <- hbrfit(y ~ xmat)
yhatint <- y - fit0$resid
}
beta0 <- solve(t(xmat) %*% xmat) %*% t(xmat) %*% yhatint
}
colbeta <- matrix(rep(0, maxstp * 2 * p), ncol = 2 * p)
istp <- 0
iflag <- 0
ic <- 1
while (ic > 0) {
apbeta0 <- getAp(xmat, center, beta0, geemod)
S <- y - apbeta0
vee <- veemat(S, center, structure, substructure, scores)
vee12 <- as.matrix(vee$v12)
vee12inv <- as.matrix(vee$v12inv)
Dmat <- getD(xmat, center, beta0, geemod)
eitb <- vee12inv %*% (y - apbeta0)
Dmatb <- vee12inv %*% Dmat
wtstuff <- wtmat(Dmatb, eitb, med, scores, hbrs)
wts <- wtstuff$wts
wts <- as.vector(wts)
msbeta0 <- wtstuff$mb * vee12 %*% one
S <- y - apbeta0 - msbeta0
veeinv <- vee12inv %*% diag(wts) %*% vee12inv
part1 <- t(Dmat) %*% veeinv %*% Dmat
inc <- solve(part1) %*% t(Dmat) %*% veeinv %*% S
beta1 <- beta0 + inc
istp <- istp + 1
colbeta[istp, ] <- c(beta0, beta1)
chk <- sum(inc^2)/(sum(beta0^2) + 1e-05)
if (chk < eps) {
iflag <- 1
ic <- 0
betag <- beta1
}
else {
if (istp >= maxstp) {
ic <- 0
betag <- beta1
}
beta0 <- beta1
}
}
apbetag <- getAp(xmat, center, betag, geemod)
S <- y - apbetag
vee <- veemat(S, center, structure, substructure, scores)
vc <- vee$vc
vee12 <- as.matrix(vee$v12)
vee12inv <- as.matrix(vee$v12inv)
Dmat <- getD(xmat, center, betag, geemod)
eitbg <- vee12inv %*% (y - apbetag)
Dmatb <- vee12inv %*% Dmat
wtstuff <- wtmat(Dmatb, eitbg, med, scores, hbrs)
wts <- wtstuff$wts
wts <- as.vector(wts)
Dmat <- getD(xmat, center, betag, geemod)
if (varcovst == "var2") {
vcmat <- varcov2(eitbg, center, scores, vee12inv, Dmat,
p, delta, hparm)
}
if (varcovst == "var1") {
vcmat <- varcov1(eitbg, center, scores, vee12inv, Dmat,
wts)
}
if (varcovst == "var3") {
vcmat <- varcov3(eitbg, scores, vee12inv, Dmat, p, delta,
hparm)
}
betahist <- colbeta[1:istp, ]
se <- sqrt(diag(vcmat))
tab <- cbind(betag, se, betag/se)
colnames(tab) <- c("Est", "SE", "t-ratio")
varcov <- vcmat
ehat1 <- y - xmat %*% betag
alphat2 <- median(ehat1) + alphat
ehat <- ehat1 - median(ehat1)
yhat <- y - ehat
list(tab = tab, betahist = betahist, varcov = varcov, stanhess = Dmatb,
istp = istp, vc = vc, beta0 = alphat2, residuals = ehat,
fitted.values = yhat)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.