R/svydfbetas.R

In svydiags: Linear Regression Model Diagnostics for Survey Data

svydfbetas <- function(mobj, stvar=NULL, clvar=NULL, z=3){
    xdat <- mobj$data[, all.vars(mobj$formula)]
    X <- model.matrix(mobj$formula, xdat)
    w <- mobj$weights
    n <- length(w)
    p <- length(mobj$coefficients)
    e <- mobj$residuals
    cutoff <- z / sqrt(n)

    XWX <- (t(X) * w) %*% X
    Ainv <- ginv(XWX)
    h <- svyhat(mobj)
    ewh <- e * w / (1-h)
    ewh <- matrix(rep(ewh,p), ncol=p, byrow=FALSE)
    Dfbeta <- Ainv %*% t(X * ewh)

    SEbeta <- sqrt(diag(mobj$cov.unscaled))
    SEbeta <- matrix(rep(SEbeta,n), ncol=n, byrow=FALSE)
    Dfbetas <- Dfbeta / SEbeta

    if (!is.null(clvar)){
        tmp <- svystdres(mobj, stvar, clvar)
        cutoff <- z / sqrt(tmp$n * (1 + (tmp$mbar-1)*tmp$rhohat))
    }

    DFBETAS <- list(Dfbeta = Dfbeta,
                    Dfbetas = Dfbetas,
                    cutoff = cutoff)
    DFBETAS
}