fit1_binary_problem.R

# problem with fit1() with model="binary", in giving negative or NAs for LOD scores

library(qtl2)

dat <- readRDS("fit1_binary_problem.rds")
fit1(dat$p, dat$y, model="binary")$lod # NA value

stopifnot( all(names(dat$y) == rownames(dat$p)) )

out_glm <- glm(dat$y ~ -1 + dat$p, family=binomial(link=logit),
               control=list(maxit=500))
out0_glm <- glm(dat$y ~ 1, family=binomial(link=logit))
(out0_glm$deviance - out_glm$deviance)/2/log(10) # negative, but should be positive

source("logit_regr.R")
logitreg(dat$y, dat$p) # loglik = -Inf at iteration 15

fit1_binary_problem.rds

logit_regr.R

# crude logistic regression function

logitreg <-
    function(y, x=NULL, w=NULL, use_intercept=FALSE,
             tol=1e-8, maxit=1000, eta_tol= 100)
{
    if(is.null(x)) {
        x <- cbind(rep(1, length(y)))
    } else if(use_intercept) {
        x <- cbind(1, x)
    }
    if(is.null(w)) {
        w <- rep(1, length(y))
    }

    stopifnot(length(y) == nrow(x), length(y) == length(w))

    pi <- (y*w+0.5)/(w+1)
    cur_ll <- logitreg_ll(y, pi, w)
    cat(0, cur_ll, "\n")

    for(i in 1:maxit) {
        eta <- log(pi/(1-pi))
        W <- (w*pi*(1-pi))
        z <- eta + (y*w - w*pi)/W
        z[pi < tol] <- 0
        z[pi > 1-tol] <- 1

        out <- lm(z ~ -1 + x, weight=W)
        eta <- out$fitted
        pi <- exp(eta)/(1+exp(eta))
        pi[eta < -eta_tol] <- 0
        pi[eta > eta_tol] <- 1
#        cat("                         ", round(c(range(eta), range(pi)), 8), "\n")
        ll <- logitreg_ll(y, pi, w)
        cat(i, ll, abs(ll-cur_ll), "\n")
        if(abs(ll-cur_ll) < tol) break
        cur_ll <- ll
    }

    pi
}


logitreg_ll <-
    function(y, pi, w=NULL)
{
    if(is.null(w)) w <- rep(1, length(y))

    stopifnot(length(y)==length(pi),
              length(y)==length(w))

    zero <- y==0
    sum( log10(w) ) + sum( log10(pi[!zero]) ) + sum(log10(1-pi[zero]) )
}