R/NpdeData.R

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

##' Class "NpdeData" representing the structure of the longitudinal data
##' 
##' A longitudinal data structure
##' 
##' @name NpdeData-class
##' @aliases NpdeData NpdeData-class show,NpdeData-method print,NpdeData-method 
##' summary,NpdeData-method
##' [,NpdeData-method [<-,NpdeData-method
##' @docType class
##' @section Objects from the Class: NpdeData objects are typically created by calls to \code{\link{npdeData}} and contain the following slots:
##' 
##' \describe{
##' \item{name.data}{character string giving the name of the dataset}
##' \item{name.group}{character string giving the name of the grouping term (ID)}
##' \item{name.predictor}{character string giving the name of the predictor (X)}
##' \item{name.response}{character string giving the name of the response (Y)}
##' \item{name.cens}{character string giving the name of the censoring indicator}
##' \item{name.mdv}{character string giving the name of the missing data indicator}
##' \item{name.covariates}{vector of character string giving the name(s) of the covariates}
##' \item{name.ipred}{character string giving the name of the individual predictions}
##' \item{units}{(optional) a list with the units for X, Y, and covariates}
##' \item{data}{a dataframe containing the data}
##' \item{N}{number of subjects}
##' \item{ntot.obs}{total number of non-missing observations}
##' \item{nind.obs}{vector of size N giving the number of non-missing observations for each subject}
##' \item{ind}{index of non-missing observations}
##' \item{icens}{index of censored observations (non-missing)}
##' \item{not.miss}{a vector of boolean indicating for each observation whether it is missing (FALSE) or available (TRUE)}
##' \item{loq}{the censoring value}
##' }
##' @section Methods:
##' \describe{
##'   \item{npdeData(name.data):}{Create a new \code{\linkS4class{NpdeData}} object from dataset name.data}
##'   \item{print(npde.data):}{Prints a summary of object npde.data}
##'   \item{show(npde.data):}{Prints a short summary of object npde.data}
##'   \item{showall(npde.data):}{Prints a detailed summary of object npde.data}
##'   \item{plot(npde.data):}{Plots the data in npde.data. More details can be found in \code{\link{plot.NpdeData}}}
##'   \item{summary(npde.data):}{Returns a summary of object npde.data in list format}
##'   \item{set.npde.plotoptions(npde.data):}{Sets options for graphs of npde.data (internal method used in plots)}
##' }
##' @seealso \code{\link{npde}}, \code{\link{autonpde}}, \code{\link{plot.NpdeData}}
##' @keywords classes
##' @examples
##' 
##' methods(class="NpdeData")
##' 
##' showClass("NpdeData")
##' 
##' @exportClass NpdeData

setClass(
  Class="NpdeData",
  representation=representation(
    name.data="character",	# name of dataset
    name.group="character",	# name of column with ID
    name.predictor="character",# name of column(s) with predictors 
    name.response="character",	# name of column with response
    name.cens="character",	# name of column with censoring information
    name.miss="character",	# name of column indicating missing data (not censored)
    name.covariates="character",# name of column(s) with covariates (can be used for graphs)
    name.ipred="character",	# name of column indicating individual predictions (if available in the dataset)
    units="list",		# units (list with components for x, y, and cov), used in plots
    data="data.frame",          # the data (data frame with columns name.group (subject id), index (id renamed to 1:N), name.predictors (predictors), name.response (possibly transformed during fit), cens (1 if data is censored; the corresponding response is the censoring value), mdv (1 if the data is considered as missing) name.covariates, name.ipred (individual predictions, if availablea))
    									# binary covariates are modified to 0/1 
    ind="numeric",		# index of the non-missing observations # ECO TODO remove ?
    icens="numeric",		# index of the censored observations (non-missing)
    not.miss="logical",		# vector of logical, TRUE if present (=not missing), FALSE for missing data
    N="numeric",		# number of subjects
    ntot.obs="numeric",		# total number of observations
    nind.obs="numeric",		# number of observations for each subject
    loq="numeric"		# LOQ
    ),
  validity=function(object){
#    cat ("--- Checking NpdeData object ---\n")
    if (length(object@name.data)==0) {
      stop ("[ NpdeData : validation ] Please provide a name for the data (dataset or datafile on disk).")
    }
#     if([email protected]>0) {
#       N<[email protected]
#       if(N!=length(unique([email protected])) | N!=length(unique([email protected])) | N!=length([email protected])) {
#         cat("Size mismatch: id, index and/or nind.obs no longer correspond to the number of subjects.\n")
#       }
#       nobs<[email protected]
#       if(nobs!=dim([email protected])[1]) {
#         cat("Check length of predictor and response.\n")
#       }
#     }
    return(TRUE)
  }
)

##' Class "NpdeSimData" representing the structure of the longitudinal data
##' 
##' A longitudinal data structure, with simulated data
##' 
##' @name NpdeSimData-class
##' @aliases NpdeSimData NpdeSimData-class show,NpdeSimData-method [,NpdeSimData-method [<-,NpdeSimData-method
##' @docType class
##' @section Objects from the Class: NpdeSimData objects are created by associating an NpdeData object with matching simulated data, and they contain the following slots.
##' 
##' \describe{
##' \item{name.simdata}{character string giving the name of the dataset}
##' \item{nrep}{number of replications)}
##' \item{datsim}{a dataframe containing the simulated data,  with columns: idsim (subject id), irsim (replication index), xsim (simulated x), ysim (simulated response). After a call to \code{\link{npde}} or \code{\link{autonpde}}, an additional column ydsim (decorrelated replicated data) will be added.}
##' }
##' @section Methods:
##' \describe{
##'   \item{print(npde.simdata):}{Prints a summary of object npde.simdata}
##'   \item{show(npde.simdata):}{Prints a short summary of object npde.simdata}
##'   \item{showall(npde.simdata):}{Prints a detailed summary of object npde.simdata}
##' }
##' @seealso \code{\link{npde}}, \code{\link{autonpde}}
##' @keywords classes
##' @examples
##' 
##' showClass("NpdeSimData")
##' 
##' @export

setClass(
  Class="NpdeSimData",
  representation=representation(
    name.simdata="character",	# name of dataset
    nrep="numeric",		# number of replications
    datsim="data.frame"	# data (a dataframe with columns: idsim (subject id), irsim (replication index), xsim (simulated x), ysim (simulated response), ydsim (decorrelated replicated data))
  ),
  validity=function(object){
#    cat ("--- Checking NpdeData object ---\n")
    if (length(object@name.simdata)==0) {
      stop ("[ NpdeData : validation ] Please provide a name for the data (dataset or datafile on disk).")
    }
    return(TRUE)
  }
)

###############################
# ECO validity ne semble pas etre appele automatiquement quand on cree un objet => il faut l'appeler dans initialize

setMethod(
  f="initialize",
  signature="NpdeData",
  definition= function (.Object,name.data,name.group,name.predictor, name.response,name.covariates,name.cens,name.miss,name.ipred,units,data){
#    cat ("--- initialising NpdeData Object --- \n")
    if(missing(name.data) || length(name.data)==0) stop ("Please provide a name for the data (dataset or datafile on disk).")
    .Object@name.data<-name.data
    if(missing(name.group)) name.group<-character()
# ECO TODO: reconnaissance automatique (avant affectation a la valeur 2) ?
    if(missing(name.predictor)) name.predictor<-character()
    if(missing(name.response)) name.response<-character()
    if(missing(name.covariates) || length(name.covariates)==0 || name.covariates[1]=="") name.covariates<-character()
    if(missing(name.cens) || length(name.cens)==0 || name.cens=="") name.cens<-character()
    if(missing(name.miss) || length(name.miss)==0 ||name.miss=="") name.miss<-character()
    if(missing(name.ipred) || length(name.ipred)==0 ||name.ipred=="") name.ipred<-character()
    .Object@name.group<-name.group
    .Object@name.predictor<-name.predictor
    .Object@name.response<-name.response
    .Object@name.covariates<-name.covariates
    .Object@name.cens<-name.cens
    .Object@name.miss<-name.miss
    .Object@name.ipred<-name.ipred
    if(missing(units)) units<-list(x="-",y="-")
    if(is.null(units$x)) units$x<-"-"
    if(is.null(units$y)) units$y<-"-"
    ncov<-length(name.covariates)
    if(ncov>0) {
      nunit<-length(units$covariates)
      if(nunit==0) units$covariates<-rep("-",ncov)
      if(nunit>ncov) units$covariates<-units$covariates[1:ncov]
      if(nunit<ncov) {
        length(units$covariates)<-ncov
        units$covariates[(nunit+1):ncov]<-"-"
      }
    }
    .Object@units<-units
    .Object@N<-0
# Object validation
    validObject(.Object)
    return (.Object )
  }
)

setMethod(
  f="initialize",
  signature="NpdeSimData",
  definition= function (.Object,name.simdata){
#    cat ("--- initialising NpdeSimData Object --- \n")
    if(missing(name.simdata)) stop ("Please provide a name for the simulated dataset (dataset or datafile on disk) as a character string.")
    .Object@name.simdata<-name.simdata
    return (.Object )
  }
)

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

##' Get/set methods for NpdeData object
##' 
##' Access slots of a NpdeData using the object["slot"] format
##' 
##' @keywords methods
##' @exportMethod [
##' @exportMethod [<-

#### NpdeData
# Getteur
setMethod(
  f ="[",
  signature = "NpdeData" ,
  definition = function (x,i,j,drop ){
  switch (EXPR=i,
    "name.data"={return(x@name.data)},
    "header"={return(x@header)},
    "sep"={return(x@sep)},
    "na"={return(x@na)},
    "name.group"={return(x@name.group)},
    "name.predictor"={return(x@name.predictor)},
    "name.response"={return(x@name.response)},
    "name.cens"={return(x@name.cens)},
    "name.miss"={return(x@name.miss)},
		"name.covariates"={return(x@name.covariates)},
    "name.ipred"={return(x@name.ipred)},
    "units"={return(x@units)},
    "loq"={return(x@loq)},
    "data"={return(x@data)},
  	"ind"={return(x@ind)},
  				"icens"={return(x@icens)},
  				"not.miss"={return(x@not.miss)},
  	"N"={return(x@N)},
  	"ntot.obs"={return(x@ntot.obs)},
    "nind.obs"={return(x@nind.obs)},
    stop("No such attribute\n")
   )
  }
)

# Setteur
setReplaceMethod(
  f ="[",
  signature = "NpdeData" ,
  definition = function (x,i,j,value){
  switch (EXPR=i,
    "name.data"={x@name.data<-value},
    "name.group"={x@name.group<-value},
    "name.predictor"={x@name.predictor<-value},
    "name.response"={x@name.response<-value},
    "name.cens"={x@name.cens<-value},
    "name.miss"={x@name.miss<-value},
  	"name.covariates"={x@name.covariates<-value},
  	"name.ipred"={x@name.ipred<-value},
    "units"={x@units<-value},
    "loq"={x@loq<-value},
    "data"={x@data<-value},
  	"ind"={x@ind<-value},	
  				"icens"={x@icens<-value},	
  				"not.miss"={x@not.miss<-value},		
  	"N"={x@N<-value},
  	"ntot.obs"={x@ntot.obs<-value},
    "nind.obs"={x@nind.obs<-value},
    stop("No such attribute\n")
   )
   validObject(x)
   return(x)
  }
)

#### Simulated data
# Getteur
setMethod(
  f ="[",
  signature = "NpdeSimData" ,
  definition = function (x,i,j,drop ){
  switch (EXPR=i,
    "name.simdata"={return(x@name.simdata)},
    "nrep"={return(x@nrep)},
    "datsim"={return(x@datsim)},
    stop("No such attribute\n")
   )
  }
)

# Setteur
setReplaceMethod(
  f ="[",
  signature = "NpdeSimData" ,
  definition = function (x,i,j,value){
  switch (EXPR=i,
    "name.simdata"={x@name.simdata<-value},
    "nrep"={x@nrep<-value},
    "datsim"={x@datsim<-value},
    stop("No such attribute\n")
   )
   validObject(x)
   return(x)
  }
)

##################################################################################
# print/show/showall
# alias in class documentation

##' @S3method print NpdeData
print.NpdeData <- function(x,nlines=10,...) {
    digits<-2;nsmall<-2
    cat("Object of class NpdeData\n")
    cat("    longitudinal data\n")
    if(length(x@name.data)>0)
    cat("Dataset",x@name.data,"\n")
    if(length(x@name.group)>0) {
    st1<-paste(x@name.response," ~ ",paste(x@name.predictor,collapse=" + ")," | ", x@name.group,sep="")
    cat("    Structured data:",st1,"\n")
    cat("    predictor:",x@name.predictor,paste("(",x@units$x,")",sep=""),"\n")
    if(length(x@name.covariates)>0) {
      cat("    covariates:",paste(paste(x@name.covariates," (",x@units$covariates,")",sep=""),collapse=", "),"\n")
    }
    if(dim(x@data)[1]>0) {
      if(nlines==0) return()
      cat("Dataset characteristics:\n")
        cat("    number of subjects:    ",x@N,"\n")
        cat("    number of non-missing observations:",x@ntot.obs,"\n")
        cat("    average/min/max nb obs:",format(mean(x@nind.obs),digits=digits, nsmall=nsmall), " / ", min(x@nind.obs)," / ",max(x@nind.obs),"\n")
        if(length(x@loq)>0) cat("      LOQ:    ",x@loq,"\n")
#    if(length([email protected])>0) print([email protected])
      if(nlines==(-1)) {
        cat("Data:\n")
        print(x@data)
      } else {
        cat("First",nlines,"lines of data:\n")
        nrowShow <- min (nlines , nrow(x@data))
        print(x@data[1:nrowShow,])
      }
    } else cat("No data.\n")
  } else cat("Empty object\n")
}

##' @S3method print NpdeSimData

print.NpdeSimData<-function(x,nlines=10,...) {
    cat("Object of class NpdeSimData\n")
    cat("    dataset:               ",x@name.simdata,"\n")
    if(length(x@nrep)>0) cat("    number of replications:",x@nrep,"\n") else cat("    x currently empty\n")
    if(length(x@datsim)>0) {
      if(nlines==0) return()
      if(nlines==(-1)) {
        cat("Data:\n")
        print(x@datsim)
      } else {
        cat("First",nlines,"lines of simulated data:\n")
        nrowShow <- min (nlines , nrow(x@datsim ))
        print(x@datsim[1:nrowShow,])
      }
    } else cat("No data.\n")
}

setMethod("show","NpdeData",
# show.NpdeData<-
	function(object) {
    cat("Object of class NpdeData\n")
    if(length(object@ntot.obs)==0) cat("    no data\n") else {
    cat("Dataset",object@name.data,"\n")
    st1<-paste(object@name.response," ~ ",paste(object@name.predictor,collapse=" + ")," | ", object@name.group,sep="")
    cat("    Structured data:",st1,"\n")
    if(length(object@name.covariates)>0) cat("    Covariates:",object@name.covariates,"\n")
      cat("This object has the following components:\n") 
      cat("     data: data\n")
      cat("     with",object@N,"subjects\n")
      cat("     ",object@ntot.obs,"observations\n")
      cat("The data has the following components\n")
      cat("     X:",object@name.predictor,"\n")
      cat("     Y:",object@name.response,"\n")
      if(length(object@name.ipred)>0) cat("     individual model predictions:", object@name.ipred,"\n")
      if(length(object@name.miss)>0) cat("     missing data:",object@name.miss," (1=missing)\n")
      if(length(object@name.cens)>0) cat("     censored data:",object@name.cens," (1=censored)\n")
      if(length(object@loq)>0) cat("      LOQ:    ",object@loq,"\n")
#      cat("     \n")
  }
  }
)

#show.NpdeSimData<-
setMethod("show","NpdeSimData",
	function(object) {
    cat("Object of class NpdeSimData\n")
    cat("    dataset:               ",object@name.simdata,"\n")
    if(length(object@nrep)>0) cat("    number of replications:",object@nrep,"\n") else cat("    object currently empty\n")
  }
)


#
#' Brief summary of an object
#'
#' Prints a brief summary of an object
#' 
#' @name showall
#' @aliases showall showall.NpdeData showall,NpdeData-method showall.NpdeSimData showall,NpdeSimData-method  showall.NpdeRes  showall.NpdeObject
#' @param object an object
#' @keywords print
##' @S3method showall NpdeData
##' @export
# Could be print, with only head of data

#setMethod("show","NpdeSimData",

#setMethod("showall","NpdeData",
showall.NpdeData <- function(object) {
    digits<-2;nsmall<-2
    cat("Object of class NpdeData\n")
    cat("    longitudinal data\n")
    if(length(object@name.data)>0) {
    cat("Dataset",object@name.data,"\n")
    st1<-paste(object@name.response," ~ ",paste(object@name.predictor,collapse=" + ")," | ", object@name.group,sep="")
    cat("    Structured data:",st1,"\n")
    cat("    subject identifier:    ",object@name.group,"\n")
    cat("    predictor:       ",object@name.predictor, paste("(",object@units$x,")",sep=""),"\n")
    cat("    response:        ",object@name.response,paste("(",object@units$y,")",sep=""),"\n")
    if(length(object@name.covariates)>0) {
    	cat("    covariates:",paste(paste(object@name.covariates," (", object@units$covariates,")",sep=""),collapse=", "),"\n")
    }
    cat("This object has the following components:\n") 
    cat("     data: data\n")
    cat("     with",object@N,"subjects\n")
    cat("     ",object@ntot.obs,"observations\n")
    cat("The data has the following components\n")
    cat("     X:",object@name.predictor,"\n")
    cat("     Y:",object@name.response,"\n")
    if(length(object@name.ipred)>0) cat("     individual model predictions:", object@name.ipred,"\n")
    if(length(object@name.miss)>0) cat("     missing data:",object@name.miss," (1=missing)\n")
    if(length(object@name.cens)>0) cat("     censored data:",object@name.cens," (1=censored)\n")
    if(length(object@loq)>0) cat("      LOQ:    ",object@loq,"\n")
    cat("Dataset characteristics:\n")
    cat("    number of subjects:    ",object@N,"\n")
    if(object@N>0) {
      cat("    number of non-missing observations:",object@ntot.obs,"\n")
      cat("    average/min/max nb obs:",format(mean(object@nind.obs),digits=digits, nsmall=nsmall), " / ", min(object@nind.obs)," / ",max(object@nind.obs),"\n")
#    if(length([email protected])>0) print([email protected])
    }
    if(dim(object@data)[1]>0) {
      cat("First lines of data:\n")
      nrowShow <- min (10 , nrow(object@data))
      print(object@data[1:nrowShow,])
    } else cat("No data.\n")
  } else cat("Empty object\n")
    }
#)

##' @S3method showall NpdeSimData

#setMethod("showall","NpdeSimData",
showall.NpdeSimData<-
					function(object) {
						cat("Object of class NpdeSimData\n")
						cat("    dataset:               ",object@name.simdata,"\n")
						if(length(object@nrep)>0) cat("    number of replications:",object@nrep,"\n") else cat("    object currently empty\n")
					}
#)

##' @S3method summary NpdeData

summary.NpdeData <- function(object, print=TRUE, ...) {
	  if(length(object@data)==0) {
			cat("Object of class NpdeData, empty.\n")
			return()
		}
	  res<-list(N=object@N,data=object@data, ntot.obs=object@ntot.obs,nind.obs=object@nind.obs)
	  if(length(object@loq)>0) res$loq<-object@loq
	  invisible(res)
}
##################################################################################

#
#' Set graph options
#'
#' This function is used to set options for graphs
#' 
##' @aliases set.npde.plotoptions set.npde.plotoptions,NpdeData-method
#' @usage set.npde.plotoptions(object, ...)
#' @param object an object of class NpdeData or NpdeObject
#' @param \dots arguments to replace default arguments (currently ignored)
#' @return a list of options for graphs
#' @details See documentation for a list of available options. 
#' @author Emmanuelle Comets <[email protected]@bichat.inserm.fr>
#' @seealso \code{\link{npde}}, \code{\link{autonpde}}
#' @keywords plot 
##' @S3method set.npde.plotoptions NpdeData
##' @export

set.npde.plotoptions.NpdeData<-function(object, ...) {
# setting default plot options
  plot.opt<-list(
# Options for plot types
    ilist=c(1:object["N"]),
    level=0:1,
    plot.loq=TRUE,			# Whether data under the LOQ should be plotted
    line.loq=TRUE,			# Whether an horizontal line should be plotted at Y=LOQ
    # General graphical options
    new=TRUE,				# whether a new page should be called
    ask=FALSE,				# whether the program should ask before creating a new page
    mfrow=c(),				# page layout (if empty, defaults to the default layout for each graph type)
    main="",				# title
    xlab="",
    ylab="",
    col="black",
    pch=20,
    lty=1,
    lwd=1,
    xlim=c(),
    ylim=c(),
    xlog=FALSE,
    ylog=FALSE,
    type="b",
    cex=1,
    cex.axis=1,
    cex.lab=1,
    cex.main=1,
    obs.pch=20,				# symbol for observed/censored data
    cens.pch=8,
    obs.col="black",		# color to plot observed/censored data
    cens.col="red",
    lcol="black",
    ablinecol="black",
    ablinelty=2,
    ablinelwd=2
    )
        
     if(is.null(plot.opt$name.X))
        plot.opt$name.X<-object["name.predictor"][1]
    plot.opt$xlab<-paste(plot.opt$name.X," (",plot.opt$units$x,")", sep="")
     if(length(object["name.response"])>0)
    plot.opt$ylab<-paste(object["name.response"]," (",object["units"]$y,")", sep="")
   return(plot.opt)
}


#
#' Replace graph options
#'
#' This function is used to replace graph options (available in the prefs slot of the NpdeData object) for plotting NpdeData objects
#' 
#' @usage replace.npde.plotoptions(plot.opt,...)
#' @param plot.opt a list of graphical preferences
#' @param \dots names and values of the options to be replaced
#' 
#' @return an updated list of options for graphs
#' @details See documentation for a list of available options. During replacement, invalid (including misspelled) options will raise warnings. 
#' @author Emmanuelle Comets <[email protected]@bichat.inserm.fr>
#' @seealso \code{\link{npde}}, \code{\link{autonpde}}
#' @keywords plot internal

replace.npde.plotoptions<-function(plot.opt,...) {
	args1<-match.call(expand.dots=TRUE)
# These arguments are used by other functions and may be passed on via "...", so we want to ignore them. Other arguments not in list will raise warnings
	legacy<-c("plot.type","namsav","namgr","loq")
	if(length(args1)>2) {
		# General arguments: col, pch
		i1<-match("col",names(args1))
		if(!is.na(i1)) {
			plot.opt$col<-eval(args1[[i1]])
			if(is.na(match("pcol",names(args1)))) plot.opt$pcol<-eval(args1[[i1]])
			if(is.na(match("lcol",names(args1)))) plot.opt$lcol<-eval(args1[[i1]])
			if(is.na(match("obs.col",names(args1)))) plot.opt$obs.col<-eval(args1[[i1]])
		}
		i1<-match("pch",names(args1))
		if(!is.na(i1)) {
			plot.opt$pch<-eval(args1[[i1]])
			if(is.na(match("obs.pch",names(args1)))) plot.opt$obs.pch<-eval(args1[[i1]])
		}
		# Other arguments
		for(i in 3:length(args1)) {
			if(match(names(args1)[i],names(plot.opt),nomatch=0)>0) {
				#    plot.opt[[names(args1)[i]]]<-args1[[i]] else {
				if(!is.null(eval(args1[[i]]))) plot.opt[[names(args1)[i]]]<-eval(args1[[i]])
				} else {
					if(is.na(match(names(args1)[i],legacy))) cat("Argument",names(args1)[i],"not available, check spelling.\n")
				}
			}
		}
	return(plot.opt)
	}

#' Plots a NpdeData object
#'
#' Plots the data in a NpdeData object
#' 
#' @param x a NpdeData object
#' @param y unused, here for compatibility with the base plot function
#' @param \dots additional graphical parameters to be passed on to the plot
#' @details The default plot is a spaghetti plot of all the data, with a line joining the observations for each subject. If censored data is present, it is shown with a different symbol and colour.
#' @references K. Brendel, E. Comets, C. Laffont, C. Laveille, and F.Mentre. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide. \emph{Pharmaceutical Research}, 23:2036--49, 2006.
##' @seealso \code{\link{set.npde.plotoptions}}
#' @keywords plot
#' @examples
#' 
#' data(theopp)
#' 
#' x<-npdeData(theopp,name.group="ID",name.predictor="Time",name.response="Conc", 
#' name.covariates=c("Wt"),units=list(x="hr",y="mg/L",covariates="kg"))
#' plot(x)
#' 
##' @importFrom graphics plot
##' @method plot NpdeData
##' @export

plot.NpdeData <- function(x, y, ...) {
# Plot the data, either as points or as lines grouped by [email protected]
    if(length(x@data)>0) {
    args1<-match.call(expand.dots=TRUE)
    i1<-match("type",names(args1))
    if(!is.na(i1)) {
      plot.type<-as.character(args1[[i1]])
      plot.type<-plot.type[plot.type!="c"]
    } else plot.type<-"b"
    plot.opt<-set.npde.plotoptions(x)
    plot.opt$new<-TRUE
    plot.opt$xlab<-paste(x@name.predictor," (",x@units$x,")",sep="")
    plot.opt$ylab<-paste(x@name.response," (",x@units$y,")",sep="")
    plot.opt<-replace.npde.plotoptions(plot.opt,...)
    logtyp<-paste(ifelse(plot.opt$xlog,"x",""),ifelse(plot.opt$ylog,"y",""),sep="")
    if(plot.opt$new) par(mfrow=c(1,1))
    tab<-x@data[x@ind,] # remove missing data
    has.cens<-(length(x@name.cens)>0)
    if(has.cens && max(tab[,x@name.cens])==0) has.cens<-FALSE
      if(plot.type=="p" | plot.type=="b") {
      	if(has.cens) {
        plot(tab[,x@name.predictor],tab[,x@name.response],xlab=plot.opt$xlab, ylab=plot.opt$ylab,col=plot.opt$col,log=logtyp,xlim=plot.opt$xlim, ylim=plot.opt$ylim,main=plot.opt$main,cex=plot.opt$cex,cex.axis=plot.opt$cex.axis, cex.lab=plot.opt$cex.lab,type="n")
        points(tab[tab[,x@name.cens]==0,x@name.predictor], tab[tab[,x@name.cens]==0,x@name.response],col=plot.opt$obs.col, pch=plot.opt$obs.pch,cex=plot.opt$cex)
        points(tab[tab[,x@name.cens]==1,x@name.predictor], tab[tab[,x@name.cens]==1,x@name.response],col=plot.opt$cens.col, pch=plot.opt$cens.pch,cex=plot.opt$cex)
      	} else 
      		plot(tab[,x@name.predictor],tab[,x@name.response],xlab=plot.opt$xlab, ylab=plot.opt$ylab,col=plot.opt$col,pch=plot.opt$pch,log=logtyp,xlim=plot.opt$xlim, ylim=plot.opt$ylim,main=plot.opt$main,cex=plot.opt$cex,cex.axis=plot.opt$cex.axis, cex.lab=plot.opt$cex.lab)
      }
      if(plot.type=="l") {
        plot(tab[,x@name.predictor],tab[,x@name.response],xlab=plot.opt$xlab, ylab=plot.opt$ylab,col=plot.opt$col,lty=plot.opt$lty,lwd=plot.opt$lwd,type="n", log=logtyp,xlim=plot.opt$xlim,ylim=plot.opt$ylim,main=plot.opt$main, cex=plot.opt$cex,cex.axis=plot.opt$cex.axis, cex.lab=plot.opt$cex.lab)
      }
      if(plot.type=="l" | plot.type=="b") {
        for(isuj in unique(tab[,x@name.group])) {
          lines(tab[tab[,x@name.group]==isuj,x@name.predictor], tab[tab[,x@name.group]==isuj,x@name.response],col=plot.opt$lcol, lty=plot.opt$lty,lwd=plot.opt$lwd)
      }
      }
    if(has.cens) abline(h=x@loq,col=plot.opt$ablinecol,lty=plot.opt$ablinelty, lwd=plot.opt$ablinelwd)

    } else cat("No data to plot.\n")
  }

##################################################################################
# @usage read.npdeData(object, header, sep, na.strings, detect)

#' Read observed data
#' 
#' Creates an NpdeData object containing the observed data, either from disk or from a dataframe
#' 
#'
#' @param header boolean indicating whether the file has a header (mandatory if 
#' detect is TRUE)
#' @param sep field separator (for files on disk)
#' @param na.strings strings to be considered as indicating NA
#' @param detect a boolean; if TRUE, automatic recognition of names will be attempted to detect necessary items (longitudinal data structure, missing data and censoring information)
#' @return an object of class \code{"\linkS4class{NpdeData}"}
#' @keywords internal

setMethod("read.npdeData","NpdeData",
  function(object,header,sep,na.strings,detect,verbose) {
    ow <- options("warn")
    options("warn"=-1)
# ce test devrait aller dans la definition de la classe
    if(class(object@name.data)!="character") {
    cat("Please provide the name of the data (data.frame or path to file on disk) as a character string.\n")
    return("Creation of npdeData failed")
  }
    if(exists(object@name.data)) {
      if(verbose) cat("Using the object called",object@name.data,"in this R session as the data.\n")
      dat<-get(object@name.data)
    } else {
      if(verbose) cat("Reading data from file",object@name.data,"\n")
      dat<-try(read.table(object@name.data,header=header,sep=sep,na=na.strings))
      if(class(dat)=="try-error") stop("The file ",object@name.data," does not exist. Please check the name and path.\n")      
      if(verbose) {
      	cat("These are the first lines of the dataset as read into R. Please check the format of the data is appropriate, if not, modify the na and/or sep items and retry:\n")
      	print(head(dat))
      }
    }
    if(dim(dat)[2]<2) {
      cat("The dataset contains only one column. To compute npde, we need at least 3 columns, with subject ID, predictor (at least one) and response. \nPlease check the field separator, currently given as:", paste("sep=\"",sep,"\"",sep=""), "\n")
      return("Creation of npdeData failed")
    }
# Automatic recognition of columns 
#    ID (one of id, subject or sujet regardless of case)
#    response (one of Y, conc, concentration, resp, response regardless of case)
#    predictors (time and/or dose, regardless of case)
# ECO TODO: improve automatic recognition ?
# check that we have at least a column id, response and X
    if(!is.na(as.integer(object@name.group))) {
# group given as a column number
      object@name.group<-colnames(dat)[as.integer(object@name.group)]
    }
    if(is.na(object@name.group) || object@name.group=="") {
    	if(!detect) {
    		cat("Missing ID column and automatic detection is OFF. Please provide a valid name for the ID column\n")
    		return("Creation of npdeData failed")
    	}
      if(verbose) cat("Missing ID column, attempting to detect it\n")
      object@name.group<-""
      i1<-match("id",tolower(colnames(dat)))
      if(length(i1)==0 | is.na(i1)) {
      	i1<-c(match(c("subject","sujet","group"),tolower(colnames(dat))))
      }
      if(length(i1)>0) {
        object@name.group<-colnames(dat)[i1[1]]
        if(verbose) cat("    no name for the group variable (ID) given, will use column --",object@name.group,"-- in the dataset.\n")
      }
    }
    if(object@name.group=="" | is.na(match(object@name.group,colnames(dat)))) {
      cat("Please provide a name for the ID column.\n")
      return("Creation of npdeData failed")
    }
# Predictors
   i1<-as.integer(object@name.predictor[!is.na(as.integer(object@name.predictor))])
    if(length(i1)>0) { 
      object@name.predictor[!is.na(as.integer(object@name.predictor))]<- colnames(dat)[i1]
    }
    if(is.na(object@name.predictor) | length(object@name.predictor)==0 | (length(object@name.predictor)==1 & object@name.predictor[1]=="")) {
    	if(!detect) {
    		cat("Missing X column and automatic detection is OFF. Please provide a valid name for the column with the predictor.\n")
    		return("Creation of npdeData failed")
    	}
    	if(verbose) cat("Missing predictor column, attempting to detect it\n")
      object@name.predictor<-""
      i1<-c(match(c("xobs","time","temps","tps","tim","x","dose"), tolower(colnames(dat))))
      i1<-i1[!is.na(i1)]
      if(length(i1)>0) {
        object@name.predictor<-colnames(dat)[i1][1]
        if(verbose) cat("    no name for the predictor variable given, will use column(s) --",object@name.predictor,"-- in the dataset.\n")
      }
    }
    id1<-match(object@name.predictor,colnames(dat),nomatch=0)
    if(length(id1[id1==0])>0) {
    	if(verbose) cat("    cannot find column(s) --",object@name.predictor[id1==0],"-- dropping them from the data.\n")
    }
    xnam<-object@name.predictor[id1>0]
    if(length(xnam)==0) object@name.predictor<-"" else object@name.predictor<-xnam
    if(length(xnam)==0) {
      cat("Please provide at least one predictor.\n")
      return("Creation of npdeData failed")
    }
# Response
    if(!is.na(as.integer(object@name.response))) { 
# response given as a column number
      object@name.response<-colnames(dat)[as.integer(object@name.response)]
    }
    if(is.na(object@name.response) || object@name.response=="") {
    	if(!detect) {
    		cat("Missing response column and automatic detection is OFF. Please provide a valid name for the column with the response.\n")
    		return("Creation of npdeData failed")
    	}
    	if(verbose) cat("Missing response column, attempting to detect it\n")
      object@name.response<-""
      i1<-match("y",tolower(colnames(dat)))
      if(length(i1)==0 | is.na(i1)) { 
        i1<-c( match(c("yobs","resp","conc"),tolower(colnames(dat))), grep("response",tolower(colnames(dat)),fixed=TRUE),grep("concentration", tolower(colnames(dat)),fixed=TRUE))
        i1<-i1[!is.na(i1)]
      }
      if(length(i1)>0) {
        object@name.response<-colnames(dat)[i1[1]]
        if(verbose) cat("    no name for the response variable given, will use column --",object@name.response,"-- in the dataset.\n")
      }
    }
    if(is.na(object@name.response)) object@name.response<-""
    if(object@name.response=="" | is.na(match(object@name.response,colnames(dat)))) {
      cat("Please provide a name for the response column.\n")
      return("Creation of npdeData failed")
    }
# ECO TODO: verifier que les colonnes existent et sinon corriger
    
# IPRED : column with individual predictions
    detect.ipred<-FALSE
    if(length(object@name.ipred)>0 && !is.na(as.integer(object@name.ipred))) # ipred given as a column number
    	object@name.ipred<-colnames(dat)[as.integer(object@name.ipred)]
    if(length(object@name.ipred)>0 && match(object@name.ipred,colnames(dat),nomatch=0)==0) {
    	if(detect & verbose) cat("Can't find a column named",object@name.ipred,"in the dataset for individual predictions, will attempt automatic detection.\n")
    	object@name.ipred<-character()
    }
    if(length(object@name.ipred)==0 || is.na(object@name.ipred)) detect.ipred<-TRUE
    if(detect.ipred) {
    	i1<-c(grep("ipred",tolower(colnames(dat)),fixed=T))
    	if(length(i1)>0) {
    		object@name.ipred<-colnames(dat)[i1[1]]
    		if(detect.ipred & verbose) cat("    assuming that individual predictions are given in column --",object@name.ipred,"-- in the dataset (to ignore this column, add the argument detect=FALSE in the call to npdeData()).\n")
    	}
    }
# CENS : column indicating censoring
    detect.cens<-FALSE
    if(length(object@name.cens)>0 && !is.na(as.integer(object@name.cens))) # cens given as a column number
    	object@name.cens<-colnames(dat)[as.integer(object@name.cens)]
    if(length(object@name.cens)>0 && match(object@name.cens,colnames(dat),nomatch=0)==0) {
    	if(detect & verbose) cat("Can't find a column named",object@name.cens,"in the dataset containing censoring, will attempt automatic detection.\n")
    	object@name.cens<-character()
    }
    if(length(object@name.cens)==0 || is.na(object@name.cens)) detect.cens<-TRUE
    if(detect.cens) {
    	i1<-c(grep("cens",tolower(colnames(dat)),fixed=T))
    	if(length(i1)>0) {
    		object@name.cens<-colnames(dat)[i1[1]]
    		if(detect.cens & verbose) cat("    assuming that censoring information is given in column --",object@name.cens,"-- in the dataset (to ignore this column, add the argument detect=FALSE in the call to npdeData()).\n")
    	}
    }
    if(length(object@name.cens)>0) { # checking validity of censoring column
    	if(!isTRUE(all.equal(sort(unique(dat[,object@name.cens]), na.last=TRUE),as.integer(c(0,1))))) {
    		if(verbose) cat("The column with censoring information should only contain 0 and 1s.\n")
    	  object@name.cens<-character()
    }}
# MDV : column indicating missing data
    detect.miss<-FALSE
    if(length(object@name.miss)>0 && !is.na(as.integer(object@name.miss))) # miss given as a column number
    	object@name.miss<-colnames(dat)[as.integer(object@name.miss)]
    if(length(object@name.miss)>0 && match(object@name.miss,colnames(dat),nomatch=0)==0) {
    	if(detect & verbose) cat("Can't find a column named",object@name.miss,"in the dataset containing missing data status, will attempt automatic detection.\n")
    	object@name.miss<-character()
    }
    if(length(object@name.miss)==0 || is.na(object@name.miss)) detect.miss<-TRUE
    if(detect.miss) {
    	i1<-c(grep("mdv",tolower(colnames(dat)),fixed=T), grep("miss",tolower(colnames(dat)),fixed=T))
    	if(length(i1)>0) {
    		object@name.miss<-colnames(dat)[i1[1]]
    		if(detect.miss & verbose) cat("    assuming that column --",object@name.miss,"-- in the dataset contains missing data information (to ignore this column, add the argument detect=FALSE in the call to npdeData()).\n")
    	}
    }
    if(length(object@name.miss)>0) { # checking validity of MDV column
    	if(!isTRUE(all.equal(sort(unique(dat[,object@name.miss]), na.last=TRUE),as.integer(c(0,1)))) & !isTRUE(all.equal(sort(unique(dat[,object@name.miss]), na.last=TRUE),as.integer(c(0)))) & !isTRUE(all.equal(sort(unique(dat[,object@name.miss]), na.last=TRUE),as.integer(c(1))))) {
    		if(verbose) cat("The column with information about missing data should only contain 0 and 1s.\n")
    		object@name.miss<-character()
    	}}
# Covariates
    if(length(object@name.covariates)>0 & object@name.covariates[1]!="") {
      i1<- as.integer(object@name.covariates[!is.na(as.integer(object@name.covariates))])
      object@name.covariates[!is.na(as.integer(object@name.covariates))]<- colnames(dat)[i1]
    }
# Saving covariates in the original format in ocov, transforming binary covariates in dat to factors - No: here we can keep covariates as factors
#     [email protected]<-dat[,[email protected],drop=FALSE]
#     for(icov in [email protected]) {
#     	if(length(unique(dat[,icov]))==2) dat[,icov]<-as.integer(factor(dat[,icov]))-1
#     }   
    
    if(nchar(object@name.group)*length(object@name.predictor)* nchar(object@name.response)<=0) {
      stop("Please check the structure of the data file and provide information concerning which columns specify the group structure (ID), the predictors (eg dose, time) and the response (eg Y, conc). See documentation for automatic recognition of column names for these elements.\n")
    }
# Data
    all.names<-c(object@name.group,object@name.predictor,object@name.response, object@name.covariates,object@name.miss,object@name.cens,object@name.ipred)
    tab<-dat[,all.names]
# Index (ID may be numbers, strings,...)
    id<-tab[,object@name.group]
    # ECO TODO: et si un sujet n'a que des donnees NA ???
    object@N<-length(unique(id))
    nind.obs.full<-tapply(id,id,length) # individual numbers of observations (1xN)
    nind.obs.full<-nind.obs.full[match(unique(id),names(nind.obs.full))]
    tab<-data.frame(index=rep(1:object@N,times=nind.obs.full),tab)
# Missing data
    if(length(object@name.miss)>0) mdv<-tab[,object@name.miss] else {
    	mdv<-rep(0,length(id))
    	object@name.miss<-"mdv"
    } 
    mdv[is.na(tab[,object@name.response])]<-1
    tab[,object@name.miss]<-mdv    
    object@data<-tab
    object@ind<-which(mdv==0)
    icens<-numeric()
    if(length(object@name.cens)>0) {
    	icens<-which(mdv==0 & dat[,object@name.cens]==1)
    } 
    object@icens<-icens
    object@not.miss<-(mdv==0)
# ECO TODO: what about missing data in covariates & predictor columns
    if(length(object@name.covariates)>0 && sum(is.na(object@data[object@not.miss,object@name.covariates]))>0) {
    	tab<-object@data
    	for(icov in object@name.covariates) {
    		for(i in 2:dim(tab)) {
    			if(is.na(tab[i,icov])) tab[i,icov]<-tab[(i-1),icov]    			
    		}    		
    	}
    	object@data<-tab
    }
#    for(i in [email protected]) 
#      [email protected][is.na([email protected][,i]),[email protected]]<-1
    tb1<-tab[tab[,object@name.miss]==0,]
    id1<-tb1[,1]
    object@ntot.obs<-dim(tb1)[1] # total number of observations
    nind.obs<-tapply(id1,id1,length) # individual numbers of observations (1xN)
    nind.obs<-nind.obs[match(unique(id1),names(nind.obs))]
    object@nind.obs<-c(nind.obs)    

#    [email protected]<-list([email protected],[email protected], [email protected], [email protected])
    options(ow) # reset
    validObject(object)
    return(object)
  }
)

#' Read simulated data
#' 
#' Creates an NpdeSimData object containing the simulated data, either from disk or from a dataframe
#' 
#' @param header boolean indicating whether the file has a header (mandatory if 
#' detect is TRUE)
#' @param sep field separator (for files on disk)
#' @param na.strings strings to be considered as indicating NA
#' @return an object of class NpdeSimData
#' @keywords internal

setMethod("read.npdeSimData","NpdeSimData",
	function(object,header,sep,na.strings,verbose) {
    ow <- options("warn")
    options("warn"=-1)
    if(exists(object@name.simdata)) {
    	if(verbose) cat("Using the object called",object@name.simdata,"in this R session as the data.\n")
      dat<-get(object@name.simdata)
    } else {
    	if(verbose) cat("Reading data from file",object@name.simdata,"\n")
      if(missing(header)) header<-TRUE
      if(missing(sep)) sep<-""
      if(missing(na.strings)) na.strings<-c(".","NA")
      dat<-try(read.table(object@name.simdata,header=header,sep=sep,na=na.strings))
      if(class(dat)=="try-error") stop("The file ",object@name.simdata," does not exist. Please check the name and path.\n")      

      x1<-unlist(strsplit(as.character(dat[1,1]),",",fixed=TRUE))
      mysep<-""
      if(length(x1)>1) mysep<-","
      x1<-unlist(strsplit(as.character(dat[1,1]),";",fixed=TRUE))
      if(length(x1)>1) mysep<-";"
      if(mysep!="") dat<-read.table(object@name.simdata,na=c(".","NA"),sep=mysep)
      if(!is.numeric(dat[1,1]))
        dat<-read.table(object@name.simdata,na=c(".","NA"),header=TRUE,sep=mysep)
      if(!is.numeric(dat[1,1])) {
        cat("The format of the file containing the simulated data is unknown.\n")
        cat("Please use a standard R table format, with or without header,\n")
        cat("and with one of the following separators: \n")
        cat("         TAB or space(s), commas (',') or semicolons (';')\n")
        cat("Also note that a dot should be used to indicate digits in numbers.\n")
        stop("Exiting npde\n")
      }
    	if(verbose) {
    		cat("These are the first lines of the dataset as read into R. Please check the format of the data is appropriate, if not, modify the na and/or sep items and retry:\n")
        print(head(dat))
    	}
    }
    colnames(dat)<-c("idsim","xsim","ysim")
    object@datsim<-dat

    options(ow) # reset
    validObject(object)
    return(object)
  }
)

##################################################################################
#
#' Creates a NpdeData object
#'
#' This function is used to create a NpdeData object, representing a longitudinal data structure, and fill it with data from a dataframe or a file on disk
#' 
#' @usage npdeData(name.data,header=TRUE,sep="",na=c(".","NA"),name.group, name.predictor,name.response, name.covariates,name.cens,name.miss,name.ipred, units=list(x="",y="",covariates=c()),detect=TRUE,verbose=FALSE)
#' @param name.data name of the file containing the observed data, or a dataframe
#' containing the observed data
#' @param sep field separator (for files on disk)
#' @param na strings to be considered as indicating NA
#' @param header boolean indicating whether the file has a header (mandatory if 
#' detect is TRUE)
#' @param name.group name/number of the column in the observed data containing the 
#' patient ID (if missing and detect is TRUE, columns named id, subject or sujet 
#' (regardless of case) will be assumed to contain this information)
#' @param name.predictor name/number of the column in the observed data containing 
#' the independent variable X (if missing and detect is TRUE, columns named xobs, 
#' time, dose, x, temps, tim (regardless of case) will be assumed to 
#' contain this information)
#' @param name.response name/number of the column in the observed data containing 
#' the dependent variable Y (if missing and detect is TRUE, columns named yobs, 
#' response, resp, conc, concentration (regardless of case) will be assumed to 
#' contain this information)
#' @param name.miss name/number of the column containing information about missing 
#' data (MDV) (if missing and detect is TRUE, column called mdv or miss 
#' (regardless of case) will be assumed to contain this information)
#' @param name.cens name/number of the column containing information about censored 
#' data (cens) (if missing and detect is TRUE, column with a name containing cens 
#' (regardless of case) will be assumed to contain this information)
#' @param name.covariates name/number of the column(s) containing covariate 
#' information (optional)
#' @param name.ipred name/number of the column(s) with individual predictions
#' (ipred)  (if missing and detect is TRUE, column with a name containing ipred 
#' (regardless of case) will be assumed to contain this information)
#' @param units a list with components x, y and cov (optional), specifying the
#' units respectively for the predictor (x), the response (y), and the covariates 
#' (a vector of length equal to the number of covariates). Units will default to (-) if not given.	
#' @param detect a boolean controlling whether automatic recognition of columns in the dataset is on, defaults to TRUE
#' @param verbose whether to print warning messages, defaults to FALSE (set to TRUE to check how data is being handled)
#' 
#' @return an object of class NpdeData
#' @author Emmanuelle Comets <[email protected]@bichat.inserm.fr>
#' @seealso \code{\link{npde}}, \code{\link{autonpde}}
#' @references K. Brendel, E. Comets, C. Laffont, C. Laveille, and F.
#' Mentre. Metrics for external model evaluation with an application to the
#' population pharmacokinetics of gliclazide. \emph{Pharmaceutical Research},
#' 23:2036--49, 2006.
#' @keywords models
#' @export
#' @examples
#' 
#' data(theopp)
#' 
#' x<-npdeData(theopp) # Automatic detection
#' print(x)
#' x<-npdeData(theopp,name.group="ID",name.predictor="Time",name.response="Conc", 
#' name.covariates=c("Wt"),units=list(x="hr",y="mg/L",covariates="kg")) # Explicit
#' print(x)
#' plot(x)

npdeData<-function(name.data,header=TRUE,sep="",na=c(".","NA"),name.group, name.predictor,name.response, name.covariates,name.cens,name.miss,name.ipred, units=list(x="",y="",covariates=c()),detect=TRUE,verbose=FALSE) {
# setting proper types for the NpdeData class
  if(missing(name.data) ||length(name.data)==0) {
    cat("Error in npdeData: please provide the name of the datafile or dataframe (between quotes)\n")
    return("Creation of NpdeData failed")
  }
  if(is.data.frame(name.data)) name.data<-deparse(substitute(name.data))
  if(missing(name.group)) name.group<-"" else name.group<-as.character(name.group)
  if(missing(name.predictor)) name.predictor<-"" else name.predictor<-as.character(name.predictor)
  if(missing(name.response)) name.response<-"" else  name.response<-as.character(name.response)
  if(missing(name.covariates) || name.covariates[1]==0) name.covariates<-character() else name.covariates<-as.character(name.covariates)
  if(missing(name.miss) || name.miss==0) name.miss<-character() else name.miss<-as.character(name.miss)
  if(missing(name.cens) || name.cens==0) name.cens<-character() else name.cens<-as.character(name.cens)
  if(missing(name.ipred) || name.ipred==0) name.ipred<-character() else name.ipred<-as.character(name.ipred)
  if(missing(detect)) detect<-TRUE
  x<-new(Class="NpdeData",name.data=name.data,name.group=name.group, name.predictor=name.predictor,name.response=name.response, name.covariates=name.covariates,name.cens=name.cens,name.miss=name.miss, name.ipred=name.ipred,units=units)
#  showall(x)
  if(detect & verbose) cat("Automatic detection of variables is ON. The program will attempt to detect both mandatory variables (ID, X, Y) and optional variables (IPRED, MDV, CENS) when they are not specifically given or when the user-specified names are not found in the dataset, by looking in the names of the columns (to override this behaviour, please use argument detect=FALSE in the call to npdeData().\n")
  x1<-read.npdeData(x,header=header,sep=sep,na=na,detect=detect, verbose=verbose)
  if(class(x1)!="character") {
  if(length(x1["name.cens"])==0) loq<-NA else {
  if(sum(x1["data"][x1["data"][,x1["name.miss"]]==0,x1["name.cens"]])>0) {
    yloq<-x1["data"][x1["data"][,x1["name.cens"]]==1 & x1["data"][,x1["name.miss"]]==0,x1["name.response"]]
    if(length(unique(yloq))==1) {
      loq<-unique(yloq)
      if(verbose) cat("Same LOQ for all missing data, loq=",loq,"\n")
    } else {
      loq<-min(unique(yloq),na.rm=TRUE)
      if(verbose)cat("There are different LOQ for different observations, setting loq to the lowest value of",loq,"\n")
    }
    }
    x1["loq"]<-loq
  }
  if(verbose) {
  	cat("\n\nThe following NpdeData object was successfully created:\n\n")
  	print(x1,nlines=0)
  }
  } else x1<-"Creation of NpdeData failed"
  return(x1)
}

npdeSimData<-function(npde.data,name.simdata,verbose=FALSE) {
	if(is.data.frame(name.simdata)) name.simdata<-deparse(substitute(name.simdata))
	ierror<-FALSE
	if(missing(npde.data)) {
		ierror<-TRUE
		cat("   Error: Missing first argument.\n")
	}
	if(!ierror) {
	  x1<-try(class(npde.data))
	  cat(x1)
	  if(class(x1)=="try-error") {
	  	ierror<-TRUE
	  	cat("   Error:", deparse(substitute(npde.data)),"does not exist.\n")
	  }
	  if(!ierror && x1!="NpdeData") {
	  	ierror<-TRUE 
	  	cat("   Error:", deparse(substitute(npde.data)),"is not a NpdeData object.\n")
	  }
	}
	if(ierror) {
			cat("Function npdeSimData requires two mandatory arguments: first, a NpdeData object created by a call to npdeData() (see help page for the syntax of that function), and the name of a matching dataset containing the simulated data (either a file on disk or a data.frame. Please refer to the documentation for details and examples.\n")
		return("Creation of NpdeSimData failed")
	}
	x1<-new(Class="NpdeSimData",name.simdata=name.simdata)
	x<-read.npdeSimData(x1,verbose=verbose)
	if(sum(npde.data["data"][,npde.data["name.miss"]])>0) {
		if(verbose) cat("There are rows with MDV=1 in the original dataset, the corresponding rows will be removed from the simulated dataset.\n")
	}
	nrep<-dim(x@datsim)[1]/dim(npde.data@data)[1]
	x@nrep<-as.integer(nrep)
	if(nrep<1000 & verbose) {
		cat("Warning: the number of simulations is",nrep,"which may be too small.\n")
		cat("We advise performing at least 1000 simulations to compute npde.\n")
	} 
	irsim<-rep(1:nrep,each=dim(npde.data@data)[1])
	x@datsim$irsim<-irsim
	return(x)
}    
##################################################################################

Try the npde package in your browser

Any scripts or data that you put into this service are public.

npde documentation built on May 2, 2019, 5:23 p.m.