Compare Duration Models on a Specific Data Set

Share:

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

1
2
3
4
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 christophe.pouzat@gmail.com

References

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

See Also

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

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
## 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)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.