Description Usage Arguments Details Value Author(s) References See Also Examples
Given two point processes, A and B, simultaneously observed on the real line, this function computes a p-value for B's dependence on A. The function can be parameterised to detect triggering ('A causes B'), correlation, inhibition or anti-correlation. Further details are available in the corresponding paper.
1 2 3 |
A |
a vector of event times. We condition on A and test B given A. |
B |
a second vector of event times. B is tested conditional on A. |
start |
the start of the observation period. Must be smaller than all of the event times in A and in B. |
end |
the end of the observation period. Must be larger than all of the event times in A and in B. |
F |
an increasing function, giving the cumulative intensity of B under
the null hypothesis. If |
alternative |
can be "causes" (testing for A causing B), "inhibits" (testing for A inhibiting B), "correlation" or "anticorrelation". |
method |
can be "timeout", "simes" or "fisher". These are different ways of testing the cumulative response times. See details. |
transform |
if |
usebeforeA1 |
if |
maxtau |
if method="timeout" then the user can set |
careful |
if |
Under the null hypothesis the events of B are assumed to follow a
non-homogeneous Poisson process. This has an intensity function
λ(t). The intensity is provided by the user via the input
F
. If F
is NULL
(the default), B follows a
homogeneous Poisson process under the null hypothesis. Otherwise
F
is interpreted as F(t) = \int_{\mbox{start}}^t
λ(x) dx. The tests are invariant to any affine transformation of the input function
F^*(t) = α+β F(t). Note that F=NULL
is
equivalent to F=function(t)t
, although the latter is not
recommended because the algorithm is faster if it knows the null
hypothesis is homogeneous. To choose F
the users are invited to
estimate the non-homogeneous intensity of B via their favourite method, or
try mkF
which uses R's in-built density
estimation. Note
that the non-homogeneous Poisson assumption is a bit stronger than needed, see
corresponding paper.
The procedure computes so-called cumulative response times. These are a
sequence of values between 0 and 1, one for each event time in B (or each B after
A[1]
in the case of causes/inhibition with
usebeforeA1=FALSE
), which should be ordered uniform variables
under the null hypothesis but ‘closer to zero’ under the
alternative.
To test the cumulative response times, the user has the choice between three methods: a timeout
test (see TMT.test
), Fisher's method (see
F.test
) and Simes's test (see simes.test
).
By default A and B are analysed as given. If
transform=TRUE
, F is first used to transform A and
B so that B is homogeneous, and then the transformed processes are tested
against a homogeneous null hypothesis.
For further details see corresponding paper.
An S3 object of class “htest” with slots
p |
the p-value. |
T |
the test statistic. |
data.name |
the input processes. |
t |
the cumulative response times. |
TMTval |
if method="timeout", the largest cumulative response time
estimated to have been affected by A, |
TMTi |
if method="timeout", the index of the largest cumulative
response time estimated to have been affected by |
Patrick Rubin-Delanchy <patrick.rubin-delanchy@bristol.ac.uk> and Nicholas A Heard <n.heard@imperial.ac.uk>
Patrick Rubin-Delanchy and Nicholas A Heard. “A test for dependence between two point processes on the real line”. arXiv:1408.3845.
mkF
, TMT.test
, F.test
, simes.test
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 43 44 45 46 47 48 49 50 51 52 | A = runif(20)
Braw=runif(20)
##around ten B events are caused by A
B1=c(Braw, sample(A, 10)+abs(rnorm(10, 0,.01))); B1 = B1[B1>0&B1<1]
##about ten B events are correlated to A
B2=c(Braw, sample(A, 10)+rnorm(10, 0,.01)); B2 = B2[B2>0&B2<1]
##the ten closest B events to A are deleted (anticorrelation)
d=sapply(Braw, function(b) min(abs(b-A))); B3=Braw[-order(d)[1:10]]
##The ten B events that are closest to but after an A event are deleted (inhibition).
##(Adding an A event at zero to be sure there are 10.)
A=c(0,A); d=sapply(Braw, function(b) b-max(A[A<b])); B4=Braw[-order(d)[1:10]]
alternatives=c("causes", "correlation", "anticorrelation", "inhibits")
ps=c()
for (B in list(B1, B2, B3, B4)){
for (alternative in alternatives){
p=corrtest(A,B, alternative=alternative)$p
ps = c(ps, p)
}
}
M=matrix(ps, nrow=4, ncol=4)
colnames(M) = c("causal_data", "correlated_data", "anticorrelated_data", "inhibited_data")
rownames(M) = c("causes_test", "correlation_test", "anticorrelation_test", "inhibition_test")
##should (hopefully!) see low p-values on the diagonal:
M
## and high p-values for opposite data versus test combinations,
## e.g. testing for A inhibiting B when in fact A causes B:
M["inhibition_test", "causal_data"]
##Now for a non-homogeneous example. A and B have a common daily pattern:
##their intensity is a sinusoidal curve lambda(t) = 1+sin(2*pi*t)
start=0; end=365 #A year
##the cumulative intensity is
F=function(t){
t%/%1+t%%1+(1-cos(2*pi*t%%1))/(2*pi)
}
##Dropping 365 A and B points according to F
A=sapply(runif(365), function(u){uniroot(function(x) F(x)/365-u, interval=c(0,365))$root})
##B is independent of A aside from common periodicity
B=sapply(runif(365), function(u){uniroot(function(x) F(x)/365-u, interval=c(0,365))$root})
##Bc also has some caused events
Bc=c(B, sample(A, 10)+abs(rnorm(10, 0,.001))); Bc = Bc[Bc>start&Bc<end]
##If we do not account for the common periodicity of A and B we get
##spuriously strong evidence for A causing B (using Fisher's method for speed):
corrtest(A,B, start=start, end=end, method="fisher")
##On the other hand with the correct F, the p-value is uniformly distributed:
corrtest(A,B, start=start, end=end, F=F, method="fisher")
##For reference, with the truly caused process Bc we get
corrtest(A,Bc, start=start, end=end, method="fisher")
corrtest(A,Bc, start=start, end=end, F=F, method="fisher")
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.