R/cpgee_exc.R

In geeCRT: Bias-Corrected GEE for Cluster Randomized Trials

cpgee_exc <- function(
    y, X, id, m, family, maxiter, epsilon, printrange,
    alpadj) {
    ##################################################################################### MODULE: BEGINEND creates two vectors that have the start and
    ##################################################################################### end points for each cluster

    # INPUT n: vector of cluster sizes

    # OUTPUT first: vector with starting row for cluster i last:
    # vector with ending row for cluster i
    BEGINEND <- function(n) {
        last <- cumsum(n)
        first <- last - n + 1
        return(cbind(first, last))
    }

    ##################################################################################### Module: IS_POS_DEF A = symmetric matrix returns 1 if A is
    ##################################################################################### positive definite 0 otherwise
    is_pos_def <- function(A) {
        return(min(eigen(A)$values) > 1e-13)
    }

    ##################################################################################### Module: mu_v_c_fun X_c = covariates beta = parameters returns
    ##################################################################################### mean and unscaled variace
    mu_v_c_fun <- function(X_c, beta) {
        mu_c <- 1 / (1 + exp(c(-X_c %*% beta)))
        v_c <- mu_c * (1 - mu_c)
        return(list(mu_c = mu_c, v_c = v_c))
    }

    ##################################################################################### MODULE: GETROWB generates a row of the -E(d(cov)/dbeta) (ie Q)
    ##################################################################################### matrix

    # INPUT mu: Marginal means for cluster i j: Indicator for mean
    # k: Indicator for mean X: Covariate matrix for cluster i y:
    # Response vector for cluster i

    # OUTPUT row of E(d(cov)/dbeta) (ie Q) matrix
    GETROWB <- function(mu, j, k, X, y) {
        row <- -(y[j] - mu[j]) * X[j, ] * mu[j] * (1 - mu[j]) - (y[k] -
            mu[k]) * X[k, ] * mu[k] * (1 - mu[k])

        return(row)
    }

    ##################################################################################### MODULE: CREATEA creates residual for beta estimating equation,
    ##################################################################################### (Y - mu)

    # INPUT mu: vector of n_i marginal means y: outcome vector for
    # ith cluster

    # OUTPUT residuals for beta estimating equation
    CREATEA <- function(mu, y) {
        return(y - mu)
    }

    ##################################################################################### MODULE: CREATEC creates derivative matrix for beta estimating
    ##################################################################################### equation, dmu/dbeta

    # INPUT X: covariate matrix for cluster i mu: vector of n_i
    # marginal means

    # OUTPUT derivative matrix for beta estimating equation
    CREATEC <- function(X, mu) {
        D_mat <- X * (mu * (1 - mu))
        return(D_mat)
    }


    ##################################################################################### MODULE: CREATEB creates covariance matrix for beta estimating
    ##################################################################################### equation, var(Y)

    # INPUT gamma: vector of s_jk between mean j and mean k for
    # cluster i n: sample size (scalar) for cluster i

    # OUTPUT covariance matrix for beta estimating equation
    CREATEB <- function(gamma, n) {
        B <- matrix(0, n, n)
        l <- 1
        for (j in 1:n) {
            for (k in j:n) {
                B[j, k] <- gamma[l]
                l <- l + 1
            }
        }
        B[lower.tri(B)] <- t(B)[lower.tri(B)]
        return(B)
    }



    ##################################################################################### MODULE: CREATED creates derivative matrix for alpha estimating
    ##################################################################################### equation, deta/dalpha

    # INPUT mu: vector of marginal means v: vector of variances
    # alpha: correlation parameters n: sample size (scalar) for
    # cluster m: vector of cluster-period sizes

    # OUTPUT derivative matrix for alpha estimating equation
    CREATED <- function(mu, v, alpha, n, m) {
        alpha0 <- alpha[1]
        D <- NULL
        # fixed scale
        for (j in 1:(n - 1)) {
            dj <- (v[j] / m[j]) * (m[j] - 1)
            D <- rbind(D, c(dj))
            for (k in (j + 1):n) {
                djk0 <- sqrt(v[j] * v[k])
                D <- rbind(D, c(djk0))
            }
        }
        D <- rbind(D, c((v[n] / m[n]) * (m[n] - 1)))
        return(D)
    }


    ##################################################################################### MODULE: gammahat Creates vector of estimated covariances

    # INPUT mu: vector of marginal means alpha: correlation
    # estimates n: cluster size (# of periods) m: vector of
    # cluster-period sizes

    # OUTPUT vector of estimated covariances
    gammahat <- function(mu, v, alpha, n, m) {
        alpha0 <- alpha[1] # within-period ICC
        gamma_c <- NULL
        for (j in 1:(n - 1)) {
            sj <- (v[j] / m[j]) * (1 + (m[j] - 1) * alpha0)
            gamma_c <- c(gamma_c, sj)
            for (k in (j + 1):n) {
                sjk <- sqrt(v[j] * v[k]) * alpha0
                gamma_c <- c(gamma_c, sjk)
            }
        }
        gamma_c <- c(gamma_c, (v[n] / m[n]) * (1 + (m[n] - 1) * alpha0))
        return(gamma_c)
    }

    ##################################################################################### MODULE: checkalpha Creates compact correlation matrix for
    ##################################################################################### range checks

    # INPUT alpha: correlation estimates n: cluster size (# of
    # periods)

    # OUTPUT A correlation matrix indexed by cluster-periods
    checkalpha <- function(alpha, n) {
        alpha0 <- alpha[1] # within-period ICC
        L <- matrix(0, n, n)
        for (j in 1:(n - 1)) {
            L[j, j] <- alpha0
            for (k in (j + 1):n) {
                L[j, k] <- alpha0
            }
        }
        L[n, n] <- alpha0
        L[lower.tri(L)] <- t(L)[lower.tri(L)]
        return(L)
    }



    ##################################################################################### MODULE: SCORE generates the score matrix for each cluster and
    ##################################################################################### approximate information to be used to estimate parameters and
    ##################################################################################### generate standard errors

    # INPUT Ustarold: initial values for information matrix beta:
    # vector of marginal mean parameters alpha: vector of marginal
    # correlation parameters y: vector of outcomes (cluster-period
    # means) X: marginal mean covariates m: vector of cluster-period
    # sizes n: vector of cluster sample sizes p: number of marginal
    # mean parameters q: number of marginal correlation parameters
    # NPSDFLAG: checks if B is positive definite, terminates if not
    # NPSDADJFLAG: checks if the alpha adjustment factor matrix is
    # positive definite, terminates if not

    # OUTPUT U: score vector UUtran: sum of U_i*U_i` across all
    # clusters Ustar: approximate information matrix

    SCORE <- function(Ustarold, beta, alpha, y, X, m, n, p, q, NPSDFLAG,
                      NPSDADJFLAG) {
        U <- rep(0, p + q)
        UUtran <- Ustar <- matrix(0, p + q, p + q)
        naiveold <- ginv(Ustarold[1:p, 1:p]) # needed for Hi1 below

        locx <- BEGINEND(n)
        locz <- BEGINEND(choose(n + 1, 2))

        alpdata0 <- NULL
        alpdata1 <- NULL

        for (i in 1:length(n)) {
            X_c <- X[locx[i, 1]:locx[i, 2], , drop = FALSE]
            y_c <- y[locx[i, 1]:locx[i, 2]]
            m_c <- m[locx[i, 1]:locx[i, 2]]

            U_c <- rep(0, p + q)
            Ustar_c <- matrix(0, p + q, p + q)

            # mean and variance under different distribution families
            mu_v_c_res <- mu_v_c_fun(X_c, beta)
            mu_c <- mu_v_c_res$mu_c
            v_c <- mu_v_c_res$v_c

            # correlation matrix elements
            gamma_c <- gammahat(mu_c, v_c, alpha, n[i], m_c)

            R <- rep(0, choose(n[i] + 1, 2))
            DB <- matrix(0, choose(n[i] + 1, 2), p)

            C <- CREATEC(X_c, mu_c)
            B <- CREATEB(gamma_c, n[i])
            A <- CREATEA(mu_c, y_c)
            D <- CREATED(mu_c, v_c, alpha, n[i], m_c)

            INVB <- ginv(B)
            CtinvB <- t(C) %*% INVB
            Hi1 <- C %*% naiveold %*% CtinvB
            G_c <- tcrossprod(y_c - mu_c)

            if (alpadj) {
                CT <- t(C)
                omega <- C %*% naiveold %*% CT
                vminomega <- B - omega
                psd_vmin <- is_pos_def(vminomega)

                if (psd_vmin == 1) {
                    Ci <- B %*% ginv(vminomega)
                    G_c <- Ci %*% G_c
                } else {
                    NPSDADJFLAG <- 1
                    stop("(V - Omega) is not positive definite")
                }
            }


            l <- 1
            for (j in 1:n[i]) {
                for (k in j:n[i]) {
                    R[l] <- G_c[j, k] - gamma_c[l]
                    l <- l + 1
                }
            }

            # Check for positive definite of B
            if (min(eigen(B)$values) <= 0) {
                NPSDFLAG <- 1
                stop(paste(
                    "Var(Y) of Cluster", i, "is not Positive-Definite;",
                    "Joint Distribution Does Not Exist and Program terminates"
                ))
            }

            U_c[1:p] <- t(C) %*% INVB %*% A
            U_c[(p + 1):(p + q)] <- t(D) %*% R
            UUtran_c <- tcrossprod(U_c)
            Ustar_c[1:p, 1:p] <- t(C) %*% INVB %*% C
            Ustar_c[(p + 1):(p + q), (p + 1):(p + q)] <- crossprod(D)

            U <- U + U_c
            UUtran <- UUtran + UUtran_c
            Ustar <- Ustar + Ustar_c

            # build dataset for within-period ICC
            r1 <- cbind(diag(G_c), v_c, m_c)
            alpdata0 <- rbind(alpdata0, r1)

            # build dataset for remaining
            for (j in 1:(n[i] - 1)) {
                for (k in (j + 1):(n[i])) {
                    r2 <- c(G_c[j, k], v_c[j] * v_c[k], j, k, i)
                    alpdata1 <- rbind(alpdata1, r2)
                }
            }
        }
        return(list(
            U = U, UUtran = UUtran, Ustar = Ustar, alpdata0 = alpdata0,
            alpdata1 = alpdata1, NPSDFLAG = NPSDFLAG, NPSDADJFLAG = NPSDADJFLAG
        ))
    }

    ##################################################################################### MODULE: INITBETA generates initial values for beta
    ##################################################################################### approximates logistic regression using Newton's method

    # INPUT y: vector of cluster-period means m: vector of
    # cluster-period sizes (# of trials) X: marginal mean covariates

    # OUTPUT beta: vector of marginal mean parameters Ustar:
    # approximate score vector

    INITBETA <- function(y, m, X) {
        fit <- glm(cbind(m * y, m * (1 - y)) ~ -1 + X, family = binomial(link = "logit"))
        beta <- as.numeric(fit$coefficients)
        u <- as.numeric(fit$fitted.values)
        v <- u * (1 - u)
        Ustar <- t(X) %*% (X * v * m)

        return(list(beta = c(beta), Ustar = Ustar))
    }

    ##################################################################################### MODULE: getalpha (getalpha0) updates values for
    ##################################################################################### correlation parameters

    # INPUT alpdata0: dataset built for estimating alpha0 alpdata1:
    # dataset built for estimating the rest

    # OUTPUT alpha0: correltion parameter

    getalpha0 <- function(alpdata0, alpdata1) {
        s0 <- alpdata0[, 1]
        v <- alpdata0[, 2]
        m <- alpdata0[, 3]
        s1 <- alpdata1[, 1]
        vcross <- alpdata1[, 2]

        d1 <- sum(((m - 1) / m) * (s0 * v - v^2 / m))
        d0 <- sum(((m - 1) / m)^2 * v^2)
        e1 <- sum(sqrt(vcross) * s1)
        e0 <- sum(vcross)
        return((d1 + e1) / (d0 + e0))
    }

    ##################################################################################### MODULE: INITALPHA generates initial values for alpha

    # INPUT y: vector of cluster-period means m: vector of
    # cluster-period sizes (# of trials) X: marginal mean covariates
    # n: vector of cluster sizes (# of periods) beta: estimated
    # regression coefficients

    # OUTPUT alpha: estimated correlation value

    INITALPHA <- function(y, X, m, n, beta) {
        alpdata0 <- NULL
        alpdata1 <- NULL
        locx <- BEGINEND(n)

        for (i in 1:length(n)) {
            X_c <- X[locx[i, 1]:locx[i, 2], , drop = FALSE]
            y_c <- y[locx[i, 1]:locx[i, 2]]
            m_c <- m[locx[i, 1]:locx[i, 2]]

            mu_v_c_res <- mu_v_c_fun(X_c, beta)
            mu_c <- mu_v_c_res$mu_c
            v_c <- mu_v_c_res$v_c

            G_c <- tcrossprod(y_c - mu_c)

            # build dataset for within-period ICC
            r1 <- cbind(diag(G_c), v_c, m_c)
            alpdata0 <- rbind(alpdata0, r1)

            # build dataset for remaining
            for (j in 1:(n[i] - 1)) {
                for (k in (j + 1):(n[i])) {
                    r2 <- c(G_c[j, k], v_c[j] * v_c[k], j, k, i)
                    alpdata1 <- rbind(alpdata1, r2)
                }
            }
        }

        # produce initial estimates
        alpha0 <- getalpha0(alpdata0, alpdata1)
        return(alpha0)
    }


    ##################################################################################### MODULE: INVBIG compute (A - mm`)^{-1}c without performing the
    ##################################################################################### inverse directly

    # INPUT ainvc: inverse of matrix A times vector c ainvm: inverse
    # of matrix A times matrix (with low no. of columns) M M: matrix
    # of eigen column vectors m1,m2, ..  c: vector start: of do loop
    # end: of do loop, rank of X

    # OUTPUT ainvc: inverse of matrix A times vector c

    INVBIG <- function(ainvc, ainvm, m, c, start, end) {
        for (i in start:end) {
            b <- ainvm[, i]
            bt <- t(b)
            btm <- bt %*% m
            btmi <- btm[, i]
            gam <- 1 - btmi
            bg <- b / gam
            ainvc <- ainvc + bg %*% (bt %*% c)
            if (i < end) {
                ainvm <- ainvm + bg %*% btm
            }
        }
        return(ainvc)
    }

    ##################################################################################### MODULE: MAKEVAR creates covariance matrix of beta and alpha.

    # INPUT beta: vector of marginal mean parameters alpha: vector
    # of marginal correlation parameters Ustarold: initial values
    # for information matrix y: vector of cluster-period means X:
    # marginal mean covariates m: vector of cluster-period sizes n:
    # vector of cluster sample sizes p: number of marginal mean
    # parameters q: number of marginal correlation parameters
    # NPSDFLAG: checks if B is positive definite, terminates if not
    # NPSDADJFLAG: checks if the alpha adjustment factor matrix is
    # positive definite, terminates if not

    # OUTPUT robust: robust covariance matrix for beta and alpha
    # naive: naive (model-based) covariance matrix for beta varMD:
    # bias-corrected variance by Mancl and Derouen (2001) varKC:
    # bias-corrected variance by Kauermann and Carroll (2001) varFG:
    # bias-corrected variance by Fay and Graubard (2001)
    MAKEVAR <- function(Ustarold, beta, alpha, y, X, m, n, p, q,
                        ROBFLAG, NPSDFLAG, NPSDADJFLAG) {
        SCORE_RES <- SCORE(
            Ustarold, beta, alpha, y, X, m, n, p,
            q, NPSDFLAG, NPSDADJFLAG
        )
        U <- SCORE_RES$U
        UUtran <- SCORE_RES$UUtran
        Ustar <- SCORE_RES$Ustar
        NPSDFLAG <- SCORE_RES$NPSDFLAG
        NPSDADJFLAG <- SCORE_RES$NPSDADJFLAG

        naive <- ginv(Ustar[1:p, 1:p])
        naivealp <- ginv(Ustar[(p + 1):(p + q), (p + 1):(p + q)])

        # new commands to compute INV(I - H1)
        eigenRES1 <- eigen(naive)
        evals1 <- eigenRES1$values
        evecs1 <- eigenRES1$vectors
        sqrevals1 <- sqrt(evals1)
        sqe1 <- evecs1 %*% diag(sqrevals1)

        # new commands to compute INV(I - H2)
        eigenRES2 <- eigen(naivealp)
        evals2 <- eigenRES2$values
        evecs2 <- eigenRES2$vectors
        sqrevals2 <- sqrt(evals2)


        sqe2 <- evecs2 %*% diag(sqrevals2, 1)

        # Bias-corrected variance
        Ustar_c_array <- UUtran_c_array <- array(0, c(
            p + q, p + q,
            length(n)
        ))
        UUtran <- UUbc <- UUbc2 <- UUbc3 <- Ustar <- inustar <- matrix(
            0,
            p + q, p + q
        )

        locx <- BEGINEND(n)
        locz <- BEGINEND(choose(n + 1, 2))

        for (i in 1:length(n)) {
            X_c <- X[locx[i, 1]:locx[i, 2], , drop = FALSE]
            y_c <- y[locx[i, 1]:locx[i, 2]]
            m_c <- m[locx[i, 1]:locx[i, 2]]


            # mean and variance under different distribution families
            mu_v_c_res <- mu_v_c_fun(X_c, beta)
            mu_c <- mu_v_c_res$mu_c
            v_c <- mu_v_c_res$v_c

            U_i <- U_c <- rep(0, p + q)
            Ustar_c <- matrix(0, p + q, p + q)
            gamma_c <- gammahat(mu_c, v_c, alpha, n[i], m_c)

            # commands for beta
            C <- CREATEC(X_c, mu_c)
            B <- CREATEB(gamma_c, n[i])
            A <- CREATEA(mu_c, y_c)
            D <- CREATED(mu_c, v_c, alpha, n[i], m_c)

            INVB <- ginv(B)
            U_i[1:p] <- t(C) %*% INVB %*% A

            CtinvB <- t(C) %*% INVB
            Hi1 <- C %*% naive %*% CtinvB

            # commands for beta
            ai1 <- INVB
            mm1 <- C %*% sqe1
            ai1A <- ai1 %*% A
            ai1m1 <- ai1 %*% mm1
            ai1A <- INVBIG(ai1A, ai1m1, mm1, A, 1, p)
            U_c[1:p] <- t(C) %*% ai1A

            # commands for alpha
            R <- rep(0, choose(n[i] + 1, 2))
            DB <- matrix(0, choose(n[i] + 1, 2), p)
            G_c <- tcrossprod(y_c - mu_c)

            if (alpadj) {
                # MAEE
                CT <- t(C)
                omega <- C %*% naive %*% CT
                vminomega <- B - omega
                psd_vmin <- is_pos_def(vminomega)

                if (psd_vmin == 1) {
                    Ci <- B %*% ginv(vminomega)
                    G_c <- Ci %*% G_c
                } else {
                    NPSDADJFLAG <- 1
                    stop("(V - Omega) is not positive definite")
                }
            }

            # RANGE CHECKS
            rangeflag <- 0
            L_c <- checkalpha(alpha, n[i])
            l <- 1
            for (j in 1:n[i]) {
                for (k in j:n[i]) {
                    if ((L_c[j, k] >= min(sqrt((mu_c[j] * (1 - mu_c[k])) / (mu_c[k] *
                        (1 - mu_c[j]))), sqrt((mu_c[k] * (1 - mu_c[j])) / (mu_c[j] *
                        (1 - mu_c[k]))))) | (L_c[j, k] <= max(-sqrt((mu_c[j] *
                        mu_c[k]) / ((1 - mu_c[j]) * (1 - mu_c[k]))), -sqrt(((1 -
                        mu_c[j]) * (1 - mu_c[k])) / (mu_c[j] * mu_c[k]))))) {
                        rangeflag <- 1
                        if (printrange) {
                            warning(cat(
                                "Range Violation Detected for Cluster",
                                i, "and Pair", j, k, "\n"
                            ))
                        }
                        break
                    }

                    DB[l, ] <- GETROWB(mu_c, j, k, X_c, y_c)

                    R[l] <- G_c[j, k] - gamma_c[l]
                    l <- l + 1
                }
            }
            # Check for positive definite of B
            if (min(eigen(B)$values) <= 0) {
                NPSDFLAG <- 1
                stop(paste(
                    "Var(Y) of Cluster", i, "is not Positive-Definite;",
                    "Joint Distribution Does Not Exist and Program terminates"
                ))
            }

            U_i[(p + 1):(p + q)] <- t(D) %*% R
            mm2 <- D %*% sqe2
            ai2R <- R
            ai2m2 <- mm2
            ai2R <- INVBIG(ai2R, ai2m2, mm2, R, 1, q)
            U_c[(p + 1):(p + q)] <- t(D) %*% ai2R

            Ustar_c[1:p, 1:p] <- t(C) %*% INVB %*% C
            Ustar_c[(p + 1):(p + q), 1:p] <- t(D) %*% DB
            Ustar_c[(p + 1):(p + q), (p + 1):(p + q)] <- crossprod(D)
            Ustar <- Ustar + Ustar_c

            UUtran_c <- tcrossprod(U_i)
            UUtran <- UUtran + UUtran_c
            UUbc_c <- tcrossprod(U_c)
            UUbc <- UUbc + UUbc_c
            UUbc_ic <- tcrossprod(U_c, U_i)
            UUbc2 <- UUbc2 + UUbc_ic

            Ustar_c_array[, , i] <- Ustar_c
            UUtran_c_array[, , i] <- UUtran_c
        }

        inustar[1:p, 1:p] <- ginv(Ustar[1:p, 1:p])
        inustar[(p + 1):(p + q), (p + 1):(p + q)] <- ginv(Ustar[(p +
            1):(p + q), (p + 1):(p + q)])
        inustar[(p + 1):(p + q), 1:p] <- inustar[
            (p + 1):(p + q),
            (p + 1):(p + q)
        ] %*% Ustar[(p + 1):(p + q), 1:p] %*%
            inustar[1:p, 1:p]

        # the minus sign above is crucial, especially for large
        # correlation;
        inustartr <- t(inustar)

        # calculating adjustment factor for BC3
        for (i in 1:length(n)) {
            Hi <- diag(1 / sqrt(1 - pmin(0.75, c(diag(Ustar_c_array[
                , ,
                i
            ] %*% inustar)))))
            UUbc3 <- UUbc3 + Hi %*% UUtran_c_array[, , i] %*% Hi
        }

        # BC0 or usual Sandwich estimator
        robust <- inustar %*% UUtran %*% inustartr

        # BC1 or Variance estimator that extends Kauermann and Carroll
        # (2001);
        varKC <- inustar %*% (UUbc2 + t(UUbc2)) %*% inustartr / 2

        # BC2 or Variance estimator that extends Mancl and DeRouen
        # (2001);
        varMD <- inustar %*% UUbc %*% inustartr

        # BC3 or Variance estimator that extends Fay and Graubard
        # (2001);
        varFG <- inustar %*% UUbc3 %*% inustartr

        # model-based variance
        naive <- inustar[1:p, 1:p]

        if (min(diag(robust)) <= 0) {
            ROBFLAG <- 1
        }
        if (min(diag(varMD)) <= 0) {
            ROBFLAG <- 1
        }
        if (min(diag(varKC)) <= 0) {
            ROBFLAG <- 1
        }
        if (min(diag(varFG)) <= 0) {
            ROBFLAG <- 1
        }

        return(list(
            robust = robust, naive = naive, varMD = varMD,
            varKC = varKC, varFG = varFG, rangeflag = rangeflag,
            ROBFLAG = ROBFLAG, NPSDFLAG = NPSDFLAG, NPSDADJFLAG = NPSDADJFLAG
        ))
    }

    ##################################################################################### MODULE: FITPRENTICE Performs the paired estimating equations
    ##################################################################################### method

    # INPUT y: vector of cluster-period means X: marginal mean
    # covariates n: Vector of cluster sizes m: vector of
    # cluster-period sizes maxiter: maximum number of iterations
    # epsilon: tolerence for convergence SINGFLAG: THE ALGORITHM
    # terminated due to singular MB covariance matrix ROBFLAG: THE
    # ALGORITHM terminated due to singular robust variance matrix
    # NPSDFLAG: checks if B is positive definite, terminates if not
    # NPSDADJFLAG: checks if the alpha adjustment factor matrix is
    # positive definite, terminates if not

    # OUTPUT beta: a p x 1 vector of marginal mean parameter
    # estimates alpha: a q x 1 vector of marginal correlation
    # parameter estimates robust: robust covariance matrix for beta
    # and alpha naive: naive (model-based) covariance matrix for
    # beta varMD: bias-corrected variance due to Mancl and DeRouen
    # (2001) varKC: bias-corrected variance due to Kauermann and
    # Carroll (2001) varFG: bias-corrected variance due to Fay and
    # Graubard (2001) niter: number of iterations required for
    # convergence converge: did the algorithm converge (0 = no, 1 =
    # yes)
    FITPRENTICE <- function(y, X, m, n, maxiter, epsilon, SINGFLAG,
                            ROBFLAG, NPSDFLAG, NPSDADJFLAG) {
        p <- ncol(X)
        converge <- 0
        rangeflag <- 0


        INITRES <- INITBETA(y, m, X)
        beta <- INITRES$beta
        Ustar <- INITRES$Ustar

        alpha <- INITALPHA(y, X, m, n, beta)
        q <- length(alpha)
        # tolerances for beta and alpha
        delta <- rep(2 * epsilon, p)
        deltaalp <- rep(2 * epsilon, q)

        # iterative training
        niter <- 1
        while ((niter <= maxiter) & (max(abs(c(delta, deltaalp))) >
            epsilon)) {
            SINGFLAG <- 0
            NPSDFLAG <- 0
            NPSDADJFLAG <- 0
            Ustarold <- Ustar
            alphaold <- alpha

            # scores
            SCORE_RES <- SCORE(
                Ustarold, beta, alpha, y, X, m, n,
                p, q, NPSDFLAG, NPSDADJFLAG
            )
            U <- SCORE_RES$U
            UUtran <- SCORE_RES$UUtran
            Ustar <- SCORE_RES$Ustar
            NPSDFLAG <- SCORE_RES$NPSDFLAG
            NPSDADJFLAG <- SCORE_RES$NPSDADJFLAG

            # update alpha
            alpdata0 <- SCORE_RES$alpdata0
            alpdata1 <- SCORE_RES$alpdata1
            alpha <- getalpha0(alpdata0, alpdata1)
            deltaalp <- alpha - alphaold


            # update beta
            psdustar <- is_pos_def(Ustar[1:p, 1:p])
            mineig <- min(eigen(Ustar[1:p, 1:p])$values)
            if (psdustar) {
                delta <- solve(Ustar[1:p, 1:p], U[1:p])
                beta <- beta + delta
                converge <- (max(abs(c(delta, deltaalp))) <= epsilon)
            } else {
                SINGFLAG <- 1
            }
            niter <- niter + 1
        }

        Ustarold <- Ustar

        # inference
        MAKEVAR_RES <- MAKEVAR(
            Ustarold, beta, alpha, y, X, m, n,
            p, q, ROBFLAG, NPSDFLAG, NPSDADJFLAG
        )
        robust <- MAKEVAR_RES$robust
        naive <- MAKEVAR_RES$naive
        varMD <- MAKEVAR_RES$varMD
        varKC <- MAKEVAR_RES$varKC
        varFG <- MAKEVAR_RES$varFG
        rangeflag <- MAKEVAR_RES$rangeflag
        ROBFLAG <- MAKEVAR_RES$ROBFLAG
        NPSDFLAG <- MAKEVAR_RES$NPSDFLAG
        NPSDADJFLAG <- MAKEVAR_RES$NPSDADJFLAG
        return(list(
            beta = beta, alpha = alpha, robust = robust,
            naive = naive, varMD = varMD, varKC = varKC, varFG = varFG,
            niter = niter, converge = converge, SINGFLAG = SINGFLAG,
            ROBFLAG = ROBFLAG, NPSDFLAG = NPSDFLAG, NPSDADJFLAG = NPSDADJFLAG
        ))
    }

    ##################################################################################### MODULE: RESULTS creates printed output to screen of parameters
    ##################################################################################### and other information

    # INPUT beta: vector of marginal mean parameters alpha: vector
    # of marginal correlation Parameters robust: robust covariance
    # matrix for beta and alpha naive: naive (model-based)
    # covariance matrix for beta niter: number of iterations until
    # convergence n: vector of cluster sample sizes

    # OUTPUT to screen
    RESULTS <- function(beta, alpha, robust, naive, varMD, varKC,
                        varFG, niter, n) {
        p <- length(beta)
        q <- length(alpha)
        K <- length(n)

        # ** the next message is structure specific **
        corstr <- "Exchangeable"



        beta_numbers <- as.matrix(seq(1:p)) - 1
        bSE <- sqrt(diag(naive))
        bSEBC0 <- sqrt(diag(robust[1:p, 1:p]))
        bSEBC1 <- sqrt(diag(varKC[1:p, 1:p]))
        bSEBC2 <- sqrt(diag(varMD[1:p, 1:p]))
        bSEBC3 <- sqrt(diag(varFG[1:p, 1:p]))

        alpha_numbers <- as.matrix(seq(1:q)) - 1
        aSEBC0 <- sqrt(robust[(p + 1):(p + q), (p + 1):(p + q)])
        aSEBC1 <- sqrt(varKC[(p + 1):(p + q), (p + 1):(p + q)])
        aSEBC2 <- sqrt(varMD[(p + 1):(p + q), (p + 1):(p + q)])
        aSEBC3 <- sqrt(varFG[(p + 1):(p + q), (p + 1):(p + q)])



        outbeta <- cbind(
            beta_numbers, beta, bSE, bSEBC0, bSEBC1,
            bSEBC2, bSEBC3
        )
        outalpha <- cbind(
            alpha_numbers, alpha, aSEBC0, aSEBC1, aSEBC2,
            aSEBC3
        )
        colnames(outbeta) <- c(
            "Beta", "Estimate", "MB-stderr", "BC0-stderr",
            "BC1-stderr", "BC2-stderr", "BC3-stderr"
        )
        colnames(outalpha) <- c(
            "Alpha", "Estimate", "BC0-stderr",
            "BC1-stderr", "BC2-stderr", "BC3-stderr"
        )


        return(list(outbeta = outbeta, outalpha = outalpha))
    }

    # reasons for non-results are identified and tallied
    SINGFLAG <- 0
    ROBFLAG <- 0
    NPSDFLAG <- 0
    NPSDADJFLAG <- 0

    id1 <- id[order(id)]
    y <- y[order(id)]
    X <- X[order(id), ]
    m <- m[order(id)]
    id <- id1

    # Fit the GEE/MAEE algorithm
    n <- as.vector(table(id))

    PRENTICE_RES <- FITPRENTICE(
        y, X, m, n, maxiter, epsilon, SINGFLAG,
        ROBFLAG, NPSDFLAG, NPSDADJFLAG
    )
    beta <- PRENTICE_RES$beta
    alpha <- PRENTICE_RES$alpha

    robust <- PRENTICE_RES$robust
    naive <- PRENTICE_RES$naive
    varMD <- PRENTICE_RES$varMD
    varKC <- PRENTICE_RES$varKC
    varFG <- PRENTICE_RES$varFG
    niter <- PRENTICE_RES$niter

    converge <- PRENTICE_RES$converge
    SINGFLAG <- PRENTICE_RES$SINGFLAG
    ROBFLAG <- PRENTICE_RES$ROBFLAG
    NPSDFLAG <- PRENTICE_RES$NPSDFLAG
    NPSDADJFLAG <- PRENTICE_RES$NPSDADJFLAG

    # Final Results
    if (SINGFLAG == 1) {
        stop("Derivative matrix for beta is singular during updates")
    }
    if (ROBFLAG == 1) {
        stop("Sandwich variance is not positive definite")
    }
    if (converge == 0 & SINGFLAG == 0) {
        stop("The algorithm did not converge")
    }
    if (converge == 1 & ROBFLAG == 0) {
        result <- RESULTS(
            beta, alpha, robust, naive, varMD, varKC,
            varFG, niter, n
        )
        outList <- list(
            outbeta = result$outbeta, outalpha = result$outalpha,
            beta = beta, alpha = alpha, MB = naive, BC0 = robust,
            BC1 = varKC, BC2 = varMD, BC3 = varFG, niter = niter
        )
        class(outList) <- "cpgeeSWD"
        return(outList)
    }
}