R/bptwin.R

In mets: Analysis of Multivariate Event Times

##' Liability-threshold model for twin data
##'
##' @aliases bptwin twinlm.time bptwin.time
##' @title Liability model for twin data
##' @seealso \code{\link{twinlm}}, \code{\link{twinlm.time}}, \code{\link{twinlm.strata}}, \code{\link{twinsim}}
##' @param x Formula specifying effects of covariates on the response.
##' @param data \code{data.frame} with one observation pr row. In
##'     addition a column with the zygosity (DZ or MZ given as a factor) of
##'     each individual much be
##'     specified as well as a twin id variable giving a unique pair of
##'     numbers/factors to each twin pair.
##' @param id The name of the column in the dataset containing the twin-id variable.
##' @param zyg The name of the column in the dataset containing the
##'     zygosity variable.
##' @param DZ Character defining the level in the zyg variable
##' corresponding to the dyzogitic twins. 
##' @param group Optional. Variable name defining group for interaction analysis (e.g., gender)
##' @param num Optional twin number variable
##' @param weights Weight matrix if needed by the chosen estimator (IPCW)
##' @param weights.fun Function defining a single weight each individual/cluster
##' @param strata Strata
##' @param messages Control amount of messages shown 
##' @param control Control argument parsed on to the optimization routine. Starting values may be parsed as '\code{start}'.
##' @param type Character defining the type of analysis to be
##'     performed. Should be a subset of "acde" (additive genetic factors, common
##'     environmental factors, dominant
##'     genetic factors, unique environmental factors).
##' @param eqmean Equal means (with type="cor")?
##' @param pairs.only Include complete pairs only?
##' @param stderr Should standard errors be calculated?
##' @param robustvar If TRUE robust (sandwich) variance estimates of the variance are used
##' @param p Parameter vector p in which to evaluate log-Likelihood and score function
##' @param indiv If TRUE the score and log-Likelihood contribution of each twin-pair
##' @param constrain Development argument
##' @param samecens Same censoring
##' @param allmarg Should all marginal terms be included
##' @param varlink Link function for variance parameters
##' @param ... Additional arguments to lower level functions
##' @author Klaus K. Holst
##' @export
##' @examples
##' data(twinstut)
##' b0 <- bptwin(stutter~sex,
##'              data=droplevels(subset(twinstut,zyg%in%c("mz","dz"))),
##'              id="tvparnr",zyg="zyg",DZ="dz",type="ae")
##' summary(b0)
bptwin <- function(x, data, id, zyg, DZ, group=NULL,
                   num=NULL,
                   weights=NULL,
                   weights.fun=function(x) ifelse(any(x<=0), 0, max(x)),
                   strata=NULL,
                   messages=1,
                   control=list(trace=0),
                   type="ace",
                   eqmean=TRUE,
                   pairs.only=FALSE,
                   samecens=TRUE,
                   allmarg=samecens&!is.null(weights),
                   stderr=TRUE,                  
                   robustvar=TRUE,
                   p, indiv=FALSE,
                   constrain,
                   varlink,
                   ...) {
    
###{{{ setup

  mycall <- match.call()
  formulaId <- unlist(Specials(x,"cluster"))
  formulaStrata <- unlist(Specials(x,"strata"))
  formulaSt <- paste("~.-cluster(",formulaId,")-strata(",paste(formulaStrata,collapse="+"),")")
  formula <- update(x,formulaSt)
  if (!is.null(formulaId)) {
    id <- formulaId
    mycall$id <- id
  }
  if (!is.null(formulaStrata)) strata <- formulaStrata
  mycall$formula <- formula
 
  if (!is.null(strata)) {
    dd <- split(data,interaction(data[,strata]))
    nn <- unlist(lapply(dd,nrow))
    dd[which(nn==0)] <- NULL
    if (length(dd)>1) {
      fit <- lapply(seq(length(dd)),function(i) {
        if (messages>0) message("Strata '",names(dd)[i],"'")
        mycall$data <- dd[[i]]
        eval(mycall)
      })
      res <- list(model=fit)
      res$strata <- names(res$model) <- names(dd)
      class(res) <- c("twinlm.strata","biprobit")
      res$coef <- unlist(lapply(res$model,coef))
      res$vcov <- blockdiag(lapply(res$model,vcov.biprobit))
      res$N <- length(dd)
      res$idx <- seq(length(coef(res$model[[1]])))
      rownames(res$vcov) <- colnames(res$vcov) <- names(res$coef)
      return(res)
    }
  }

##################################################
### No strata
  if (is.null(control$method)) {
      if (!samecens & !is.null(weights)) {
        control$method <- "bhhh"
      } else {
        if (requireNamespace("ucminf",quietly=TRUE)) {
          control$method <- "gradient"
        } else control$method <- "nlminb"
      }
  }
  if (tolower(type)=="cor") type <- "u"
  if (length(grep("flex",tolower(type)))>0) { type <- "u"; eqmean <- FALSE }

  yvar <- paste(deparse(formula[[2]]),collapse="")
  data <- data[order(data[,id]),]
  idtab <- table(data[,id])
  if (sum(idtab>2)) stop("More than two individuals with the same id ")

  if (pairs.only) {
    data <- data[as.character(data[,id])%in%names(idtab)[idtab==2],]
    idtab <- table(data[,id])
  }
  if (is.logical(data[,yvar])) data[,yvar] <- data[,yvar]*1
  if (is.factor(data[,yvar])) data[,yvar] <- as.numeric(data[,yvar])-1  

  idx2 <- NULL
  if (missing(DZ)) {
    DZ <- levels(as.factor(data[,zyg]))[1]
    message("Using '",DZ,"' as DZ",sep="")
  }
  OS <- NULL
  OSon <- FALSE
  if (!is.null(OS)) {
    idx2 <- which(data[,zyg]==OS)
    OSon <- TRUE
    if (length(idx2)==0) {
      warning("No OS twins found")
      OSon <- FALSE
    }
  }
  idx1 <- which(data[,zyg]==DZ) ## DZ
  if (length(idx1)==0) stop("No DZ twins found")
  idx0 <- which(!(data[,zyg]%in%c(DZ,OS))) ## MZ
  if (length(idx1)==0) stop("No MZ twins found")
  zyg2 <- rep(1,nrow(data)); zyg2[idx0] <- 0; zyg2[idx2] <- 2
  data[,zyg] <- zyg2 ## MZ=0, DZ=1, OS=2

  
  opt <- options(na.action="na.pass")
  Data <- cbind(model.matrix(formula,data),data[,c(yvar,id,zyg,weights,num)])
  options(opt)

  rnames1 <- setdiff(colnames(Data),c(yvar,id,weights,zyg,num))
  nx <- length(rnames1) 
  if (nx==0) stop("Zero design not allowed")

  bidx0 <- seq(nx)
  midx0 <- bidx0; midx1 <- midx0+nx
  dS0. <- rbind(rep(1,4),rep(1,4),rep(1,4)) ## MZ
  dS1. <- rbind(c(1,.5,.5,1),rep(1,4),c(1,.25,.25,1)) ## DZ
  dS2.  <- rbind(c(1,0,0,1),rep(1,4),c(1,0,0,1),c(0,1,1,0))
  
  ##mytr <- function(x) x; dmytr <- function(x) 1
  ##mytr <- function(x) x^2; dmytr <- function(x) 2*x
  ##mytr <- function(z) 1/(1+exp(-z)); dmytr <- function(z) exp(-z)/(1+exp(-z))^2
  ACDU <- sapply(c("a","c","d","e","u"),function(x) length(grep(x,tolower(type)))>0)
  ACDU <- c(ACDU,os=OSon)

  if (missing(varlink) || (!is.null(varlink) && varlink%in%"log")) {
      mytr <- exp; dmytr <- exp; myinvtr <- log
      trname <- "exp"; invtrname <- "log"
  } else {
      mytr <- myinvtr <- identity; dmytr <- function(x) rep(1,length(x))
      trname <- ""; invtrname <- ""
  }
  
  dmytr2 <- function(z) 4*exp(2*z)/(exp(2*z)+1)^2
  mytr2 <- tanh;  myinvtr2 <- atanh
  trname2 <- "tanh"; invtrname2 <- "atanh"

  if (OSon) {
    ## logit <- function(p) log(p/(1-p))
    ## tigol <- function(z) 1/(1+exp(-z))
    ## dlogit <- function(p) 1/(p*(1-p))
    ## dtigol <- function(z) tigol(z)^2*exp(-z)    
    ## mytr <- function(p) c(exp(p[-length(p)]),tigol(p[length(p)]))
    ## myinvtr <- function(z) c(log(z[-length(z)]),logit(z[length(z)]))
    ## dmytr <- function(p) c(exp(p[-length(p)]),dtigol(p[length(p)]))
    mytr <- function(x) c(exp(x[-length(x)]),mytr2(x[length(x)]))
    myinvtr <- function(z) c(log(z[-length(z)]),myinvtr2(z[length(z)]))
    dmytr <- function(x) c(exp(x[-length(x)]),dmytr2(x[length(x)]))
  }

  if (ACDU["u"]) {
    ##      datanh <- function(r) 1/(1-r^2)
    dmytr <- dmytr2
    mytr <- mytr2;  myinvtr <- myinvtr2
    trname <- trname2; invtrname <- invtrname2
    dS0 <- rbind(c(0,1,1,0))
    vidx0 <- 1
    vidx1 <- 2
    vidx2 <- 3
    dS2 <- dS1 <- dS0
    nvar <- length(vidx0)+length(vidx1)
    if (OSon) nvar <- nvar+length(vidx2)
  } else {
    nvar <- sum(ACDU[1:3])
    vidx0 <- vidx1 <- seq(nvar); vidx2 <- seq(nvar+1)
    if (OSon) nvar <- nvar+1
    dS0 <- dS0.[ACDU[1:3],,drop=FALSE]
    dS1 <- dS1.[ACDU[1:3],,drop=FALSE]
    dS2 <- dS2.[which(c(ACDU[1:3],TRUE)),,drop=FALSE]
  }  
  if (eqmean) {
    bidx2 <- bidx1 <- bidx0
  } else {
    bidx1 <- bidx0+nx
    bidx2 <- bidx1+nx    
    if (OSon) nx <- 3*nx else nx <- 2*nx;
  }
  
  vidx0 <- vidx0+nx; vidx1 <- vidx1+nx; vidx2 <- vidx2+nx
  vidx <- nx+seq_len(nvar)
  midx <- seq_len(nx)
  plen <- nx+nvar

  Am <- matrix(c(1,.5,.5,1),ncol=2)
  Dm <- matrix(c(1,.25,.25,1),ncol=2)
  Vm <- matrix(c(1,0,0,1),ncol=2)
  Rm <- matrix(c(0,1,1,0),ncol=2)

##################################################

  ## system.time(Wide <- reshape(as.data.frame(Data),idvar=c(id,zyg),timevar=time,direction="wide"))
  ##  system.time(Wide <- as.data.frame(fast.reshape(Data,id=c(id),sep=".")))
  Wide <- as.data.frame(fast.reshape(Data,id=c(id,zyg),sep=".",idcombine=FALSE,labelnum=TRUE))
  yidx <- paste(yvar,1:2,sep=".")
  rmidx <- c(id,yidx,zyg)
    
  W0 <- W1 <- W2 <- NULL
  if (!is.null(weights)) {
    widx <- paste(weights,1:2,sep=".")
    rmidx <- c(rmidx,widx)
    W0 <- as.matrix(Wide[which(Wide[,zyg]==0),widx,drop=FALSE])
    W1 <- as.matrix(Wide[which(Wide[,zyg]==1),widx,drop=FALSE])
    W2 <- as.matrix(Wide[which(Wide[,zyg]==2),widx,drop=FALSE])
  }
  XX <- as.matrix(Wide[,setdiff(colnames(Wide),rmidx)])
  XX[is.na(XX)] <- 0

  Y0 <- as.matrix(Wide[which(Wide[,zyg]==0),yidx,drop=FALSE])
  Y1 <- as.matrix(Wide[which(Wide[,zyg]==1),yidx,drop=FALSE])
  Y2 <- as.matrix(Wide[which(Wide[,zyg]==2),yidx,drop=FALSE])
  XX0 <- XX[which(Wide[,zyg]==0),,drop=FALSE]
  XX1 <- XX[which(Wide[,zyg]==1),,drop=FALSE]
  XX2 <- XX[which(Wide[,zyg]==2),,drop=FALSE]
  
##################################################

###}}} setup

###{{{ Mean/Var function

  ##Marginals etc.
  MyData0 <- ExMarg(Y0,XX0,W0,dS0,eqmarg=TRUE,allmarg=allmarg)
  MyData1 <- ExMarg(Y1,XX1,W1,dS1,eqmarg=TRUE,allmarg=allmarg)
  MyData2 <- ExMarg(Y2,XX2,W2,dS2,eqmarg=TRUE,allmarg=allmarg)

  N <- cbind(length(idx0),length(idx1),length(idx2)); 
  N <- cbind(N,
             2*nrow(MyData0$Y0)+if (!pairs.only) NROW(MyData0$Y0_marg) else 0, 
             2*nrow(MyData1$Y0)+if (!pairs.only) NROW(MyData1$Y0_marg) else 0,
             2*nrow(MyData2$Y0)+if (!pairs.only) NROW(MyData2$Y0_marg) else 0,
             NROW(MyData0$Y0),NROW(MyData1$Y0),NROW(MyData2$Y0))

  
  colnames(N) <- c("Total.MZ","Total.DZ","Total.OS","Complete.MZ","Complete.DZ","Complete.OS","Complete pairs.MZ","Complete pairs.DZ","Complete pairs.OS")
  rownames(N) <- rep("",nrow(N))
  if (!OSon) N <- N[,-c(3,6,9),drop=FALSE]
  
  if (samecens & !is.null(weights)) {
    MyData0$W0 <- cbind(apply(MyData0$W0,1,weights.fun))
    if (!is.null(MyData0$Y0_marg))
      MyData0$W0_marg <- cbind(apply(MyData0$W0_marg,1,weights.fun))

    MyData1$W0 <- cbind(apply(MyData1$W0,1,weights.fun))
    if (!is.null(MyData1$Y0_marg))
        MyData1$W0_marg <- cbind(apply(MyData1$W0_marg,1,weights.fun))

    MyData2$W0 <- cbind(apply(MyData2$W0,1,weights.fun))
    if (!is.null(MyData2$Y0_marg))
        MyData2$W0_marg <- cbind(apply(MyData2$W0_marg,1,weights.fun))
  }

  rm(Y0,XX0,W0,Y1,XX1,W1,Y2,XX2,W2)
  
  Sigma <- function(p0) {
    Sigma2 <- NULL
    p0[vidx] <- mytr(p0[vidx])    
    if (ACDU["u"]) {
      pos0 <- ifelse(OSon, plen-2, plen-1)
      Sigma0 <- diag(2) + p0[pos0]*Rm
      Sigma1 <- diag(2) + p0[pos0+1]*Rm
      if (OSon) Sigma2 <- diag(2) + p0[pos0+2]*Rm
    } else {
      ii <- ACDU; ii[4:5] <- FALSE
      pv <- ACDU*1;  pv[ii] <- p0[vidx]
      Sigma0 <- Vm*pv["e"] + pv["a"] + pv["c"] + pv["d"]
      Sigma1 <- Vm*pv["e"] + pv["a"]*Am + pv["c"] + pv["d"]*Dm
      Sigma2 <- Vm*pv["e"] + pv["c"] + (pv["a"] + pv["d"])*Vm +
        pv["os"]*(pv["a"] + pv["d"])*Rm
      if (OSon) {
        dS2 <- dS2.
        dS2[c(1,3),2:3] <- pv["os"]
        dS2[4,2:3] <- pv["a"]+pv["d"]
        dS2 <- dS2[which(c(ACDU[1:3],TRUE)),]
      }
    }
    return(list(Sigma0=Sigma0,Sigma1=Sigma1,Sigma2=Sigma2,dS2=dS2))
  }

  env <- new.env(parent=baseenv())
  #environment(mytr) <- baseenv()
  for (v in c("ACDU","mytr","vidx","OSon","plen",
              "Rm","Am","Vm","Dm","dS2.","dS2")) {
    assign(v, get(v), envir=env)
  }
  environment(Sigma) <- env

###}}} Mean/Var function
  
###{{{ U  

  p0 <- rep(-1,plen); ##p0[vidx] <- 0
  if (!missing(varlink) && is.null(varlink)) p0 <- rep(0.5,plen)
  if (OSon) p0[length(p0)] <- 0.3
  if (type=="u")
    p0[vidx] <- 0.3
  if (!is.null(control$start)) {
    p0 <- control$start
    control$start <- NULL
  } else {
    X <- rbind(MyData0$XX0[,midx0,drop=FALSE],MyData0$XX0[,midx1,drop=FALSE])
    Y <- rbind(MyData0$Y0[,1,drop=FALSE],MyData0$Y0[,2,drop=FALSE])
    g <- suppressWarnings(glm(Y~-1+X,family=binomial(probit)))
    p0[midx] <- coef(g)
  }
  
  U <- function(p,indiv=FALSE) {
    b0 <- cbind(p[bidx0])
    b1 <- cbind(p[bidx1])
    b2 <- cbind(p[bidx2])
    b00 <- b0; b11 <- b1; b22 <- b2
    ## if (bound) p[vidx] <- min(p[vidx],20)
    S <- Sigma(p)
    lambda <- eigen(S$Sigma0)$values
    if (any(lambda<1e-12 | lambda>1e9)) stop("Variance matrix out of bounds")
    
    Mu0 <- with(MyData0, cbind(XX0[,midx0,drop=FALSE]%*%b00,
                               XX0[,midx1,drop=FALSE]%*%b00))
    U0 <- with(MyData0, .Call("biprobit0",
                             Mu0,
                             S$Sigma0,dS0,Y0,XX0,W0,!is.null(W0),samecens,
                             PACKAGE="mets"))

    if (!is.null(MyData0$Y0_marg) &&!pairs.only) {
      mum <- with(MyData0, XX0_marg%*%b00)
      dSmarg <- dS0[,1,drop=FALSE]
       U_marg <- with(MyData0, .Call("uniprobit",
                                   mum,XX0_marg,
                                   S$Sigma0[1,1],t(dSmarg),Y0_marg,
                                   W0_marg,!is.null(W0_marg),TRUE,
                                   PACKAGE="mets"))
      U0$score <- rbind(U0$score,U_marg$score)
      U0$loglik <- c(U0$loglik,U_marg$loglik)
    }

    Mu1 <- with(MyData1, cbind(XX0[,midx0,drop=FALSE]%*%b11,
                               XX0[,midx1,drop=FALSE]%*%b11))

    U1 <- with(MyData1, .Call("biprobit0",
                             Mu1,
                             S$Sigma1,dS1,Y0,XX0,W0,!is.null(W0),samecens,
                             PACKAGE="mets"))
    if (!is.null(MyData1$Y0_marg) &&!pairs.only) {
      mum <- with(MyData1, XX0_marg%*%b11)
      dSmarg <- dS1[,1,drop=FALSE]
      U_marg <- with(MyData1, .Call("uniprobit",
                                   mum,XX0_marg,
                                   S$Sigma1[1,1],t(dSmarg),Y0_marg,
                                   W0_marg,!is.null(W0_marg),TRUE,
                                   PACKAGE="mets"))
      U1$score <- rbind(U1$score,U_marg$score)
      U1$loglik <- c(U1$loglik,U_marg$loglik)
    }

    U2 <- val2 <- NULL
    if (OSon) {
      Mu2 <- with(MyData2, cbind(XX0[,midx0,drop=FALSE]%*%b22,
                                 XX0[,midx1,drop=FALSE]%*%b22))
      U2 <- with(MyData2, .Call("biprobit0",
                               Mu2,
                               S$Sigma2,S$dS2,Y0,XX0,W0,!is.null(W0),samecens,
                               PACKAGE="mets"))
      if (!is.null(MyData2$Y0_marg) &&!pairs.only) {
        mum <- with(MyData2, XX0_marg%*%b22)
        dSmarg <- S$dS2[,1,drop=FALSE]
        U_marg <- with(MyData2, .Call("uniprobit",
                                     mum,XX0_marg,
                                     S$Sigma2[1,1],t(dSmarg),Y0_marg,
                                     W0_marg,!is.null(W0_marg),TRUE,
                                     PACKAGE="mets"))
        U2$score <- rbind(U2$score,U_marg$score)
        U2$loglik <- c(U2$loglik,U_marg$loglik)
      }
    }

    
    if (indiv) {
      ll0 <- U0$loglik
      ll1 <- U1$loglik
      val0 <- U0$score[MyData0$id,,drop=FALSE]
      val1 <- U1$score[MyData1$id,,drop=FALSE]
      
      N0 <- length(MyData0$id)
      idxs0 <- seq_len(N0)
      if (length(MyData0$margidx)>0) {
        for (i in seq_len(N0)) {
          idx0 <- which((MyData0$idmarg)==(MyData0$id[i]))+N0
          idxs0 <- c(idxs0,idx0)
          val0[i,] <- val0[i,]+colSums(U0$score[idx0,,drop=FALSE])
        }
        val0 <- rbind(val0, U0$score[-idxs0,,drop=FALSE])
        ll0 <- c(ll0,ll0[-idxs0])
      }
      
      N1 <- length(MyData1$id)
      idxs1 <- seq_len(N1)
      if (length(MyData1$margidx)>0) {
        for (i in seq_len(N1)) {
          idx1 <- which((MyData1$idmarg)==(MyData1$id[i]))+N1
          idxs1 <- c(idxs1,idx1)
          val1[i,] <- val1[i,]+colSums(U1$score[idx1,,drop=FALSE])
        }
        val1 <- rbind(val1, U1$score[-idxs1,,drop=FALSE])
        ll1 <- c(ll1,ll1[-idxs1])        
      }

      if (OSon) {
        ll2 <- U2$loglik
        val2 <- U2$score[MyData2$id,,drop=FALSE]
        N2 <- length(MyData2$id)
        idxs2 <- seq_len(N2)
        if (length(MyData2$margidx)>0) {
          for (i in seq_len(N2)) {
            idx2 <- which((MyData2$idmarg)==(MyData2$id[i]))+N2
            idxs2 <- c(idxs2,idx2)
            val2[i,] <- val2[i,]+colSums(U2$score[idx2,,drop=FALSE])
          }
          val2 <- rbind(val2, U2$score[-idxs2,,drop=FALSE])
          ll2 <- c(ll2,ll2[-idxs2])        
        }
      } 
      
      val <- matrix(0,ncol=plen,nrow=nrow(val0)+nrow(val1) + NROW(val2))
      val[seq_len(nrow(val0)),c(bidx0,vidx0)] <- val0
      val[nrow(val0)+seq_len(nrow(val1)),c(bidx1,vidx1)] <- val1
      if (OSon) {
        val[nrow(val0)+nrow(val1)+seq_len(nrow(val2)),c(bidx2,vidx2)] <- val2
      }
      trp <- dmytr(p[vidx])
      for (i in seq(length(vidx))) {
        val[,vidx[i]] <- val[,vidx[i]]*trp[i]
      }
      attributes(val)$logLik <- c(U0$loglik,U1$loglik,U2$loglik)
      return(val)
      
    }
    val <- numeric(plen)
    val[c(bidx0,vidx0)] <- colSums(U0$score)
    val[c(bidx1,vidx1)] <- val[c(bidx1,vidx1)]+colSums(U1$score)
    if (OSon) val[c(bidx2,vidx2)] <- val[c(bidx2,vidx2)]+colSums(U2$score)
    val[vidx] <- val[vidx]*dmytr(p[vidx])
    attributes(val)$logLik <- sum(U0$loglik)+sum(U1$loglik)+sum(U2$loglik)
    return(val)
  }

###}}} U

###{{{ optim

  if (!missing(p)) return(U(p,indiv=indiv))


  f <- function(p) crossprod(U(p))[1]
  f0 <- function(p) -sum(attributes(U(p))$logLik)
  g0 <- function(p) -as.numeric(U(p))
  h0 <- function(p) crossprod(U(p,indiv=TRUE))

  
  if (!missing(constrain)) {
    freeidx <- is.na(constrain)
    f <- function(p) {      
      p1 <- constrain; p1[freeidx] <- p
      res <- U(p1)[freeidx]
      crossprod(res)[1]
    }
    f0 <- function(p) {
      p1 <- constrain; p1[freeidx] <- p
      -sum(attributes(U(p1))$logLik)
    }
    g0 <- function(p) {
      p1 <- constrain; p1[freeidx] <- p
      -as.numeric(U(p1)[freeidx])
    }
    p0 <- p0[is.na(constrain)]
  }


  ## Derivatives, Sanity check 
  ## ff <- function(p) attributes(U(p,indiv=FALSE))$logLik
  ## pp <- c(0,-.1,.1,0.5)
  ## numDeriv::grad(ff,pp)
  ## U(pp,indiv=FALSE)

  
  controlstd <- list(hessian=0)
  controlstd[names(control)] <- control
  control <- controlstd
  
  nlminbopt <- intersect(names(control),c("eval.max","iter.max","trace","abs.tol","rel.tol","x.tol","step.min"))
  ucminfopt <- intersect(names(control),c("trace","grtol","xtol","stepmax","maxeval","grad","gradstep","invhessian.lt"))
  optimopt <- names(control) 

  if (control$method=="optimx") {
    if (!requireNamespace("optimx")) stop("Package `optimx` required.")
    if (is.null(control$lower)) control$lower <- -Inf
    if (is.null(control$upper)) control$upper <- Inf
    if (is.null(control$optimx.method)) control$optimx.method <- c("BFGS","Nelder-Mead")
  }

  op <- switch(tolower(control$method),
               nlminb=nlminb(p0,f0,gradient=g0,control=control[nlminbopt]),
               optim=optim(p0,fn=f0,gr=g0,control=control[ucminfopt]),
               optimx=optimx::optimx(p0, f0, g0, ...,
                                     lower=control$lower,
                                     upper=control$upper,
                                     method=control$optimx.method,
                                     control=control$optimx.control),
               ucminf=,
               quasi=,
               gradient=ucminf::ucminf(p0,fn=f0,gr=g0,control=control[ucminfopt],hessian=0),
               ## ,
               ## bhhh={
               ##   controlnr <- list(stabil=FALSE,
               ##                     gamma=0.1,
               ##                     gamma2=1,
               ##                     ngamma=5,
               ##                     iter.max=200,
               ##                     epsilon=1e-12,
               ##                     tol=1e-9,
               ##                     trace=1,
               ##                     stabil=FALSE)
               ##   controlnr[names(control)] <- control
               ##   lava:::NR(start=p0,NULL,g0, h0,control=controlnr)
##               },
               ##                 op <- switch(mycontrol$method,
               ##                              ucminf=ucminf(p0,f,control=mycontrol[ucminfopt],hessian=F),
               ##                optim=optim(p0,f,control=mycontrol[ucminfopt],...),
                 nlminb(p0,f,control=control[nlminbopt]))

  if (inherits(op, "optimx")) {
    val <- op$value
    if (length(val)==1) {
      pp <- as.numeric(coef(op))
    } else {
      idx <- which.min(op$value)
      pp <- as.numeric(coef(op)[idx,])
    }
    op <- list(par=pp, opt=op)
  }
  if (stderr) {
    UU <- U(op$par,indiv=TRUE)    
    I <- -numDeriv::jacobian(U,op$par)
    tol <- 1e-15
    iI <- Inverse(I,tol)
    V <- iI
    sqrteig <- attributes(V)$sqrteig
    J <- NULL
    if (robustvar) {
      J <- crossprod(UU)
      V <- iI%*%J%*%iI
    }
    if (any(sqrteig<tol)) warning("Near-singular covariance matrix (pseudo-inverse used)")
  } else {
    UU <- matrix(NA,ncol=length(op$par),nrow=1)
    I <- J <- V <- matrix(NA,ncol=length(op$par),nrow=length(op$par))
  }

###}}} optim

###{{{ return

  suppressWarnings(cc <- cbind(op$par,sqrt(diag(V))))
  cc <- cbind(cc,cc[,1]/cc[,2],2*(pnorm(abs(cc[,1]/cc[,2]),lower.tail=FALSE)))
  colnames(cc) <- c("Estimate","Std.Err","Z","p-value")
  vnames1 <- NULL
  trnam <- " "

  if (!eqmean) {
    rnam <- rnames1
    rnames1 <- c(paste(rnam,"MZ",sep=trnam),paste(rnam,"DZ",sep=trnam))
    if (OSon) rnames1 <- c(rnames1,paste(rnam,"OS",sep=trnam))

  }

  if (ACDU["u"]) {
    groups <- c("MZ","DZ"); if (OSon) groups <- c(groups,"OS")
    rnames <- c(rnames1,paste(c("atanh(rho)","atanh(rho)"),groups,sep=trnam))
  } else {
      rnames0 <- c("var(A)","var(C)","var(D)")[ACDU[1:3]]
      if (invtrname!="")
          rnames0 <- paste(invtrname,"(",rnames0,")",sep="")
      rnames <- c(rnames1,rnames0)
      if (OSon) rnames <- c(rnames,"atanh(rho(G[OS]))")
  }
  if (!missing(constrain)) rnames <- rnames[freeidx]
  rownames(cc) <- rnames
  rownames(V) <- colnames(V) <- rnames
  S <- Sigma(op$par)

  npar <- list(intercept=attributes(terms(formula))$intercept,
               pred=nrow(attributes(terms(formula))$factor)-1,
               var=sum(ACDU[-4]),
               ACDU=ACDU[-4]*1)
  
  npar[unlist(lapply(npar,length))==0] <- 0

  environment(mytr) <- baseenv()
  environment(myinvtr) <- baseenv()
  environment(dmytr) <- baseenv()

  val <- list(coef=cc,vcov=V,bread=iI,I=I,score=UU,logLik=attributes(UU)$logLik,opt=op,
              id=Wide[,id], model.frame=Wide,
              Sigma0=S$Sigma0, Sigma1=S$Sigma1, Sigma2=S$Sigma2,
              dS0=dS0, dS1=dS1, dS2=dS2,
              N=N,
              midx0=midx0, midx1=midx1,
              vidx0=vidx0, vidx1=vidx1, vidx2=vidx2,
              eqmean=eqmean, I=I,J=J, robustvar=robustvar,
              transform=list(tr=mytr, invtr=myinvtr, dtr=dmytr,
                  name=trname, invname=invtrname),
              SigmaFun=Sigma,
              npar=npar,
              OS=OSon
              )
  class(val) <- c("bptwin","biprobit")
  return(val)
}

###}}} return

##' @export
model.frame.bptwin <- function(formula,...) {
    formula$model.frame
}