inst/doc/free1way.R

In free1way.docreg: Additional Documentation and Regression Tests for `stats::free1way'

### R code from vignette source 'free1way.Rnw'

###################################################
### code chunk number 1: citation
###################################################
yr <- format(dt <- as.Date(packageDescription("free1way.docreg")$Date), "%Y")
vs <- packageDescription("free1way.docreg")$Version
title <- "Semiparametrically Efficient Population and Permutation Inference in 
       Distribution-free Stratified $K$-sample Oneway Layouts"
DOI <- paste0("10.32614/CRAN.package.", packageDescription("free1way.docreg")$Package)


###################################################
### code chunk number 2: localfun
###################################################
Nsim <- 100
options(digits = 5)
.table2list <- function(x) 
{

    
dx <- dim(x)
if (length(dx) == 1L)
    stop("incorrect dimensions")
if (length(dx) == 2L)
    x <- as.table(array(x, dim = c(dx, 1)))
dx <- dim(x)
if (length(dx) < 3L)
    stop("incorrect dimensions")
C <- dim(x)[1L]
K <- dim(x)[2L]
B <- dim(x)[3L]
if (C < 2L)
    stop("at least two response categories required")
if (K < 2L)
    stop("at least two groups required")
xrc <- NULL
if (length(dx) == 4L) {
    if (dx[4] == 2L) {
        xrc <- array(x[,,,"FALSE", drop = TRUE], dim = dx[1:3])
        x <- array(x[,,,"TRUE", drop = TRUE], dim = dx[1:3])
    } else {
        stop(gettextf("%s currently only allows independent right-censoring",
                      "free1way"),
             domain = NA)
    }
}


xlist <- xrclist <- vector(mode = "list", length = B)

for (b in seq_len(B)) {
    xb <- matrix(x[,,b, drop = TRUE], ncol = K)
    xw <- rowSums(abs(xb)) > 0
    if (sum(xw) > 1L) {
        ### do not remove last parameter if there are corresponding
        ### right-censored observations
        wm <- which(xw)[sum(xw)]
        if (!is.null(xrc) && any(xrc[wm:dx[1],,b,drop = TRUE] > 0))
            xw[length(xw)] <- TRUE
        xlist[[b]] <- xb[xw,,drop = FALSE]
        Cidx <- rep.int(1L, times = C)
        Cidx[xw] <- Cidx[xw] + seq_len(sum(xw))
        attr(xlist[[b]], "idx") <- Cidx
        if (!is.null(xrc)) {
            ### count right-censored observations between distinct event
            ### times
            cs <- apply(xrc[,,b,drop = TRUE] * (!xw), 2, function(x) 
                diff(c(0, cumsum(x)[xw])))
            xrclist[[b]] <- matrix(xrc[xw,,b,drop = TRUE], ncol = K) + cs
            idx <- seq_len(C)[xw]
            idx <- rep(seq_len(sum(xw)), times = c(idx[1], diff(idx)))
            Cidx <- rep.int(1L, times = C)
            Cidx[seq_along(idx)] <- Cidx[seq_along(idx)] + idx
            attr(xrclist[[b]], "idx") <- Cidx
        }
    }
}
### remove empty blocks
strata <- !vapply(xlist, is.null, NA)
xlist <- xlist[strata]
xrclist <- xrclist[strata]



    ret <- list(xlist = xlist)
    if (!is.null(xrc))
        ret$xrclist <- xrclist
    ret$strata <- strata
    ret
}

.nll <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    bidx <- seq_len(ncol(x) - 1L)
    delta <- c(0, mu + parm[bidx])
    intercepts <- c(-Inf, parm[- bidx], Inf)
    tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                      ncol = ncol(x),
                                      byrow = TRUE)
    Ftmb <- F(tmb)
    if (rightcensored) {
        prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
    } else {
        prb <- Ftmb[- 1L, , drop = FALSE] - 
               Ftmb[- nrow(Ftmb), , drop = FALSE]
    } 
    
    if (any(prb < .Machine$double.eps^10)) 
        return(Inf)
    return(- sum(x * log(prb)))
}


.nsc <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    bidx <- seq_len(ncol(x) - 1L)
    delta <- c(0, mu + parm[bidx])
    intercepts <- c(-Inf, parm[- bidx], Inf)
    tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                      ncol = ncol(x),
                                      byrow = TRUE)
    Ftmb <- F(tmb)
    if (rightcensored) {
        prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
    } else {
        prb <- Ftmb[- 1L, , drop = FALSE] - 
               Ftmb[- nrow(Ftmb), , drop = FALSE]
    } 
    

    
    ftmb <- f(tmb)
    zu <- x * ftmb[- 1, , drop = FALSE] / prb
    if (rightcensored) zu[] <- 0 ### derivative of a constant
    zl <- x * ftmb[- nrow(ftmb), , drop = FALSE] / prb
    

    ret <- numeric(length(parm))
    ret[bidx] <- .colSums(zl, m = nrow(zl), n = ncol(zl))[-1L] -
                 .colSums(zu[-nrow(zu),,drop = FALSE], 
                          m = nrow(zu) - 1L, n = ncol(zu))[-1L]
    ret[- bidx] <- .rowSums(zu[-nrow(zu),,drop = FALSE] - 
                            zl[-1,,drop = FALSE], 
                            m = nrow(zu) - 1L, n = ncol(zu))
    return(- ret)
}


.nsr <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    bidx <- seq_len(ncol(x) - 1L)
    delta <- c(0, mu + parm[bidx])
    intercepts <- c(-Inf, parm[- bidx], Inf)
    tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                      ncol = ncol(x),
                                      byrow = TRUE)
    Ftmb <- F(tmb)
    if (rightcensored) {
        prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
    } else {
        prb <- Ftmb[- 1L, , drop = FALSE] - 
               Ftmb[- nrow(Ftmb), , drop = FALSE]
    } 
    

    
    ftmb <- f(tmb)
    zu <- x * ftmb[- 1, , drop = FALSE] / prb
    if (rightcensored) zu[] <- 0 ### derivative of a constant
    zl <- x * ftmb[- nrow(ftmb), , drop = FALSE] / prb
    

    ret <- .rowSums(zl - zu, m = nrow(zl), n = ncol(zl)) / 
           .rowSums(x, m = nrow(x), n = ncol(x))
    ret[!is.finite(ret)] <- 0
    return(- ret)
}


.hes <- function(parm, x, mu = 0, rightcensored = FALSE, full = FALSE) 
{

    
    bidx <- seq_len(ncol(x) - 1L)
    delta <- c(0, mu + parm[bidx])
    intercepts <- c(-Inf, parm[- bidx], Inf)
    tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                      ncol = ncol(x),
                                      byrow = TRUE)
    Ftmb <- F(tmb)
    if (rightcensored) {
        prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
    } else {
        prb <- Ftmb[- 1L, , drop = FALSE] - 
               Ftmb[- nrow(Ftmb), , drop = FALSE]
    } 
    

    
    ftmb <- f(tmb)
    fptmb <- fp(tmb)

    dl <- ftmb[- nrow(ftmb), , drop = FALSE]
    du <- ftmb[- 1, , drop = FALSE]
    if (rightcensored) du[] <- 0
    dpl <- fptmb[- nrow(ftmb), , drop = FALSE]
    dpu <- fptmb[- 1, , drop = FALSE]
    if (rightcensored) dpu[] <- 0
    dlm1 <- dl[,-1L, drop = FALSE]
    dum1 <- du[,-1L, drop = FALSE]
    dplm1 <- dpl[,-1L, drop = FALSE]
    dpum1 <- dpu[,-1L, drop = FALSE]
    prbm1 <- prb[,-1L, drop = FALSE]

    i1 <- length(intercepts) - 1L
    i2 <- 1L
    

    
    Aoffdiag <- - .rowSums(x * du * dl / prb^2, m = nrow(x), n = ncol(x))[-i2]
    Aoffdiag <- Aoffdiag[-length(Aoffdiag)]
    
    
    Adiag <- - .rowSums((x * dpu / prb)[-i1,,drop = FALSE] - 
                        (x * dpl / prb)[-i2,,drop = FALSE] - 
                        ((x * du^2 / prb^2)[-i1,,drop = FALSE] + 
                         (x * dl^2 / prb^2)[-i2,,drop = FALSE] ), 
                        m = nrow(x) - length(i1), n = ncol(x)
                       )
                      
    
    
    xm1 <- x[,-1L,drop = FALSE] 
    X <- ((xm1 * dpum1 / prbm1)[-i1,,drop = FALSE] - 
          (xm1 * dplm1 / prbm1)[-i2,,drop = FALSE] - 
          ((xm1 * dum1^2 / prbm1^2)[-i1,,drop = FALSE] - 
           (xm1 * dum1 * dlm1 / prbm1^2)[-i2,,drop = FALSE] -
           (xm1 * dum1 * dlm1 / prbm1^2)[-i1,,drop = FALSE] +
           (xm1 * dlm1^2 / prbm1^2)[-i2,,drop = FALSE]
          )
         )

    Z <- - .colSums(xm1 * (dpum1 / prbm1 - 
                           dplm1 / prbm1 -
                           (dum1^2 / prbm1^2 - 
                            2 * dum1 * dlm1 / prbm1^2 +
                            dlm1^2 / prbm1^2
                           )
                          ),
                    m = nrow(xm1), n = ncol(xm1)
                    )
    if (length(Z) > 1L) Z <- diag(Z)
    

    if (length(Adiag) > 1L) {
        if (!isFALSE(full)) {
            A <- list(Adiag = Adiag, Aoffdiag = Aoffdiag)
        } else {
            A <- Matrix::bandSparse(length(Adiag), 
                k = 0:1, diagonals = list(Adiag, Aoffdiag), 
                symmetric = TRUE)
        }
    } else {
        if (!isFALSE(full)) {
            A <- list(Adiag = Adiag, Aoffdiag = NULL)
        } else {
            A <- matrix(Adiag)
        }
    }
    return(list(A = A, X = X, Z = Z))
}


.snll <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    C <- vapply(x, NROW, 0L) ### might differ by stratum
    K <- unique(do.call("c", lapply(x, ncol))) ### the same
    B <- length(x)
    sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                   levels = seq_len(B))
    bidx <- seq_len(K - 1L)
    delta <- parm[bidx]
    intercepts <- split(parm[-bidx], sidx)
    

    ret <- 0
    for (b in seq_len(B))
        ret <- ret + .nll(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                          rightcensored = rightcensored)
    return(ret)
}


.snsc <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    C <- vapply(x, NROW, 0L) ### might differ by stratum
    K <- unique(do.call("c", lapply(x, ncol))) ### the same
    B <- length(x)
    sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                   levels = seq_len(B))
    bidx <- seq_len(K - 1L)
    delta <- parm[bidx]
    intercepts <- split(parm[-bidx], sidx)
    

    ret <- numeric(length(bidx))
    for (b in seq_len(B)) {
        nsc <- .nsc(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                    rightcensored = rightcensored)
        ret[bidx] <- ret[bidx] + nsc[bidx]
        ret <- c(ret, nsc[-bidx])
    }
    return(ret)
}


.shes <- function(parm, x, mu = 0, xrc = NULL, full = FALSE, 
                  retMatrix = FALSE) 
{

    
    C <- vapply(x, NROW, 0L) ### might differ by stratum
    K <- unique(do.call("c", lapply(x, ncol))) ### the same
    B <- length(x)
    sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                   levels = seq_len(B))
    bidx <- seq_len(K - 1L)
    delta <- parm[bidx]
    intercepts <- split(parm[-bidx], sidx)
    

    if (!isFALSE(ret <- full)) {
        
        for (b in seq_len(B)) {
            H <- .hes(c(delta, intercepts[[b]]), x[[b]], mu = mu, full = full)
            if (!is.null(xrc)) {
                Hrc <- .hes(c(delta, intercepts[[b]]), xrc[[b]], mu = mu, 
                            rightcensored = TRUE, full = full)
                H$X <- H$X + Hrc$X
                H$A$Adiag <- H$A$Adiag + Hrc$A$Adiag
                H$A$Aoffdiag <- H$A$Aoffdiag + Hrc$A$Aoffdiag
                H$Z <- H$Z + Hrc$Z
            }
            if (b == 1L) {
                Adiag <- H$A$Adiag
                Aoffdiag <- H$A$Aoffdiag
                X <- H$X
                Z <- H$Z
            } else {
                Adiag <- c(Adiag, H$A$Adiag)
                Aoffdiag <- c(Aoffdiag, 0, H$A$Aoffdiag)
                X <- rbind(X, H$X)
                Z <- Z + H$Z
            }
        }

        if (length(Adiag) > 1L) {
            A <- Matrix::bandSparse(length(Adiag),
                                    k = 0:1, diagonals = list(Adiag, Aoffdiag),
                                    symmetric = TRUE)
        } else {
            A <- matrix(Adiag)
        }

        ret <- cbind(Z, t(X))
        ret <- rbind(ret, cbind(X, A))
        if (retMatrix) return(ret)
        return(as.matrix(ret))
        
    }
    ret <- matrix(0, nrow = length(bidx), ncol = length(bidx))
    for (b in seq_len(B)) {
        H <- .hes(c(delta, intercepts[[b]]), x[[b]], mu = mu)
        if (!is.null(xrc)) {
            Hrc <- .hes(c(delta, intercepts[[b]]), xrc[[b]], mu = mu, 
                        rightcensored = TRUE)
            H$X <- H$X + Hrc$X
            H$A <- H$A + Hrc$A
            H$Z <- H$Z + Hrc$Z
        }
        sAH <- tryCatch(Matrix::solve(H$A, H$X), error = function(e) NULL)
        if (is.null(sAH))
            stop(gettextf("error computing the Hessian in %s",
                          "free1way"),
                 domain = NA)
        ret <- ret + (H$Z - crossprod(H$X, sAH))
    }
    as.matrix(ret)
}


.snsr <- function(parm, x, mu = 0, rightcensored = FALSE) 
{

    
    C <- vapply(x, NROW, 0L) ### might differ by stratum
    K <- unique(do.call("c", lapply(x, ncol))) ### the same
    B <- length(x)
    sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                   levels = seq_len(B))
    bidx <- seq_len(K - 1L)
    delta <- parm[bidx]
    intercepts <- split(parm[-bidx], sidx)
    

    ret <- c()
    for (b in seq_len(B)) {
        idx <- attr(x[[b]], "idx")
        ### idx == 1L means zero residual, see definition of idx
        sr <- c(0, .nsr(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                        rightcensored = rightcensored))
        ret <- c(ret, sr[idx])
    }
    return(ret)
}


.free1wayML <- function(x, link, mu = 0, start = NULL, fix = NULL, 
                        residuals = TRUE, score = TRUE, hessian = TRUE, 
                        MPL_Jeffreys = FALSE,
                        ### use nlminb for small sample sizes
                        dooptim = c(".NewtonRaphson", "nlminb")[1 + (sum(x) < 20)],                         
                        control = list(
                            "nlminb" = list(trace = trace, iter.max = 200,
                                            eval.max = 200, rel.tol = 1e-10,
                                            abs.tol = 1e-20, xf.tol = 1e-16),
                            ".NewtonRaphson" = list(iter.max = 200, trace = trace, 
                                             objtol = 5e-4, 
                                             gradtol = 1e-5 * sum(x) / 1000, 
                                             paramtol = 1e-5, minstepsize = 1e-2, 
                                             tolsolve = .Machine$double.eps)
                        )[dooptim],
                        trace = FALSE, 
                        tol = sqrt(.Machine$double.eps), ...) 
{

    ### convert to three-way table
    xt <- x
    if (!is.table(x))
        stop(gettextf("invalid argument '%s'", "x"), domain = NA) # 'y' in free1way ...
    dx <- dim(x)
    dn <- dimnames(x)
    if (length(dx) == 2L) {
        x <- as.table(array(c(x), dim = dx <- c(dx, 1L)))
        dimnames(x) <- dn <- c(dn, list(A = "A"))
    }

    ### short-cuts for link functions
    F <- function(q) .p(link, q = q)
    Q <- function(p) .q(link, p = p)
    f <- function(q) .d(link, x = q)
    fp <- function(q) .dd(link, x = q)

    if(!suppressPackageStartupMessages(requireNamespace("Matrix")))
        stop(gettextf("%s needs package 'Matrix' correctly installed",
                      ".free1wayML"),
                 domain = NA)

    

    dx <- dim(x)
    if (length(dx) == 1L)
        stop("incorrect dimensions")
    if (length(dx) == 2L)
        x <- as.table(array(x, dim = c(dx, 1)))
    dx <- dim(x)
    if (length(dx) < 3L)
        stop("incorrect dimensions")
    C <- dim(x)[1L]
    K <- dim(x)[2L]
    B <- dim(x)[3L]
    if (C < 2L)
        stop("at least two response categories required")
    if (K < 2L)
        stop("at least two groups required")
    xrc <- NULL
    if (length(dx) == 4L) {
        if (dx[4] == 2L) {
            xrc <- array(x[,,,"FALSE", drop = TRUE], dim = dx[1:3])
            x <- array(x[,,,"TRUE", drop = TRUE], dim = dx[1:3])
        } else {
            stop(gettextf("%s currently only allows independent right-censoring",
                          "free1way"),
                 domain = NA)
        }
    }


    xlist <- xrclist <- vector(mode = "list", length = B)

    for (b in seq_len(B)) {
        xb <- matrix(x[,,b, drop = TRUE], ncol = K)
        xw <- rowSums(abs(xb)) > 0
        if (sum(xw) > 1L) {
            ### do not remove last parameter if there are corresponding
            ### right-censored observations
            wm <- which(xw)[sum(xw)]
            if (!is.null(xrc) && any(xrc[wm:dx[1],,b,drop = TRUE] > 0))
                xw[length(xw)] <- TRUE
            xlist[[b]] <- xb[xw,,drop = FALSE]
            Cidx <- rep.int(1L, times = C)
            Cidx[xw] <- Cidx[xw] + seq_len(sum(xw))
            attr(xlist[[b]], "idx") <- Cidx
            if (!is.null(xrc)) {
                ### count right-censored observations between distinct event
                ### times
                cs <- apply(xrc[,,b,drop = TRUE] * (!xw), 2, function(x) 
                    diff(c(0, cumsum(x)[xw])))
                xrclist[[b]] <- matrix(xrc[xw,,b,drop = TRUE], ncol = K) + cs
                idx <- seq_len(C)[xw]
                idx <- rep(seq_len(sum(xw)), times = c(idx[1], diff(idx)))
                Cidx <- rep.int(1L, times = C)
                Cidx[seq_along(idx)] <- Cidx[seq_along(idx)] + idx
                attr(xrclist[[b]], "idx") <- Cidx
            }
        }
    }
    ### remove empty blocks
    strata <- !vapply(xlist, is.null, NA)
    xlist <- xlist[strata]
    xrclist <- xrclist[strata]
    
    
    ## allow specification of start = delta and fix = 1:K
    ## for evaluating the likelihood at given delta parameters
    ## without having to specify all intercept parameters
    if (is.null(start))
        start <- rep.int(0, K - 1L)
    NS <- length(start) == (K - 1L)
    lwr <- rep(-Inf, times = K - 1L)
    for (b in seq_len(length(xlist))) {
        bC <- nrow(xlist[[b]]) - 1L
        lwr <- c(lwr, -Inf, rep.int(0, times = bC - 1L))
        if (NS) {
            ecdf0 <- cumsum(rowSums(xlist[[b]]))
            ### ensure that 0 < ecdf0 < 1 such that quantiles exist
            ecdf0 <- pmax(1, ecdf0[-length(ecdf0)]) / (max(ecdf0) + 1)
            Qecdf <- Q(ecdf0)
            start <- c(start, Qecdf)
        }
    }
    
    
    .nll <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        bidx <- seq_len(ncol(x) - 1L)
        delta <- c(0, mu + parm[bidx])
        intercepts <- c(-Inf, parm[- bidx], Inf)
        tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                          ncol = ncol(x),
                                          byrow = TRUE)
        Ftmb <- F(tmb)
        if (rightcensored) {
            prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
        } else {
            prb <- Ftmb[- 1L, , drop = FALSE] - 
                   Ftmb[- nrow(Ftmb), , drop = FALSE]
        } 
        
        if (any(prb < .Machine$double.eps^10)) 
            return(Inf)
        return(- sum(x * log(prb)))
    }
    
    
    .nsc <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        bidx <- seq_len(ncol(x) - 1L)
        delta <- c(0, mu + parm[bidx])
        intercepts <- c(-Inf, parm[- bidx], Inf)
        tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                          ncol = ncol(x),
                                          byrow = TRUE)
        Ftmb <- F(tmb)
        if (rightcensored) {
            prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
        } else {
            prb <- Ftmb[- 1L, , drop = FALSE] - 
                   Ftmb[- nrow(Ftmb), , drop = FALSE]
        } 
        

        
        ftmb <- f(tmb)
        zu <- x * ftmb[- 1, , drop = FALSE] / prb
        if (rightcensored) zu[] <- 0 ### derivative of a constant
        zl <- x * ftmb[- nrow(ftmb), , drop = FALSE] / prb
        

        ret <- numeric(length(parm))
        ret[bidx] <- .colSums(zl, m = nrow(zl), n = ncol(zl))[-1L] -
                     .colSums(zu[-nrow(zu),,drop = FALSE], 
                              m = nrow(zu) - 1L, n = ncol(zu))[-1L]
        ret[- bidx] <- .rowSums(zu[-nrow(zu),,drop = FALSE] - 
                                zl[-1,,drop = FALSE], 
                                m = nrow(zu) - 1L, n = ncol(zu))
        return(- ret)
    }
    
    
    .nsr <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        bidx <- seq_len(ncol(x) - 1L)
        delta <- c(0, mu + parm[bidx])
        intercepts <- c(-Inf, parm[- bidx], Inf)
        tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                          ncol = ncol(x),
                                          byrow = TRUE)
        Ftmb <- F(tmb)
        if (rightcensored) {
            prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
        } else {
            prb <- Ftmb[- 1L, , drop = FALSE] - 
                   Ftmb[- nrow(Ftmb), , drop = FALSE]
        } 
        

        
        ftmb <- f(tmb)
        zu <- x * ftmb[- 1, , drop = FALSE] / prb
        if (rightcensored) zu[] <- 0 ### derivative of a constant
        zl <- x * ftmb[- nrow(ftmb), , drop = FALSE] / prb
        

        ret <- .rowSums(zl - zu, m = nrow(zl), n = ncol(zl)) / 
               .rowSums(x, m = nrow(x), n = ncol(x))
        ret[!is.finite(ret)] <- 0
        return(- ret)
    }
    
    
    .hes <- function(parm, x, mu = 0, rightcensored = FALSE, full = FALSE) 
    {

        
        bidx <- seq_len(ncol(x) - 1L)
        delta <- c(0, mu + parm[bidx])
        intercepts <- c(-Inf, parm[- bidx], Inf)
        tmb <- intercepts - matrix(delta, nrow = length(intercepts),  
                                          ncol = ncol(x),
                                          byrow = TRUE)
        Ftmb <- F(tmb)
        if (rightcensored) {
            prb <- 1 - Ftmb[- nrow(Ftmb), , drop = FALSE]
        } else {
            prb <- Ftmb[- 1L, , drop = FALSE] - 
                   Ftmb[- nrow(Ftmb), , drop = FALSE]
        } 
        

        
        ftmb <- f(tmb)
        fptmb <- fp(tmb)

        dl <- ftmb[- nrow(ftmb), , drop = FALSE]
        du <- ftmb[- 1, , drop = FALSE]
        if (rightcensored) du[] <- 0
        dpl <- fptmb[- nrow(ftmb), , drop = FALSE]
        dpu <- fptmb[- 1, , drop = FALSE]
        if (rightcensored) dpu[] <- 0
        dlm1 <- dl[,-1L, drop = FALSE]
        dum1 <- du[,-1L, drop = FALSE]
        dplm1 <- dpl[,-1L, drop = FALSE]
        dpum1 <- dpu[,-1L, drop = FALSE]
        prbm1 <- prb[,-1L, drop = FALSE]

        i1 <- length(intercepts) - 1L
        i2 <- 1L
        

        
        Aoffdiag <- - .rowSums(x * du * dl / prb^2, m = nrow(x), n = ncol(x))[-i2]
        Aoffdiag <- Aoffdiag[-length(Aoffdiag)]
        
        
        Adiag <- - .rowSums((x * dpu / prb)[-i1,,drop = FALSE] - 
                            (x * dpl / prb)[-i2,,drop = FALSE] - 
                            ((x * du^2 / prb^2)[-i1,,drop = FALSE] + 
                             (x * dl^2 / prb^2)[-i2,,drop = FALSE] ), 
                            m = nrow(x) - length(i1), n = ncol(x)
                           )
                          
        
        
        xm1 <- x[,-1L,drop = FALSE] 
        X <- ((xm1 * dpum1 / prbm1)[-i1,,drop = FALSE] - 
              (xm1 * dplm1 / prbm1)[-i2,,drop = FALSE] - 
              ((xm1 * dum1^2 / prbm1^2)[-i1,,drop = FALSE] - 
               (xm1 * dum1 * dlm1 / prbm1^2)[-i2,,drop = FALSE] -
               (xm1 * dum1 * dlm1 / prbm1^2)[-i1,,drop = FALSE] +
               (xm1 * dlm1^2 / prbm1^2)[-i2,,drop = FALSE]
              )
             )

        Z <- - .colSums(xm1 * (dpum1 / prbm1 - 
                               dplm1 / prbm1 -
                               (dum1^2 / prbm1^2 - 
                                2 * dum1 * dlm1 / prbm1^2 +
                                dlm1^2 / prbm1^2
                               )
                              ),
                        m = nrow(xm1), n = ncol(xm1)
                        )
        if (length(Z) > 1L) Z <- diag(Z)
        

        if (length(Adiag) > 1L) {
            if (!isFALSE(full)) {
                A <- list(Adiag = Adiag, Aoffdiag = Aoffdiag)
            } else {
                A <- Matrix::bandSparse(length(Adiag), 
                    k = 0:1, diagonals = list(Adiag, Aoffdiag), 
                    symmetric = TRUE)
            }
        } else {
            if (!isFALSE(full)) {
                A <- list(Adiag = Adiag, Aoffdiag = NULL)
            } else {
                A <- matrix(Adiag)
            }
        }
        return(list(A = A, X = X, Z = Z))
    }
    
    
    .snll <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        C <- vapply(x, NROW, 0L) ### might differ by stratum
        K <- unique(do.call("c", lapply(x, ncol))) ### the same
        B <- length(x)
        sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                       levels = seq_len(B))
        bidx <- seq_len(K - 1L)
        delta <- parm[bidx]
        intercepts <- split(parm[-bidx], sidx)
        

        ret <- 0
        for (b in seq_len(B))
            ret <- ret + .nll(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                              rightcensored = rightcensored)
        return(ret)
    }
    
    
    .snsc <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        C <- vapply(x, NROW, 0L) ### might differ by stratum
        K <- unique(do.call("c", lapply(x, ncol))) ### the same
        B <- length(x)
        sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                       levels = seq_len(B))
        bidx <- seq_len(K - 1L)
        delta <- parm[bidx]
        intercepts <- split(parm[-bidx], sidx)
        

        ret <- numeric(length(bidx))
        for (b in seq_len(B)) {
            nsc <- .nsc(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                        rightcensored = rightcensored)
            ret[bidx] <- ret[bidx] + nsc[bidx]
            ret <- c(ret, nsc[-bidx])
        }
        return(ret)
    }
    
    
    .shes <- function(parm, x, mu = 0, xrc = NULL, full = FALSE, 
                      retMatrix = FALSE) 
    {

        
        C <- vapply(x, NROW, 0L) ### might differ by stratum
        K <- unique(do.call("c", lapply(x, ncol))) ### the same
        B <- length(x)
        sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                       levels = seq_len(B))
        bidx <- seq_len(K - 1L)
        delta <- parm[bidx]
        intercepts <- split(parm[-bidx], sidx)
        

        if (!isFALSE(ret <- full)) {
            
            for (b in seq_len(B)) {
                H <- .hes(c(delta, intercepts[[b]]), x[[b]], mu = mu, full = full)
                if (!is.null(xrc)) {
                    Hrc <- .hes(c(delta, intercepts[[b]]), xrc[[b]], mu = mu, 
                                rightcensored = TRUE, full = full)
                    H$X <- H$X + Hrc$X
                    H$A$Adiag <- H$A$Adiag + Hrc$A$Adiag
                    H$A$Aoffdiag <- H$A$Aoffdiag + Hrc$A$Aoffdiag
                    H$Z <- H$Z + Hrc$Z
                }
                if (b == 1L) {
                    Adiag <- H$A$Adiag
                    Aoffdiag <- H$A$Aoffdiag
                    X <- H$X
                    Z <- H$Z
                } else {
                    Adiag <- c(Adiag, H$A$Adiag)
                    Aoffdiag <- c(Aoffdiag, 0, H$A$Aoffdiag)
                    X <- rbind(X, H$X)
                    Z <- Z + H$Z
                }
            }

            if (length(Adiag) > 1L) {
                A <- Matrix::bandSparse(length(Adiag),
                                        k = 0:1, diagonals = list(Adiag, Aoffdiag),
                                        symmetric = TRUE)
            } else {
                A <- matrix(Adiag)
            }

            ret <- cbind(Z, t(X))
            ret <- rbind(ret, cbind(X, A))
            if (retMatrix) return(ret)
            return(as.matrix(ret))
            
        }
        ret <- matrix(0, nrow = length(bidx), ncol = length(bidx))
        for (b in seq_len(B)) {
            H <- .hes(c(delta, intercepts[[b]]), x[[b]], mu = mu)
            if (!is.null(xrc)) {
                Hrc <- .hes(c(delta, intercepts[[b]]), xrc[[b]], mu = mu, 
                            rightcensored = TRUE)
                H$X <- H$X + Hrc$X
                H$A <- H$A + Hrc$A
                H$Z <- H$Z + Hrc$Z
            }
            sAH <- tryCatch(Matrix::solve(H$A, H$X), error = function(e) NULL)
            if (is.null(sAH))
                stop(gettextf("error computing the Hessian in %s",
                              "free1way"),
                     domain = NA)
            ret <- ret + (H$Z - crossprod(H$X, sAH))
        }
        as.matrix(ret)
    }
    
    
    .snsr <- function(parm, x, mu = 0, rightcensored = FALSE) 
    {

        
        C <- vapply(x, NROW, 0L) ### might differ by stratum
        K <- unique(do.call("c", lapply(x, ncol))) ### the same
        B <- length(x)
        sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                       levels = seq_len(B))
        bidx <- seq_len(K - 1L)
        delta <- parm[bidx]
        intercepts <- split(parm[-bidx], sidx)
        

        ret <- c()
        for (b in seq_len(B)) {
            idx <- attr(x[[b]], "idx")
            ### idx == 1L means zero residual, see definition of idx
            sr <- c(0, .nsr(c(delta, intercepts[[b]]), x[[b]], mu = mu,
                            rightcensored = rightcensored))
            ret <- c(ret, sr[idx])
        }
        return(ret)
    }
    
    

    fn <- function(par) 
    {
        ret <- .snll(par, x = xlist, mu = mu)
        if (!is.null(xrc))
            ret <- ret + .snll(par, x = xrclist, mu = mu, 
                               rightcensored = TRUE)
        return(ret)
    }
    gr <- function(par) 
    {
        ret <- .snsc(par, x = xlist, mu = mu)
        if (!is.null(xrc))
            ret <- ret + .snsc(par, x = xrclist, mu = mu, 
                               rightcensored = TRUE)
        return(ret)
    }

    ### allocate memory for hessian
    Hess <- Matrix::Matrix(0, nrow = length(start), ncol = length(start))

    he <- function(par) 
    {
        if (!is.null(xrc)) {
            ret <- .shes(par, x = xlist, mu = mu, xrc = xrclist, full = Hess, 
                         retMatrix = names(control)[1L] == ".NewtonRaphson")
        } else {
            ret <- .shes(par, x = xlist, mu = mu, full = Hess, 
                         retMatrix = names(control)[1L] == ".NewtonRaphson")
        }
        return(ret)
    }
    

    .profile <- function(start, fix = seq_len(K - 1)) 
    {
        if (!all(fix %in% seq_len(K - 1)))
            stop(gettextf("invalid argument '%s'", "fix"), domain = NA)
        delta <- start[fix]
        opargs <- list(start = start[-fix], 
                         objective = function(par) {
                             p <- numeric(length(par) + length(fix))
                             p[fix] <- delta
                             p[-fix] <- par
                             fn(p)
                         },
                         gradient = function(par) {
                             p <- numeric(length(par) + length(fix))
                             p[fix] <- delta
                             p[-fix] <- par
                             gr(p)[-fix]
                         },
                         hessian = function(par) {
                             p <- numeric(length(par) + length(fix))
                             p[fix] <- delta
                             p[-fix] <- par
                             he(p)[-fix, -fix, drop = FALSE]
                         })
        opargs$control <- control[[1L]]
        MPL_Jeffreys <- FALSE ### turn off Jeffreys penalisation in .profile

        
        maxit <- control[[1L]]$iter.max
        while(maxit < 10001) {
           ret <- do.call(names(control)[[1L]], opargs)
           maxit <- 5 * maxit
           if (ret$convergence > 0) {
               opargs$control$eval.max <- maxit
               opargs$control$iter.max <- maxit
               opargs$start <- ret$par
           } else {
               break()
           }
        }

        if (isTRUE(MPL_Jeffreys)) {
            
            .pll_Jeffreys <- function(cf, start) 
            {
                fix <- seq_along(cf)
                start[fix] <- cf
                ### compute profile likelihood w/o warnings
                ret <- suppressWarnings(.profile(start, fix = fix))
                Hfull <- he(ret$par)
                Hfix <- as.matrix(solve(solve(Hfull)[fix, fix]))
                return(ret$value - 
                       .5 * determinant(Hfix, logarithm = TRUE)$modulus)
            }
            if (K == 2) {
                MLcf <- ret$par[seq_len(K - 1)]
                Fret <- optim(MLcf, fn = .pll_Jeffreys, start = ret$par,
                              method = "Brent", lower = MLcf - 5, 
                              upper = MLcf + 5)
            } else {
                ### Nelder-Mead
                Fret <- optim(ret$par[seq_len(K - 1)], fn = .pll_Jeffreys, 
                              start = ret$par)
            }
            if (Fret$convergence == 0) {
                start <- ret$par
                start[seq_len(K - 1)] <- Fret$par
                ret <- .profile(start, fix = seq_len(K - 1))
                ret$objective <- ret$value
            }
            
        } else {
            if (ret$convergence > 0) {
                if (is.na(MPL_Jeffreys)) { ### only after failure
                    warning(gettextf("Jeffreys penalisation was applied in %s because initial optimisation failed with:",
                                     "free1way"),
                            "\n  ", ret$message, domain = NA)
                    MPL_Jeffreys <- TRUE
                    
                    .pll_Jeffreys <- function(cf, start) 
                    {
                        fix <- seq_along(cf)
                        start[fix] <- cf
                        ### compute profile likelihood w/o warnings
                        ret <- suppressWarnings(.profile(start, fix = fix))
                        Hfull <- he(ret$par)
                        Hfix <- as.matrix(solve(solve(Hfull)[fix, fix]))
                        return(ret$value - 
                               .5 * determinant(Hfix, logarithm = TRUE)$modulus)
                    }
                    if (K == 2) {
                        MLcf <- ret$par[seq_len(K - 1)]
                        Fret <- optim(MLcf, fn = .pll_Jeffreys, start = ret$par,
                                      method = "Brent", lower = MLcf - 5, 
                                      upper = MLcf + 5)
                    } else {
                        ### Nelder-Mead
                        Fret <- optim(ret$par[seq_len(K - 1)], fn = .pll_Jeffreys, 
                                      start = ret$par)
                    }
                    if (Fret$convergence == 0) {
                        start <- ret$par
                        start[seq_len(K - 1)] <- Fret$par
                        ret <- .profile(start, fix = seq_len(K - 1))
                        ret$objective <- ret$value
                    }
                    
                }
           }
        }
        if (ret$convergence > 0)
            warning(gettextf("unsuccessful optimisation in %s", "free1way"),
                    ": ", ret$message, domain = NA)

        ret$MPL_Jeffreys <- MPL_Jeffreys
        ret$value <- ret$objective
        ret$objective <- NULL
        

        p <- numeric(length(start))
        p[fix] <- delta
        p[-fix] <- ret$par
        ret$par <- p
        ret
    }
    
    
    if (!length(fix)) {
        opargs <- list(start = start, 
                       objective = fn, 
                       gradient = gr,
                       hessian = he)
        opargs$control <- control[[1L]]
        
        maxit <- control[[1L]]$iter.max
        while(maxit < 10001) {
           ret <- do.call(names(control)[[1L]], opargs)
           maxit <- 5 * maxit
           if (ret$convergence > 0) {
               opargs$control$eval.max <- maxit
               opargs$control$iter.max <- maxit
               opargs$start <- ret$par
           } else {
               break()
           }
        }

        if (isTRUE(MPL_Jeffreys)) {
            
            .pll_Jeffreys <- function(cf, start) 
            {
                fix <- seq_along(cf)
                start[fix] <- cf
                ### compute profile likelihood w/o warnings
                ret <- suppressWarnings(.profile(start, fix = fix))
                Hfull <- he(ret$par)
                Hfix <- as.matrix(solve(solve(Hfull)[fix, fix]))
                return(ret$value - 
                       .5 * determinant(Hfix, logarithm = TRUE)$modulus)
            }
            if (K == 2) {
                MLcf <- ret$par[seq_len(K - 1)]
                Fret <- optim(MLcf, fn = .pll_Jeffreys, start = ret$par,
                              method = "Brent", lower = MLcf - 5, 
                              upper = MLcf + 5)
            } else {
                ### Nelder-Mead
                Fret <- optim(ret$par[seq_len(K - 1)], fn = .pll_Jeffreys, 
                              start = ret$par)
            }
            if (Fret$convergence == 0) {
                start <- ret$par
                start[seq_len(K - 1)] <- Fret$par
                ret <- .profile(start, fix = seq_len(K - 1))
                ret$objective <- ret$value
            }
            
        } else {
            if (ret$convergence > 0) {
                if (is.na(MPL_Jeffreys)) { ### only after failure
                    warning(gettextf("Jeffreys penalisation was applied in %s because initial optimisation failed with:",
                                     "free1way"),
                            "\n  ", ret$message, domain = NA)
                    MPL_Jeffreys <- TRUE
                    
                    .pll_Jeffreys <- function(cf, start) 
                    {
                        fix <- seq_along(cf)
                        start[fix] <- cf
                        ### compute profile likelihood w/o warnings
                        ret <- suppressWarnings(.profile(start, fix = fix))
                        Hfull <- he(ret$par)
                        Hfix <- as.matrix(solve(solve(Hfull)[fix, fix]))
                        return(ret$value - 
                               .5 * determinant(Hfix, logarithm = TRUE)$modulus)
                    }
                    if (K == 2) {
                        MLcf <- ret$par[seq_len(K - 1)]
                        Fret <- optim(MLcf, fn = .pll_Jeffreys, start = ret$par,
                                      method = "Brent", lower = MLcf - 5, 
                                      upper = MLcf + 5)
                    } else {
                        ### Nelder-Mead
                        Fret <- optim(ret$par[seq_len(K - 1)], fn = .pll_Jeffreys, 
                                      start = ret$par)
                    }
                    if (Fret$convergence == 0) {
                        start <- ret$par
                        start[seq_len(K - 1)] <- Fret$par
                        ret <- .profile(start, fix = seq_len(K - 1))
                        ret$objective <- ret$value
                    }
                    
                }
           }
        }
        if (ret$convergence > 0)
            warning(gettextf("unsuccessful optimisation in %s", "free1way"),
                    ": ", ret$message, domain = NA)

        ret$MPL_Jeffreys <- MPL_Jeffreys
        ret$value <- ret$objective
        ret$objective <- NULL
        
    } else if (length(fix) == length(start)) {
        ret <- list(par = start, 
                    value = fn(start))
    } else {
        ret <- .profile(start, fix = fix)
    }
     
    
    if (is.null(fix) || (length(fix) == length(start)))
        parm <- seq_len(K - 1)
    else 
        parm <- fix
    if (any(parm >= K)) return(ret)

    ret$coefficients <- ret$par[parm]
    dn2 <- dimnames(xt)[2L]
    names(ret$coefficients) <- cnames <- paste0(names(dn2), dn2[[1L]][1L + parm])

    par <- ret$par
    intercepts <- function(parm, x) 
    {

        
        C <- vapply(x, NROW, 0L) ### might differ by stratum
        K <- unique(do.call("c", lapply(x, ncol))) ### the same
        B <- length(x)
        sidx <- factor(rep(seq_len(B), times = pmax(0, C - 1L)), 
                       levels = seq_len(B))
        bidx <- seq_len(K - 1L)
        delta <- parm[bidx]
        intercepts <- split(parm[-bidx], sidx)
        

        return(intercepts)
    }
    ret$intercepts <- intercepts(par, x = xlist)

    if (score) {
        ret$negscore <- .snsc(par, x = xlist, mu = mu)[parm]
        if (!is.null(xrc))
            ret$negscore <- ret$negscore + .snsc(par, x = xrclist, mu = mu, 
                                                 rightcensored = TRUE)[parm]
    }
    if (hessian) {
        if (!is.null(xrc)) {
            ret$hessian <- .shes(par, x = xlist, mu = mu, xrc = xrclist)
        } else {
            ret$hessian <- .shes(par, x = xlist, mu = mu)
        }
        ret$vcov <- solve(ret$hessian)
        if (length(parm) != nrow(ret$hessian))
           ret$hessian <- solve(ret$vcov <- ret$vcov[parm, parm, drop = FALSE])
        rownames(ret$vcov) <- colnames(ret$vcov) <- rownames(ret$hessian) <-
            colnames(ret$hessian) <-  cnames
    }
    if (residuals) {
        ret$negresiduals <- .snsr(par, x = xlist, mu = mu)
        if (!is.null(xrc)) {
            rcr <- .snsr(par, x = xrclist, mu = mu, rightcensored = TRUE)
            ret$negresiduals <- c(rbind(matrix(ret$negresiduals, nrow = C),
                                        matrix(rcr, nrow = C)))
         }
    }
    ret$profile <- function(start, fix)
        .free1wayML(xt, link = link, mu = mu, start = start, fix = fix, tol = tol, 
                    ...) 
    ret$table <- xt

    ret$strata <- strata
    ret$mu <- mu
    if (length(ret$mu) == 1) {
        names(ret$mu) <- link$parm
    } else {
        names(ret$mu) <- c(paste(link$parm, cnames[1L], sep = ":"), cnames[-1L])
    }
    

    class(ret) <- "free1wayML"
    ret
}


.SW <- function(res, xt) 
{

    if (length(dim(xt)) == 3L) {
        res <- matrix(res, nrow = dim(xt)[1L], ncol = dim(xt)[3])
        STAT <-  Exp <- Cov <- 0
        for (b in seq_len(dim(xt)[3L])) {
            sw <- .SW(res[,b, drop = TRUE], xt[,,b, drop = TRUE])
            STAT <- STAT + sw$Statistic
            Exp <- Exp + sw$Expectation
            Cov <- Cov + sw$Covariance
        }
        return(list(Statistic = STAT, Expectation = as.vector(Exp),
                    Covariance = Cov))
    }

    Y <- matrix(res, ncol = 1, nrow = length(xt))
    weights <- c(xt)
    x <- gl(ncol(xt), nrow(xt))
    X <- model.matrix(~ x, data = data.frame(x = x))[,-1L,drop = FALSE]

    w. <- sum(weights)
    wX <- weights * X
    wY <- weights * Y
    ExpX <- colSums(wX)
    ExpY <- colSums(wY) / w.
    CovX <- crossprod(X, wX)
    Yc <- t(t(Y) - ExpY)
    CovY <- crossprod(Yc, weights * Yc) / w.
    Exp <- kronecker(ExpY, ExpX)
    Cov <- w. / (w. - 1) * kronecker(CovY, CovX) -
           1 / (w. - 1) * kronecker(CovY, tcrossprod(ExpX))
    STAT <- crossprod(X, wY)
    list(Statistic = STAT, Expectation = as.vector(Exp),
         Covariance = Cov)
}


.resample <- function(res, xt, B = 10000) 
{

    if (length(dim(xt)) == 2L)
        xt <- as.table(array(xt, dim = c(dim(xt), 1)))

    res <- matrix(res, nrow = dim(xt)[1L], ncol = dim(xt)[3L])
    stat <- 0
    ret <- .SW(res, xt)
    if (dim(xt)[2L] == 2L) {
        ret$testStat <- c((ret$Statistic - ret$Expectation) / 
                          sqrt(c(ret$Covariance)))
    } else {
        ES <- ret$Statistic - ret$Expectation
        ret$testStat <- sum(ES * solve(ret$Covariance, ES))
    }
    ret$DF <- dim(xt)[2L] - 1L

    if (B) {
        for (j in 1:dim(xt)[3L]) {
           rt <- r2dtable(B, r = rowSums(xt[,,j]), c = colSums(xt[,,j]))
           stat <- stat + vapply(rt, 
               function(x) .colSums(x[,-1L, drop = FALSE] * res[,j], 
                                    m = nrow(x), n = ncol(x) - 1L), 
                                 FUN.VALUE = rep(0, dim(xt)[[2L]] - 1L))
        }
        if (dim(xt)[2L] == 2L) {
             ret$permStat <- (stat - ret$Expectation) / 
                              sqrt(c(ret$Covariance))
        } else {
            ES <- matrix(stat, ncol = B) - ret$Expectation
            ret$permStat <- .colSums(ES * solve(ret$Covariance, ES), 
                                     m = dim(xt)[[2L]] - 1L, n = B)
        }
    }
    ret
}


### distribution function
.p <- function(link, q, ...)
    link$linkinv(q = q, ...)

### quantile function
.q <- function(link, p, ...)
    link$link(p = p, ...)

### density function
.d <- function(link, x, ...)
    link$dlinkinv(x = x, ...)

### derivative of density function
.dd <- function(link, x, ...)
    link$ddlinkinv(x = x, ...)

### 2nd derivative of density function
.ddd <- function(link, x, ...)
    link$dddlinkinv(x = x, ...)

### ratio of derivative of density to
### density function
.dd2d <- function(link, x, ...)
    link$dd2dlinkinv(x = x, ...)

### constructor
linkfun <- function(name,	### nickname
                    alias,	### char 
                    model, 	### char, semiparametric model name
                    parm, 	### char, parameter name
                    link,      	### quantile function
                    linkinv,   	### distribution function
                    dlinkinv,  	### density function
                    ddlinkinv, 	### derivative of density function
                    ...) 
{

    ret <- list(name = name, 
                alias = alias,
                model = model,
                parm = parm,
                link = link,
                linkinv = linkinv,
                dlinkinv = dlinkinv,
                ddlinkinv = ddlinkinv)
    if (is.null(ret$dd2d)) 
        ret$dd2d <- function(x) 
            ret$ddlinkinv(x) / ret$dlinkinv(x)
    ret <- c(ret, list(...))
    class(ret) <- "linkfun"
    ret
}


logit <- function()
    linkfun(name = "Logit", 
            alias = c("Wilcoxon", "Kruskal-Wallis"),
            model = "proportional odds", 
            parm = "log-odds ratio",
            link = qlogis,
            linkinv = plogis,
            dlinkinv = dlogis,
            ddlinkinv = function(x) {
                p <- plogis(x)
                p * (1 - p)^2 - p^2 * (1 - p)
            },
            dddlinkinv = function(x) {
                ex <- exp(x)
                ifelse(is.finite(x), (ex - 4 * ex^2 + ex^3) / (1 + ex)^4, 0.0)
            },
            dd2d = function(x) {
                ex <- exp(x)
                (1 - ex) / (1 + ex)
            },
            parm2PI = function(x) {
               OR <- exp(x)
               ret <- OR * (OR - 1 - x)/(OR - 1)^2
               ret[abs(x) < .Machine$double.eps] <- 0.5
               return(ret)
            },
            PI2parm = function(p) {
               f <- function(x, PI)
                   x + (exp(-x) * (PI + 
                                   exp(2 * x) * (PI - 1) + 
                                   exp(x) * (1 - 2 * PI)))
               ret <- vapply(p, function(p) 
                   uniroot(f, PI = p, interval = 50 * c(-1, 1))$root, 0)
               return(ret)
            },
            parm2OVL = function(x) 2 * plogis(-abs(x / 2))
    )


probit <- function()
    linkfun(name = "Probit",
            alias = "van der Waerden normal scores",
            model = "latent normal shift", 
            parm = "generalised Cohen's d",
            link = qnorm,
            linkinv = pnorm,
            dlinkinv = dnorm,
            ddlinkinv = function(x) 
                ifelse(is.finite(x), -dnorm(x = x) * x, 0.0), 
            dddlinkinv = function(x) 
                ifelse(is.finite(x), dnorm(x = x) * (x^2 - 1), 0.0),
            dd2d = function(x) -x,
            parm2PI = function(x) pnorm(x, sd = sqrt(2)),
            PI2parm = function(p) qnorm(p, sd = sqrt(2)),
            parm2OVL = function(x) 2 * pnorm(-abs(x / 2))
    )


cloglog <- function()
    linkfun(name = "Complementary Log-log",
            alias = "Savage",
            model = "proportional hazards", 
            parm = "log-hazard ratio",
            link = function(p, log.p = FALSE) {
                if (log.p) p <- exp(p)
                log(-log1p(- p))
            },
            linkinv = function(q, lower.tail = TRUE, log.p = FALSE) {
                ### p = 1 - exp(-exp(q))
                ret <- exp(-exp(q))
                if (log.p) {
                    if (lower.tail)
                        return(log1p(-ret))
                    return(-exp(q))
                }
                if (lower.tail)
                    return(-expm1(-exp(q)))
                return(ret)
            },
            dlinkinv = function(x) 
                ifelse(is.finite(x), exp(x - exp(x)), 0.0),
            ddlinkinv = function(x) {
                ex <- exp(x)
                ifelse(is.finite(x), (ex - ex^2) / exp(ex), 0.0)
            },
            dddlinkinv = function(x) {
                ex <- exp(x)
                ifelse(is.finite(x), (ex - 3*ex^2 + ex^3) / exp(ex), 0.0)
            },
            dd2d = function(x)
               -expm1(x),
            parm2PI = plogis,
            PI2parm = qlogis,
            parm2OVL = function(x) {
                x <- abs(x)
                ret <- exp(x / (exp(-x) - 1)) - exp(-x / (exp(x) - 1)) + 1 
                ret[abs(x) < .Machine$double.eps] <- 1
                x[] <- ret
                return(x)
            }
    )


loglog <- function()
    linkfun(name = "Log-log",
            alias = "Lehmann", 
            model = "Lehmann", 
            parm = "log-reverse time hazard ratio",
            link = function(p, log.p = FALSE) {
                if (!log.p) p <- log(p)
                -log(-p)
            },
            linkinv = function(q, lower.tail = TRUE, log.p = FALSE) {
                ### p = exp(-exp(-q))
                if (log.p) {
                    if (lower.tail)
                        return(-exp(-q))
                    return(log1p(-exp(-exp(-q))))
                }
                if (lower.tail)
                    return(exp(-exp(-q)))
                -expm1(-exp(-q))
            },
            dlinkinv = function(x) 
                ifelse(is.finite(x), exp(- x - exp(-x)), 0.0),
            ddlinkinv = function(x) {
               ex <- exp(-x)
               ifelse(is.finite(x), exp(-ex - x) * (ex - 1.0), 0.0)
            },
            dddlinkinv = function(x) {
               ex <- exp(-x)
               ifelse(is.finite(x), exp(-x - ex) * (ex - 1)^2 - 
                                    exp(-ex - 2 * x), 
                                    0.0)
            },
            dd2d = function(x) 
                expm1(-x),
            parm2PI = plogis,
            PI2parm = qlogis,
            parm2OVL = function(x) {
                x <- abs(x)
                rt <- exp(-x / (exp(x) - 1))
                ret <- rt^exp(x) + 1 - rt
                ret[abs(x) < .Machine$double.eps] <- 1
                x[] <- ret
                return(x)
            }
    )


### adopted from rms:::lrm.fit
.NewtonRaphson <- function(start, objective, gradient, hessian, 
                           control = list(iter.max = 150, trace = trace, 
                                          objtol = 5e-4, gradtol = 1e-5, 
                                          paramtol = 1e-5, minstepsize = 1e-2, 
                                          tolsolve = .Machine$double.eps),
                           trace = FALSE)
{

    theta  <- start # Initialize the parameter vector
    oldobj <- Inf
    objthe <- objective(theta)
    if (!is.finite(objthe)) {
        msg <- "Infeasible starting values"
        return(list(par = theta, objective = objthe, convergence = 1, 
                    message = msg)) 
    }

    if(!suppressPackageStartupMessages(requireNamespace("Matrix")))
        stop(gettextf("%s needs package 'Matrix' correctly installed",
                      ".NewtonRaphson"),
                 domain = NA)

    for (iter in seq_len(control$iter.max)) {

        
        gradthe <- gradient(theta)     # Compute the gradient vector
        hessthe <- hessian(theta)      # Compute the Hessian matrix

        delta <- Matrix::solve(hessthe, gradthe, tol = control$tolsolve)

        if (control$trace)
            cat(iter, ': ', theta, "\n", sep = "")

        step_size <- 1L                # Initialize step size for step-halving
        

        # Step-halving loop
        while (TRUE) {
            
            new_theta <- theta - step_size * delta # Update parameter vector
            objnew_the <- objective(new_theta)

            if (control$trace)
                cat("Old, new, old - new objective:", 
                    objthe, objnew_the, objthe - objnew_the, "\n")

            # Objective function failed to be reduced or is infinite
            if (!is.finite(objnew_the) || (objnew_the > objthe + 1e-6)) {
                step_size <- step_size / 2         # Reduce the step size

                if (control$trace) 
                    cat("Step size reduced to", step_size, "\n")

                if (step_size <= control$minstepsize) {
                    msg <- paste("Step size ", step_size, 
                                 " has reduced below minstepsize")
                    return(list(par = theta, objective = objthe, convergence = 1, 
                                message = msg)) 
                }
            } else {
                theta  <- new_theta	# accept the new parameter vector
                oldobj <- objthe
                objthe <- objnew_the
                break
            }
            
        }

        
        # Convergence check - must meet 3 criteria
        if ((objthe <= oldobj + 1e-6 && (oldobj - objthe < control$objtol)) &&
            (max(abs(gradthe)) < control$gradtol) &&
            (max(abs(delta)) < control$paramtol))

            return(list(par            = theta,
                        objective      = objthe,
                        convergence    = 0,
                        message        = "Normal convergence"))
         
    }

    msg <- paste("Reached", control$iter.max, "iterations without convergence")
    return(list(par = theta, objective = objthe, convergence = 1, message = msg)) 
}



###################################################
### code chunk number 3: glm
###################################################
library("free1way.docreg")
(x <- matrix(c(10, 5, 7, 11, 8, 9), nrow = 2))
d <- expand.grid(y = relevel(gl(2, 1), "2"), t = gl(3, 1))
d$x <- c(x)
m <- glm(y ~ t, data = d, weights = x, family = binomial())
(cf <- coef(m))


###################################################
### code chunk number 4: glm-op
###################################################
F <- plogis
f <- dlogis
op <- optim(par = c("mt2" = 0, "mt3" = 0, "(Intercept)" = 0), 
            fn = .nll, gr = .nsc, 
            x = x, method = "BFGS", hessian = TRUE)
cbind(glm = c(cf[-1] * -1, cf[1]), free1way = op$par)
logLik(m)
-op$value


###################################################
### code chunk number 5: glm-H
###################################################
fp <- function(x) 
{
    p <- plogis(x)
    p * (1 - p)^2 - p^2 * (1 - p)
}
H <- .hes(op$par, x)
### analytical covariance of parameters
solve(H$Z - crossprod(H$X, Matrix::solve(H$A, H$X)))
### numerical covariance
solve(op$hessian)[1:2,1:2]
### from glm
vcov(m)[-1,-1]


###################################################
### code chunk number 6: glm-free1way
###################################################
obj <- .free1wayML(as.table(x), link = logit())
obj$coefficients
-obj$value
### analytical covariance
obj$vcov


###################################################
### code chunk number 7: glm-stratum
###################################################
(x <- as.table(array(c(10, 5, 7, 11, 8, 9,
                        9, 4, 8, 15, 5, 4), dim = c(2, 3, 2))))
d <- expand.grid(y = relevel(gl(2, 1), "2"), t = gl(3, 1), s = gl(2, 1))
d$x <- c(x)
m <- glm(y ~ 0 + s + t, data = d, weights = x, family = binomial())
logLik(m)
(cf <- coef(m))


###################################################
### code chunk number 8: glm-op-stratum
###################################################
xl <- .table2list(x)$xlist
op <- optim(par = c("mt2" = 0, "mt3" = 0, "(Intercept 1)" = 0, "(Intercept 2)" = 0), 
            fn = .snll, gr = .snsc, 
            x = xl, 
            method = "BFGS", 
            hessian = TRUE)
cbind(glm = c(cf[-(1:2)] * -1, cf[1:2]), free1way = op$par)
logLik(m)
-op$value


###################################################
### code chunk number 9: glm-H-stratum
###################################################
### analytical covariance of parameters
solve(.shes(op$par, xl))
### numerical covariance
solve(op$hessian)[1:2,1:2]
### from glm
vcov(m)[-(1:2),-(1:2)]


###################################################
### code chunk number 10: glm-free1way-strata
###################################################
obj <- .free1wayML(as.table(x), link = logit())
obj$coefficients
-obj$value
### analytical covariance
obj$vcov


###################################################
### code chunk number 11: Newton
###################################################



###################################################
### code chunk number 12: Newton-test
###################################################
N <- 10000
P <- 30
X <- matrix(rnorm(N * P), ncol = P)
y <- X %*% runif(P) + rnorm(nrow(X))

f <- function(par) sum((y - X %*% par)^2)
g <- function(par) colSums(- 2 * c(y - X %*% par) * X)
h <- function(par) 2 * crossprod(X)

start <- runif(P)

cf <- .NewtonRaphson(start = start, objective = f, gradient = g, hessian = h)

cf2 <- coef(m <- lm(y ~ 0 + X))
all.equal(sum((y - fitted(m))^2), cf$objective)
all.equal(unname(cf$par), unname(cf2))


###################################################
### code chunk number 13: workhorse
###################################################
N <- 10
a <- matrix(c(5, 6, 4,
                    3, 5, 7,
                    3, 4, 5,
                    3, 5, 6,
                    0, 0, 0,
                    4, 6, 5), ncol = 3, byrow = TRUE)
x <- as.table(array(c(a[1:3,], a[-(1:3),]), dim = c(3, 3, 2)))
x
ret <- .free1wayML(x, logit())
ret[c("value", "par")]
cf <- ret$par
cf[1:2] <- cf[1:2] + .5
### new2old parameterisation
c(cf[1:2], cf[3], log(cf[4] - cf[3]), cf[5])
### profile for cf[1:2]
.free1wayML(x, logit(), start = cf, fix = 1:2)[c("value", "par")]
### profile for cf[2]
.free1wayML(x, logit(), start = cf, fix = 2)[c("value", "par")]
### evaluate log-likelihood at cf
.free1wayML(x, logit(), start = cf, 
            fix = seq_along(ret$par))[c("value", "par")]


###################################################
### code chunk number 14: SW
###################################################
set.seed(29)
w <- gl(2, 15)
(s <- .SW(r <- rank(u <- runif(length(w))), model.matrix(~ 0 + w)))
ps <- .resample(r, model.matrix(~ 0 + w), B = 100000)
ps$testStat^2
mean(abs(ps$permStat) > abs(ps$testStat) - .Machine$double.eps)
pchisq(ps$testStat^ifelse(ps$DF == 1, 2, 1), df = ps$DF, lower.tail = FALSE)
### exactly the same
kruskal.test(u ~ w)
library("coin")
### almost the same
kruskal_test(u ~ w, distribution = approximate(100000))


###################################################
### code chunk number 15: Wexact
###################################################
wilcox_test(u ~ w, distribution = "exact")
free1way(u ~ w, exact = TRUE)


###################################################
### code chunk number 16: Wexact-le
###################################################
wilcox_test(u ~ w, distribution = "exact", alternative = "less")
print(free1way(u ~ w, exact = TRUE), alternative = "greater")


###################################################
### code chunk number 17: Wexact-gr
###################################################
wilcox_test(u ~ w, distribution = "exact", alternative = "greater")
print(free1way(u ~ w, exact = TRUE), alternative = "less")


###################################################
### code chunk number 18: formula (eval = FALSE)
###################################################
## y ~ groups | blocks


###################################################
### code chunk number 19: free
###################################################
x
### asymptotic permutation test
(ft <- free1way(x))
coef(ft)
vcov(ft)
### Wald per parameter
summary(ft)
library("multcomp")
summary(glht(ft), test = univariate())

### global Wald
summary(ft, test = "Wald")
summary(glht(ft), test = Chisqtest())

### Rao score, Permutation score, LRT
summary(ft, test = "Rao")
summary(ft, test = "Permutation")
summary(ft, test = "LRT")

### Wald confidence intervals, unadjusted
confint(glht(ft), calpha = univariate_calpha())
confint(ft, test = "Wald")

### Rao and LRT intervals
confint(ft, test = "Rao")
confint(ft, test = "LRT")


###################################################
### code chunk number 20: wilcox-formula
###################################################
N <- 25
w <- gl(2, N)
y <- rlogis(length(w), location = c(0, 1)[w])

#### link = logit is default
ft <- free1way(y ~ w)

### Wald 
summary(ft)

### Permutation test
wilcox.test(y ~ w, alternative = "greater", correct = FALSE, exact = FALSE)$p.value
pvalue(wilcox_test(y ~ w, alternative = "greater"))
summary(ft, test = "Permutation", alternative = "less")$p.value
wilcox.test(y ~ w, alternative = "less", correct = FALSE, exact = FALSE)$p.value
pvalue(wilcox_test(y ~ w, alternative = "less"))
summary(ft, test = "Permutation", alternative = "greater")$p.value
wilcox.test(y ~ w, correct = FALSE, exact = FALSE)$p.value
kruskal.test(y ~ w)$p.value
pvalue(wilcox_test(y ~ w))
summary(ft, test = "Permutation")$p.value

### Wald tests
summary(ft, test = "Wald", alternative = "less")
summary(ft, test = "Wald", alternative = "greater")
summary(ft, test = "Wald")

### Rao score tests
summary(ft, test = "Rao", alternative = "less")
summary(ft, test = "Rao", alternative = "greater")
summary(ft, test = "Rao")

### LRT (only two-sided)
summary(ft, test = "LRT")

### confidence intervals for log-odds ratios
confint(ft, test = "Permutation")
confint(ft, test = "LRT")
confint(ft, test = "Wald")
confint(ft, test = "Rao")

### confidence interval for "Wilcoxon Parameter" = PI = AUC
confint(ft, test = "Rao", what = "AUC")

### comparison with rms::orm
library("rms")
rev(coef(or <- orm(y ~ w)))[1]
coef(ft)
logLik(or)
logLik(ft)
vcov(or)[2,2]
vcov(ft)
ci <- confint(or)
ci[nrow(ci),]
confint(ft, test = "Wald")


###################################################
### code chunk number 21: mh
###################################################
example(mantelhaen.test, echo = FALSE)
mantelhaen.test(UCBAdmissions, correct = FALSE)
ft <- free1way(UCBAdmissions)
summary(ft, test = "Wald")
exp(coef(ft))
exp(confint(ft, test = "Wald"))
exp(sapply(dimnames(UCBAdmissions)[[3L]], function(dept)
       confint(free1way(UCBAdmissions[,,dept]), test = "Permutation")))
sapply(dimnames(UCBAdmissions)[[3L]], function(dept)
       fisher.test(UCBAdmissions[,,dept], conf.int = TRUE)$conf.int)


###################################################
### code chunk number 22: pt
###################################################
prop.test(UCBAdmissions[,,1], correct = FALSE)
summary(free1way(UCBAdmissions[,,1]), test = "Rao")


###################################################
### code chunk number 23: kw
###################################################
example(kruskal.test, echo = FALSE)
kruskal.test(x ~ g)
free1way(x ~ g)


###################################################
### code chunk number 24: sw
###################################################
library("survival")
N <- 10
nd <- expand.grid(g = gl(3, N), s = gl(3, N))
nd$tm <- rexp(nrow(nd))
nd$ev <- TRUE
cm <- coxph(Surv(tm, ev) ~ g + strata(s), data = nd)

(ft <- free1way(tm ~ g | s, data = nd, link = "cloglog"))
coef(cm)
coef(ft)
vcov(cm)
vcov(ft)
### Rao score tests
summary(cm)$sctest
summary(ft, test = "Rao")
### likelihood ratio tests
summary(cm)$logtest
summary(ft, test = "LRT")
### Wald tests
summary(cm)$waldtest
summary(ft, test = "Wald")
### asymptotic permutation tests
survdiff(Surv(tm, ev) ~ g + strata(s), data = nd, rho = 0)[c("chisq", "pvalue")]
summary(ft, test = "Permutation")
library("coin")
independence_test(Surv(tm, ev) ~ g | s, data = nd, ytrafo = function(...)
                  trafo(..., numeric_trafo = logrank_trafo, block = nd$s), 
                  teststat = "quad")


###################################################
### code chunk number 25: Peto
###################################################
survdiff(Surv(tm, ev) ~ g + strata(s), data = nd, rho = 1)[c("chisq", "pvalue")]
(ft <- free1way(tm ~ g | s, data = nd, link = "logit"))
summary(ft)
summary(ft, test = "Rao")
summary(ft, test = "LRT")
summary(ft, test = "Wald")
summary(ft, test = "Permutation")


###################################################
### code chunk number 26: sw
###################################################
library("survival")
data("GBSG2", package = "TH.data")
cm <- coxph(Surv(time, cens) ~ horTh + strata(tgrade), data = GBSG2)

ft <- with(GBSG2, free1way(Surv(time, cens) ~ horTh | tgrade, 
                           link = "cloglog"))
coef(cm)
coef(ft)
vcov(cm)
vcov(ft)
### Rao score tests
summary(cm)$sctest
summary(ft, test = "Rao")
### likelihood ratio tests
summary(cm)$logtest
summary(ft, test = "LRT")
### Wald tests
summary(cm)$waldtest
summary(ft, test = "Wald")
### asymptotic permutation tests
survdiff(Surv(time, cens) ~ horTh + strata(tgrade), data = GBSG2, rho = 0)[c("chisq", "pvalue")]
summary(ft, test = "Permutation")
independence_test(Surv(time, cens) ~ horTh | tgrade, data = GBSG2, ytrafo = function(...)
                  trafo(..., numeric_trafo = logrank_trafo, block = GBSG2$tgrade), 
                  teststat = "quad")


###################################################
### code chunk number 27: Peto
###################################################
survdiff(Surv(time, cens) ~ horTh + strata(tgrade), data = GBSG2, rho = 1)[c("chisq", "pvalue")]
(ft <- with(GBSG2, free1way(Surv(time, cens) ~ horTh | tgrade, 
                            link = "logit")))
summary(ft, test = "Rao")
summary(ft, test = "LRT")
summary(ft, test = "Wald")
summary(ft, test = "Permutation")


###################################################
### code chunk number 28: sw
###################################################
library("survival")
GBSG2$str <- cut(GBSG2$tsize, breaks = c(0, 1:9 * 10, Inf))
cm <- coxph(Surv(time, cens) ~ horTh + strata(str), data = GBSG2)

ft <- with(GBSG2, free1way(Surv(time, cens) ~ horTh | str, 
                           link = "cloglog"))
coef(cm)
coef(ft)
vcov(cm)
vcov(ft)
### Rao score tests
summary(cm)$sctest
summary(ft, test = "Rao")
### likelihood ratio tests
summary(cm)$logtest
summary(ft, test = "LRT")
### Wald tests
summary(cm)$waldtest
summary(ft, test = "Wald")
### asymptotic permutation tests
survdiff(Surv(time, cens) ~ horTh + strata(str), data = GBSG2, rho = 0)[c("chisq", "pvalue")]
summary(ft, test = "Permutation")


###################################################
### code chunk number 29: Peto
###################################################
survdiff(Surv(time, cens) ~ horTh + strata(str), data = GBSG2, rho = 1)[c("chisq", "pvalue")]
(ft <- with(GBSG2, free1way(Surv(time, cens) ~ horTh | str, 
                            link = "logit")))
summary(ft, test = "Rao")
summary(ft, test = "LRT")
summary(ft, test = "Wald")
summary(ft, test = "Permutation")


###################################################
### code chunk number 30: normal
###################################################
nd$y <- rnorm(nrow(nd))
free1way(y ~ g | s, data = nd, link = "probit")
independence_test(y ~ g | s, data = nd, ytrafo = function(...)
                  trafo(..., numeric_trafo = normal_trafo, block = nd$s), 
                  teststat = "quad")


###################################################
### code chunk number 31: friedman-data
###################################################
example(friedman.test, echo = FALSE)

### Myles Hollander & Wolfe (2014, Example 7.1, page 294)
boxplot(RoundingTimes, xlab = "Method", ylab = "Rounding-First-Base Time", 
        las = 1)
matplot(t(RoundingTimes), add = TRUE, type = "l", 
        lty = 1, lwd = 2, col = rgb(.1, .1, .1, .1))

me <- colnames(RoundingTimes)
d <- expand.grid(me = factor(me, labels = me, levels = me),
                 id = factor(seq_len(nrow(RoundingTimes))))
d$time <- c(t(RoundingTimes))


###################################################
### code chunk number 32: friedman
###################################################
friedman.test(RoundingTimes)
(ft <- free1way(time ~ me | id, data = d))


###################################################
### code chunk number 33: friedman-add
###################################################
summary(ft)
library("multcomp")
glht(ft, linfct = mcp(me = "Tukey"))


###################################################
### code chunk number 34: friedman-link
###################################################
logLik(ft)
logLik(free1way(time ~ me | id, data = d, link = "probit"))
logLik(free1way(time ~ me | id, data = d, link = "cloglog"))
logLik(free1way(time ~ me | id, data = d, link = "loglog"))


###################################################
### code chunk number 35: McNemar
###################################################
example(mcnemar.test, echo = FALSE)
# set-up data frame with survey outcomes for voters
s <- gl(2, 1, labels = dimnames(Performance)[[1L]])
survey <- gl(2, 1, labels = c("1st", "2nd"))
nvoters <- c(Performance)
x <- expand.grid(survey = survey, voter = factor(seq_len(sum(nvoters))))
x$performance <- c(rep(s[c(1, 1)], nvoters[1]), rep(s[c(2, 1)], nvoters[2]),
                   rep(s[c(1, 2)], nvoters[3]), rep(s[c(2, 2)], nvoters[4]))
# note that only those voters changing their minds are relevant
mcn <- free1way(xtabs(~ performance + survey + voter, data = x))
# same result as mcnemar.test w/o continuity correction
print(mcn)
# X^2 statistic
summary(mcn, test = "Permutation")$statistic^2
mcnemar.test(Performance, correct = FALSE)
# Wald inference
summary(mcn)
confint(mcn, test = "Wald")
### because the model is saturated, the link function doesn't affect the
### p-value (but the coefficients are of course different)
free1way(xtabs(~ performance + survey + voter, data = x), link = "probit")


###################################################
### code chunk number 36: incomplete
###################################################
data("taste", package = "daewr")
### highly discrete
table(taste$score)
summary(free1way(score ~ recipe | panelist, data = taste))


###################################################
### code chunk number 37: Tukey
###################################################
tk <- free1way(Ozone ~ Month, data = airquality)
library("multcomp")
confint(glht(tk, linfct = mcp(Month = "Tukey")))


###################################################
### code chunk number 38: plot-ex
###################################################
tk <- free1way(Ozone ~ Month, data = airquality)
plot(tk, las = 1)


###################################################
### code chunk number 39: plot-ex-cdf
###################################################
plot(tk, cdf = TRUE, model = FALSE, las = 1)


###################################################
### code chunk number 40: plot-ex-ecdf
###################################################
aq <- subset(airquality, !is.na(Ozone))
aq$r <- match(aq$Ozone, sort(unique(aq$Ozone)))
library("latticeExtra")
plot(ecdfplot(~ r, data = aq, groups = Month, col = 1:5))


###################################################
### code chunk number 41: ppplot
###################################################
y <- rlogis(50)
x <- rlogis(50, location = 3)
ppplot(y, x, conf.level = .95, las = 1)


###################################################
### code chunk number 42: ppplot-savage
###################################################
ppplot(y, x, conf.args = list(link = "cloglog", type = "Wald", 
                              col = NA, border = NULL),
       conf.level = .95, las = 1)


###################################################
### code chunk number 43: rfree1way
###################################################
(logOR <- c(log(1.5), log(2)))
nd <- rfree1way(150, delta = logOR)
coef(ft <- free1way(y ~ groups, data = nd))
sqrt(diag(vcov(ft)))
logLik(ft)
nd$y <- qchisq(nd$y, df = 3)
coef(ft <- free1way(y ~ groups, data = nd))
sqrt(diag(vcov(ft)))
logLik(ft)
N <- 25
pvals <- replicate(Nsim, 
{
  nd <- rfree1way(n = N, blocks = 2, delta = c(.25, .5), alloc_ratio = 2)
  summary(free1way(y ~ groups | blocks, data = nd), test = "Permutation")$p.value
})

power.free1way.test(n = N, blocks = 2, delta = c(.25, .5), alloc_ratio = 2)
mean(pvals < .05)


###################################################
### code chunk number 44: rfree1waysurv
###################################################
N <- 1000
nd <- rfree1way(N, delta = 1, link = "cloglog")
nd$C <- rfree1way(n = N, delta = 1, offset = -c(qlogis(.25), qlogis(.5)), 
                  link = "cloglog")$y
nd$y <- Surv(pmin(nd$y, nd$C), nd$y < nd$C)
### check censoring probability
1 - tapply(nd$y[,2], nd$groups, mean)
summary(free1way(y ~ groups, data = nd, link = "cloglog"))
summary(coxph(y ~ groups, data = nd))


###################################################
### code chunk number 45: power.prop.test
###################################################
delta <- log(1.5)
power.prop.test(n = 25, p1 = .5, p2 = plogis(qlogis(.5) - delta))
power.free1way.test(n = 25, prob = c(.5, .5), delta = delta)


###################################################
### code chunk number 46: power.odds.test
###################################################
prb <- matrix(c(.25, .25, .25, .25,
                .10, .20, .30, .40), ncol = 2)
colnames(prb) <- c("s1", "s2")
power.free1way.test(n = 20, prob = prb, 
                    strata_ratio = 2,
                    alloc_ratio = c(1.5, 2, 2), 
                    delta = log(c("low" = 1.25, "med" = 1.5, "high" = 1.75)))


###################################################
### code chunk number 47: wilcox
###################################################

delta <- log(3)
N <- 15
w <- gl(2, N)
pw <- numeric(Nsim)
for (i in seq_along(pw)) {
    y <- rlogis(length(w), location = c(0, delta)[w])
    pw[i] <- wilcox.test(y ~ w)$p.value
}
mean(pw < .05)

power.free1way.test(n = N, delta = delta)

### approximate formula in Hmisc::popower
library("Hmisc")
popower(p = rep(1 / N, N), odds.ratio = exp(delta), n = 2 * N)


###################################################
### code chunk number 48: kruskal
###################################################
delta <- c("B" = log(2), "C" = log(3))
N <- 15
w <- gl(3, N)
pw <- numeric(Nsim)
for (i in seq_along(pw)) {
    y <- rlogis(length(w), location = c(0, delta)[w])
    pw[i] <- kruskal.test(y ~ w)$p.value
}
mean(pw < .05)

power.free1way.test(n = N, delta = delta)


###################################################
### code chunk number 49: table
###################################################
prb <- rep.int(1, 4) / 4
pw <- numeric(Nsim)
cf <- matrix(0, nrow = Nsim, ncol = length(delta))
colnames(cf) <- names(delta)
for (i in seq_along(pw)) {
    nd <- rfree1way(n = N, prob = prb, delta = delta)
    ft <- free1way(y ~ groups, data = nd)
    cf[i,] <- coef(ft)
    pw[i] <- summary(ft, test = "Permutation")$p.value
}
mean(pw < .05)
boxplot(cf, las = 1, ylab = expression(hat(delta)))
points(c(1:2), delta, pch = 19, col = "red")
power.free1way.test(n = N, prob = prb, delta = delta)


###################################################
### code chunk number 50: stable
###################################################
prb <- cbind(S1 = rep(1, 4), 
             S2 = c(1, 2, 1, 2), 
             S3 = 1:4)
dimnames(prb) <- list(Ctrl = paste0("i", seq_len(nrow(prb))),
                      Strata = colnames(prb))

pw <- numeric(Nsim)
cf <- matrix(0, nrow = Nsim, ncol = length(delta))
colnames(cf) <- names(delta)
for (i in seq_along(pw)) {
    nd <- rfree1way(n = N, prob = prb, delta = delta)
    ft <- free1way(y ~ groups | blocks, data = nd)
    cf[i,] <- coef(ft)
    pw[i] <- summary(ft, test = "Permutation")$p.value
}
mean(pw < .05)
boxplot(cf, las = 1, ylab = expression(hat(delta)))
points(c(1:2), delta, pch = 19, col = "red")


###################################################
### code chunk number 51: powertest
###################################################
power.free1way.test(n = N, prob = prb, delta = delta, seed = 3)
power.free1way.test(power = .8, prob = prb, delta = delta, seed = 3)
power.free1way.test(n = 19, prob = prb, delta = delta, seed = 3)


###################################################
### code chunk number 52: Jeffreys
###################################################
N <- 20
w <- gl(2, N)
y <- rnorm(length(w), mean = c(-2, 3)[w])

x <- free1way(y ~ w, link = "probit")
coef(x)
logLik(x)

pll <- function(cf) {

    start <- x$par
    start[1] <- cf
    x$profile(start, fix = 1)
}

### https://doi.org/10.1111/j.0006-341X.2001.00114.x
### https://doi.org/10.1111/j.1467-9876.2012.01057.x
### https://doi.org/10.1186/s12874-017-0313-9
### https://files.osf.io/v1/resources/fet4d_v3/providers/osfstorage/682fb176db88f967facacb5a?format=pdf&action=download&direct&version=1
### https://doi.org/10.1002/sim.6537
### https://doi.org/10.1007/s11222-023-10217-3
### https://arxiv.org/abs/2510.06465
fun <- function(cf) {
    ret <- pll(cf)
    ret$value - .5 * determinant(ret$hessian, logarithm = TRUE)$modulus
}

ci <- confint(x, level = .99, test = "Wald")
grd <- seq(from = ci[1], to = ci[2], length.out = 50)

optim(coef(x), fn = fun, method = "Brent", 
      lower = min(grd), upper = max(grd))[c("par", "value")]


###################################################
### code chunk number 53: MPL_Jeffreys
###################################################
free1way(y ~ w, link = "probit", MPL_Jeffreys = TRUE)