STAR (version 0.3-7)

compModels: Compare Duration Models on a Specific Data Set

Description

Fit duration models with the maximum likelihood method to a given duration data set. The data can be censored. The models should be among the following list: inverse Gaussian, log normal, log logistic, gamma, Weibull, refractory exponential. The Akaike information criterion (AIC) is used to produce a numerical output. Diagnostic QQ or survival plots can also be generated.

Usage

compModels(yi, ni = numeric(length(yi)) + 1, si = numeric(length(yi)) + 1, models = c("invgauss","lnorm","gamma","weibull","llogis","rexp"), type = c("qq","s"), log = TRUE, plot = TRUE)

Arguments

yi

vector of (possibly binned) observations or a spikeTrain object.

ni

vector of counts for each value of yi; default: numeric(length(yi))+1.

si

vector of counts of uncensored observations for each value of yi; default: numeric(length(yi))+1.

models

a character vector whose elements are selected among: "invgauss", "lnorm", "gamma", "weibull", "llogis", "rexp".

type

should a QQ plot ("qq") or a survival plot ("s") be generated?

log

should a log scale be used?

plot

should a plot be generated?

Value

A vector whose component are nammed according to the model used and ordered along increasing AIC values.if argument plot is set to TRUE (the default), a plot is generated as a side effect.

Details

Fits are performed by maximizing the likelihood.

References

Lindsey, J.K. (2004) The Statistical Analysis of Stochastic Processes in Time. CUP.

See Also

qqDuration, invgaussMLE, lnormMLE, llogisMLE, rexpMLE, gammaMLE, weibullMLE

Examples

## Not run: 
# ## load spontaneous data of 4 putative projection neurons
# ## simultaneously recorded from the cockroach (Periplaneta
# ## americana) antennal lobe
# data(CAL1S)
# ## convert data into spikeTrain objects
# CAL1S <- lapply(CAL1S,as.spikeTrain)
# ## look at the individual trains
# ## first the "raw" data
# CAL1S[["neuron 1"]]
# ## next some summary information
# summary(CAL1S[["neuron 1"]])
# ## next the renewal tests
# renewalTestPlot(CAL1S[["neuron 1"]])
# ## It does not look too bad so let fit simple models
# compModels(CAL1S[["neuron 1"]])
# 
# ## Simulate a sample with 100 events from an inverse Gaussian
# set.seed(1102006,"Mersenne-Twister")
# mu.true <- 0.075
# sigma2.true <- 3
# sampleSize <- 100
# sampIG <- rinvgauss(sampleSize,mu=mu.true,sigma2=sigma2.true)
# 
# ## Compare models and display QQ plot
# compModels(sampIG,type="qq")
# 
# ## Compare models and display survival plot
# compModels(sampIG,type="s")
# 
# 
# ## Generate a censored sample using an exponential distribution
# sampEXP <- rexp(sampleSize,1/(2*mu.true))
# sampIGtime <- pmin(sampIG,sampEXP)
# sampIGstatus <- as.numeric(sampIG <= sampEXP)
# ## Compare models and display QQ plot
# ## WARNING with censored data like here the QQ plot is misleading
# compModels(yi=sampIGtime,si=sampIGstatus,type="qq")
# ## Compare models and display survival plot
# compModels(yi=sampIGtime,si=sampIGstatus,type="s")
# ## End(Not run)