R/fitPF.R

In LuchTiarna/PsyMetFuns: Package for extrapolating psychometric information

#'fitPF
#'
#'Provides the frame for psycometrics functions. Combines the sigmoid and core function.
#'@param gamma sets the loves boundary of function
#'@param lambda sets the highes boundary of function
#'@param sigmoid determines the sigmoid of the fuction
#'@param core dermines the core of the function
#'@param x the vector of level values
#'@param ... specifies the parametres or core function
#'@param type specifies, whether function is CDF of PDF type
#'@param inverse specifies, whether to compute the inverse function
#'
#'@return vector of return values
#'@export
fitPF <- function(formula, data, sigmoid, core, gamma=NULL, lambda=NULL, split_by=NULL, par=NULL, fn=NULL, gr=NULL, type="yes/no", ...,
                  method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
                  lower = -Inf, upper = Inf,
                  control = list(), hessian = FALSE){

  dd <- data.frame(levels=data[[formula[[3]]]], yes=data[[formula[[2]][[2]]]],no=data[[formula[[2]][[3]]]]) # creating a suiting representation of data
  if(tolower(type)=="hitpercentage"){
    #conversing to yes/no arrangement
    dd$yes <- dd$yes*dd$no
    dd$no <- dd$no - dd$yes
  }else if (tolower(type) != "yes/no"){
    warning("Unknown data organization. Use \"hitPercentage\" or \"yes/no\" notation.\n");return(NULL)
  }

  if(!all((dd$yes) == round(dd$yes)) || !all(dd$no == round(dd$no))){
    if(tolower(type) == "hitpercentage"){warning("The multiplication product of hitPercengate and obsNumbers should be an integer, because Yes and No responses counts should be integers.\n")}
    else{ warning("Yes and No responses counts should be integers.\n")}
  }

  #spliting data according to additional parameters
  if(is.null(split_by)){
    dd <- list(dd)
  }else{
    tryDD <- tryCatch({  split(x=dd, f=split_by, drop=TRUE, lex.order = FALSE) })

    if(!is.list(tryDD) || nrow(tryDD[[1]]) == nrow(dd)){
      error <- tryDD
      tryDD <- tryCatch({
        split_by <- as.list(data[unlist(split_by)])
        split(x=dd, f=split_by, drop = TRUE, lex.order = FALSE)
      })
    }
    if(!is.list(tryDD)){ stop(error, tryDD) }
    dd <- tryDD
  }

  #fits the data using different fixed parameters
  if(is.null(gamma) && is.null(lambda)){
    model <- mapply(FUN=fitPF.PF, dd, SIMPLIFY=FALSE, MoreArgs=list(sigmoid=sigmoid, core=core, par=par, fn=fn, gr=gr, ...,
                                                                    method=method,lower=lower, upper=upper, control=control, hessian=hessian))  }
  else if(!is.null(gamma) && is.null(lambda)){
    model <- mapply(FUN=fitPF.PF_fixed_gamma, dd, SIMPLIFY=FALSE, MoreArgs=list(sigmoid=sigmoid, core=core, par=par, fn=fn, gr=gr, gamma=gamma, ...,
                                                                                method=method,lower=lower, upper=upper, control=control, hessian=hessian)) }
  else if(is.null(gamma) && !is.null(lambda)){
    model <- mapply(FUN=fitPF.PF_fixed_lambda, dd, SIMPLIFY=FALSE, MoreArgs=list(sigmoid=sigmoid, core=core, par=par, fn=fn, gr=gr, lambda=lambda, ...,
                                                                             method=method,lower=lower, upper=upper, control=control, hessian=hessian))  }
  else {
    model <- mapply(FUN=fitPF.PF_fixed_gamma_lambda, dd, SIMPLIFY=FALSE, MoreArgs=list(sigmoid=sigmoid, core=core, par=par, fn=fn, gr=gr, gamma=gamma, lambda=lambda, ...,
                                                                                   method=method,lower=lower, upper=upper, control=control, hessian=hessian))  }

  if(is.null(split_by)){model <- model[[1]]} # when data is not splited returns only PF, not list of PF
  return(model)
}

fitPF.PF <- function(data, sigmoid, core, par=NULL, fn=NULL, gr=NULL, ...,
                     method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
                     lower = -Inf, upper = Inf,
                     control = list(), hessian = FALSE){
  model <- NULL
  if(is.null(fn)){ default_fn=TRUE
  fn=function(params){
    if(params[1]<0 || params[2]<0 || (params[1]+params[2] > 1)) {return(-Inf)}
    x <- PFunction(sigmoid, core, data$levels, params[1], params[2], params[-c(1,2)])
    if(is.nan(x)){return(-Inf)}

    if(length(x) != length(data$yes) || length(data$yes) != length(data$no) )
    {warning("All vectors must have the same length."); return(NaN)}
    pe <- data$yes*base::log(x)
    pe <- pe + data$no*base::log(1-x)

    return(-sum(pe))
  }
  }else{default_fn=FALSE} #if fn is not specified the most-likelihood function is used
  if(is.null(gr)){
    gr <- function(params){
      positiveCoef <- data$yes /      PFunction(sigmoid, core, data$levels, params[1], params[2], params[-c(1,2)])
      negativeCoef <- data$no  / (1 - PFunction(sigmoid, core, data$levels, params[1], params[2], params[-c(1,2)]))
      coef <- positiveCoef - negativeCoef

      gammaSlope <- sum(coef * PFunction(sigmoid, core, data$levels, params[1],  params[2], params[-c(1,2)], type="pdf_g"))
      lambdaSlope <- sum(coef * PFunction(sigmoid, core, data$levels, params[1], params[2], params[-c(1,2)], type="pdf_l"))

      gradient <- c(gammaSlope, lambdaSlope)
      for(n in 1:(length(params)-2)){
        parName <- paste("pdf_p", as.character(n), sep="")
        paramSlope <- sum(coef * PFunction(sigmoid,core, data$levels, params[1], params[2], params[-c(1,2)], type=parName))
        gradient <- c(gradient,paramSlope)
      }
      return(-gradient)
    }
  }
  if(is.null(par)){
    dat <- data.frame(level=data$level, hitPer=(data$yes)/(data$yes+data$no))
    ga <- min(dat$hitPer)
    la <- 1 - max(dat$hitPer)

    par=c(ga,la,rep(1, 2))
  } #TODO

  fit <- NULL
  fit <- tryCatch({ stats::optim(par=par, fn=fn, gr=gr, ..., method=method, lower=lower, upper=upper, control=control)})

  if(!is.list(fit)){return(fit)}

  model <- append(fit, list(sigmoid=sigmoid, core=core, gamma=fit$par[1], lambda=fit$par[2], params=c(fit$par[-c(1,2)])))
  model$par <- NULL
  class(model) <- c(class(model),"PF")
  return(model)
}

fitPF.PF_fixed_gamma <- function(data, sigmoid, core, gamma, par=NULL, fn=NULL, ...,
                                 method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
                                 lower = -Inf, upper = Inf,
                                 control = list(), hessian = FALSE){
  model <- NULL
  if(gamma<0){stop("Gamma must be in interval [0,1).")}

  if(is.null(fn)){ default_fn=TRUE
  fn=function(params){
    if(params[1]<0 || (gamma + params[1] > 1)) {return(NaN)}
    x <- PFunction(sigmoid, core, data$levels, gamma, params[1], params[-c(1)])

    if(length(x) != length(data$yes) || length(data$yes) != length(data$no) )
    {warning("All vectors must have the same length."); return(NaN)}
    pe <- data$yes*base::log(x)
    pe <- pe + (1-data$yes)*base::log(1-x)
    return(-sum(pe))
  }}else{default_fn=FALSE} #if fn is not specified the most-likelihood function is used

  if(is.null(fn)){ fn=function(params){
    if(params[1]<0 || (gamma + params[1] > 1)) {return(NaN)}
    x <- PFunction(sigmoid, core, data$levels, gamma, params[1], params[-c(1)])

    if(length(x) != length(data$yes) || length(data$yes) != length(data$no) )
    {warning("All vectors must have the same length."); return(NaN)}

    pe <- data$yes*base::log(x)
    pe <- pe + data$no*base::log(1-x)

    return(-sum(pe))
  }

  } #if fn is not specified the most-likelihood function is used
  if(is.null(par)){ par=c(0.05,rep(1, 2))} #TODO

  fit <- NULL
  fit <- tryCatch({ stats::optim(par, fn, ...) })

  if(!is.list(fit)){return(NULL)}

  model <- append(fit, list(sigmoid=sigmoid, core=core, gamma=gamma, lambda=fit$par[1], params=c(fit$par[-c(1)])))
  model$par <- NULL
  class(model) <- c(class(model),"PF")
  return(model)
  }

fitPF.PF_fixed_lambda <- function(data, sigmoid, core, lambda, par=NULL, fn=NULL, ...,
                                  method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
                                  lower = -Inf, upper = Inf,
                                  control = list(), hessian = FALSE){
  model <- NULL

  if(lambda<0){stop("Lambda must be in interval [0,1).")}

  if(is.null(fn)){ default_fn=TRUE
  fn=function(params){
    if(params[1]<0 || (lambda + params[1] > 1)) {return(NaN)}
    x <- PFunction(sigmoid, core, data$levels, params[1], lambda, params[-c(1)])

    if(length(x) != length(data$yes) || length(data$yes) != length(data$no) )
    {warning("All vectors must have the same length."); return(NaN)}
    pe <- data$yes*base::log(x)
    pe <- pe + (1-data$yes)*base::log(1-x)

    return(-sum(pe))
  }}else{default_fn=FALSE} #if fn is not specified the most-likelihood function is used

  if(is.null(par)){ par=c(0.05,rep(1, 2))} #TODO

  fit <- NULL
  fit <- tryCatch({ stats::optim(par, fn, ...) })

  if(!is.list(fit)){return(NULL)}

  model <- append(fit, list(sigmoid=sigmoid, core=core, gamma=fit$par[1], lambda=lambda, params=c(fit$par[-c(1)])))
  model$par <- NULL
  class(model) <- c(class(model),"PF")
  return(model)
}

fitPF.PF_fixed_gamma_lambda <- function(data, sigmoid, core, gamma, lambda, par=NULL, fn=NULL, ...,
                                        method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
                                        lower = -Inf, upper = Inf,
                                        control = list(), hessian = FALSE){
  model <- NULL
  if(gamma<0){stop("Gamma must be in interval [0,1).")}
  if(lambda<0){ stop("Gamma must be in interval [0,1).")}
  if(gamma + lambda > 1){stop("Summ of gamma and lambda must be in interval lesser than 1.")}

  if(is.null(fn)){ default_fn=TRUE
  fn=function(params){
    x <- PFunction(sigmoid, core, data$levels, gamma, lambda, params)

    if(length(x) != length(data$yes) || length(data$yes) != length(data$no) )
    {warning("All vectors must have the same length."); return(NaN)}
    pe <- data$yes*base::log(x)
    pe <- pe + (1-data$hitPercentage)*base::log(1-x)

    return(-sum(pe))
  }}else{default_fn=FALSE} #if fn is not specified the most-likelihood function is used

    if(is.null(par)){ par=rep(1, 2)} #TODO

  fit <- NULL
  fit <- tryCatch({ stats::optim(par, fn, ...) })

  if(!is.list(fit)){return(NULL)}

  model <- append(fit, list(sigmoid=sigmoid, core=core, gamma=gamma, lambda=lambda, params=c(fit$par)))
  model$par <- NULL
  class(model) <- c(class(model),"PF")
  return(model)
}