compModels: Compare Duration Models on a Specific Data Set

In STAR: Spike Train Analysis with R

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?

Details

Fits are performed by maximizing the likelihood.

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.

Author(s)

Christophe Pouzat [e-mail]

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)

Example output

Loading required package: survival
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
Loading required package: R2HTML
Loading required package: gss
Loading required package: codetools
A spike train with 195 events, starting at: 0.318 and ending at: 30.379 (s).
The mean ISI is: 0.155 and its SD is: 0.283 (s).
The mean log(ISI) is: -2.838 and its SD is: 1.291
The shortest interval is: 0.007
 and the longest is: 1.99 (s).
 invgauss     lnorm    llogis   weibull     gamma      rexp 
-477.1587 -448.4228 -444.7260 -387.0661 -364.5763 -348.6731 
Warning messages:
1: In rep(yi[si > 0], each = ni[si > 0]) :
  first element used of 'each' argument
2: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
3: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
4: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
5: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
6: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
7: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
 invgauss     lnorm    llogis     gamma   weibull      rexp 
-394.7610 -394.6475 -392.1374 -388.5818 -372.4361 -372.3161 
Warning messages:
1: In rep(yi[si > 0], each = ni[si > 0]) :
  first element used of 'each' argument
2: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
3: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
4: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
5: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
6: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
7: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
 invgauss     lnorm    llogis     gamma   weibull      rexp 
-394.7610 -394.6475 -392.1374 -388.5818 -372.4361 -372.3161 
Warning messages:
1: In rep(yi[si > 0], each = ni[si > 0]) :
  first element used of 'each' argument
2: In rep(yi[ni > 0], each = ni[ni > 0]) :
  first element used of 'each' argument
 invgauss     lnorm     gamma    llogis   weibull      rexp 
-169.4718 -169.2137 -168.0349 -167.1806 -163.5449 -146.0264 
Warning messages:
1: In rep(yi[si > 0], each = ni[si > 0]) :
  first element used of 'each' argument
2: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
3: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
4: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
5: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
6: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
7: In rep(yi[si > 0], each = si[si > 0]) :
  first element used of 'each' argument
 invgauss     lnorm     gamma    llogis   weibull      rexp 
-169.4718 -169.2137 -168.0349 -167.1806 -163.5449 -146.0264 
Warning messages:
1: In rep(yi[si > 0], each = ni[si > 0]) :
  first element used of 'each' argument
2: In rep(yi[ni > 0], each = ni[ni > 0]) :
  first element used of 'each' argument

STAR documentation built on May 2, 2019, 4:53 p.m.