Bayesian inference of grid-cell firing patterns using hierarchical Poisson point process models, based on arXiv 2303.17217 preprint:
Papastathopoulos, I., Auld, G., Lindgren, F., Gerlei, K. and Nolan, N, (2023). ‘Bayesian inference of grid cell firing patterns using Poisson point process models with latent oscillatory Gaussian random fields’.
- load_data.R
spde.osc.temporalcontains the main code for processing data into a pipeline for computing integration weights and fiting Poisson point process models with latent Gaussian effects (M0, M1 and M2). Models M0 and M1 may be fit fast on a standard computer. Model M2 demands memory.predictions_M0.M1.M2.Rcontains code used in plotting summaries of the posterior of intensity functions.
Packages and helper functions, and setting of options for optimization
in INLA software.
library(tidyverse)
library(dplyr)
library(purrr)
library(sp)
library(INLA)
library(inlabru)
source("R/Functions.R")We load the raw data (session id M14_2018-05-16_11-29-05, gridness score 0.913)
load(url("https://www.maths.ed.ac.uk/~ipapasta/mouse_data.RData"))
ls()The firing events are stored in data frame Y and the positional data in data frame X.
##
## Firing events and trajectory
##
mycoords <- SpatialPoints(Y[,c("position_x", "position_y")])
Pl <- Polygon(cbind(X$position_x, X$position_y)); ID <- "[0,1]x[0,1]"
Pls <- Polygons(list(Pl), ID=ID); SPls <- SpatialPolygons(list(Pls))
trajectory <- SpatialPolygonsDataFrame(SPls, data.frame(value=1, row.names=ID))
## plot(trajectory)
points(mycoords, col=2, pch=16, cex=0.5)
Next we create a spatial and and a circular mesh. The parameter k in
the code below controls the resolution of the spatial mesh. Smaller
values for k yield finer spatial meshes. The spatial mesh is
constructed in such a way that the resolution is finer for regions
where we observe more spikes and coarser otherwise. Also, since
oscillating Gaussian random fields exhibit strong resonance effects,
it is important to triangulate the domain outside the arena where the
animal runs. The properties of the mesh in the region that is exterior
to the arena can be taken to coarser. For the circular mesh, we choose
a regular grid of of 29 vertices on the circle (the 30th is folded in
the 1st due to cyclicity in the construction below, see mesh.hd$n).
k <- 5
mesh <- inla.mesh.2d(loc=mycoords, max.edge=c(k, 25*k), offset=c(0.03, 120), cutoff=k/2)
p <- mesh$n
p.theta <- 30
theta.nodes <- seq(0, 2*pi, len=p.theta)
mesh.hd <- inla.mesh.1d(theta.nodes, boundary="cyclic", degree=1)We plot the spatial and circular meshes obtained from the code above.
plot(mesh, asp=1)
x1 <- cos(mesh.hd$loc); y1 <- sin(mesh.hd$loc)
plot(x1, y1, pch=16, asp=1, axes=FALSE, xlab="", ylab="", main="circular mesh")
abline(h=0); abline(v=0)
text(-.05,1, paste("1")); text(1,-.05, paste("1"))Illustrated next are the meshes used for the constructions of the latent spatial and latent head-directional processes
To facilitate a stable numerical integration for the integral in the Poisson point process likelihood, all line segments from the trajectory need to be further split so that each segment falls precisely
- in one and only one triangle of the spatial mesh and;
- in one and only one arc out of the possible 29 arcs of the circular mesh.
The code below uses the wrapper function
split.segments.wrapper.function. First, this function uses
split.arcs to split the path segments so that each new segment is
boxed by an arc on the circular mesh. Then, the function
inlabru::split_lines is used to split again the segments obtained
from of split.arcs so that every line segment is also boxed by a triangle.
Ypos.ls <- split.segments.wrapper.function(X=X, mesh=mesh, mesh.hd=mesh.hd)
Ypos <- Ypos.tmp.ls$Ypos
filter.index <- Ypos.tmp.ls$filter.index
line.segments <- Ypos.ls$line.segmentsThis output is stored in Ypos and the plots below illustrate the
effect of line splitting. The new line segments are plotted over a
subregion of the spatial domain on the left panel below. Each newly
derived line segment is given by the start and end of an arrow whilst
the raw segments by contiguous blue dots. Similarly for the head
direction on the right panel below, each newly derived line segment is
given by contiguous black circles whilst the raw segments by blue
dots. If an line segment was initially boxed by a triangle and by an
arc segment, then this line segment remains unchanged. When there is a
split, however, the values of the covariate (location, head direction
and time) are imputed with linear interpolation.
par(mfrow=c(1,2))
plot(mesh, xlim=c(46, 54), ylim=c(46,54), asp=1, main="")
arrows(do.call("rbind", Ypos$sp)[,1], do.call("rbind", Ypos$sp)[,2],
do.call("rbind", Ypos$ep)[,1], do.call("rbind", Ypos$ep)[,2], col=2, lwd=1, length=0.05)
points(X$position_x, X$position_y, col="blue", pch=16, cex=0.5)
plot(Ypos$time[1:50,], Ypos$hd[1:50], type="b", cex=1, xlab="time", ylab="head direction")
points(X$synced_time, X$hd, col="blue", cex=0.5, pch=16)
abline(h = seq(0, 2*pi, len=30), lty=2, lwd=.5)
Lastly, we save the covariate in objects coords.trap, HD.data and T.data.
coords.trap <- rbind(do.call("rbind",Ypos$sp)[filter.index,], tail(do.call("rbind",Ypos$ep),1))
HD.data <- c(do.call("c", (Ypos %>% mutate(HD=lapply(HDi, function(x) attr(x, "data"))))$HD), tail(Ypos$hd.lead, 1))
T.data <- c(do.call("c", (Ypos %>% mutate(T=lapply(Ti, function(x) attr(x, "data"))))$T), tail(Ypos$time.lead, 1))Additionally to the spatial and circular meshes, a temporal mesh is
also required when the effect of time on the variation of spikes is
included in models. The temporal mesh is constructed by thinning the
grid of times that are obtained from the newly derived line
segments. Below we choose the fixed value 300 for thinning so that,
if the consecutive times at the starting positions of the line segment
(obtained via splitting) are denoted by t[1],t[2], .., then the
temporal mesh is taken as t[1],t[300],t[600] ... With this
construction, each line segment is also boxed by contiguous time
points in the temporal mesh.
par(mfrow=c(1,2))
plot(trajectory)
points(coords.trap[seq(1, length(T.data), by = 300), 1], coords.trap[seq(1, length(T.data), by = 300), 2], pch=4, col=2, cex=.8, asp=1)
plot(mesh1d$loc[1:50], rep(0,50), pch=16, cex=0.6, xlab="time", ylab="", axes=FALSE); axis(1)
abline(h=0)The temporal mesh is shown by red crosses on the path in the spatial domain in the left panel of the Figure below and by black dots on the time axis in the right panel.
Aosc <- inla.mesh.projector(mesh, loc=coords.trap)$proj$A
Ahd <- inla.mesh.projector(mesh.hd, loc=HD.data)$proj$A
A <- inla.row.kron(Ahd, Aosc)
Atilde <- inla.mesh.projector(mesh1d, loc=T.data)$proj$AAosc.obs <- inla.spde.make.A(mesh=mesh, loc=as.matrix(data$Y %>% dplyr:: select(position_x, position_y)))
Ahd.obs <- inla.spde.make.A(mesh=mesh.hd, data$Y$hd)
Aobs <- inla.row.kron(Ahd.obs, Aosc.obs)
Atildeobs <- inla.spde.make.A(mesh=mesh1d, data$Y$firing_times)par(mfrow=c(1,3))
image(Aosc, lwd=2, asp=1)
image(A, lwd=2, asp=1)
image(Atilde, lwd=2, asp=1)First, for each line segment of the path, we need to know how much distance was traveled and how much time it took for the animal to traverse the segment.
dGamma <- c(do.call("c", Ypos$Li))
dT <- diff(T.data)Second, to compute the integration weights, we need to store all basis function evaluations (i.e., spatial, spatial-directional and temporal basis functions) at starting points of each line segment of the path.
## spatial
Aosctmp <- as(Aosc, "dgTMatrix")
Aosc.indices <- cbind(cbind(Aosctmp@i+1, Aosctmp@j+1), Aosctmp@x) # (i,j, A[i,j]) for which A[i,j] is non-zero (Omega x Theta)
Aosc.indices <- Aosc.indices[order(Aosc.indices[,1]),] %>% as.data.frame #
## spatial-directional
Atmp <- as(A, "dgTMatrix")
A.indices <- cbind(cbind(Atmp@i+1, Atmp@j+1), Atmp@x) # (i,j, A[i,j]) for which A[i,j] is non-zero (Omega x Theta)
A.indices <- A.indices[order(A.indices[,1]),] %>% as.data.frame #
## temporal
Attmp <- as(Atilde, "dgTMatrix")
At.indices <- cbind(cbind(Attmp@i+1, Attmp@j+1), Attmp@x) # (i,j, A[i,j]) for which Atilde[i,j] is non-zero (Time)
At.indices <- At.indices[order(At.indices[,1]),] %>% as.data.frame
##
names(Aosc.indices) <- c("tk", "i", "psi.o") #ot: omega
names(A.indices) <- c("tk", "i", "psi.ot") #ot: omega x theta
names(At.indices) <- c("tk", "l", "psi.t")A.indices and At.indices: first column is renamed to tk which
stands for the index of the line/time/arc segment. For example, for
the spatio-directional basis functions, each tk appears 6 times,
i.e., length(A.indices[,1])/6 = N, where N is the number of
line/time/arc segments. In A.indices the second column is renamed to
i which stands for the index of the spatio-directional basis
function. In At.indices the second row is renamed to l which
stands for the index of the temporal basis function
We check dim(A)[1] = dim(Atilde)[1]= is TRUE, both matrices are
basis function evaluations at the starting coordinates and head
directions (A) and the starting times (Atilde) of the line
segments, that is, each row stores basis function evaluations for a
line. For example, for the matrix of spatio-directional basis
functions, for each starting point of a line segment, there are 6
spatio-temporal basis functions that give a non-zero contribution,
that is 3 knots of a triangle * 2 knots of an arc whilst for the
matrix of temporal basis function, there are 2 temporal basis
functions that give non-zero contributions, that is, 2 time interval
knots.
Below we use helper functions df.prism.M0.wrapper and
df.prism.M0.M1.wrapper to compute quantities that are necessary for
the integration weights. For example, the helper function
df.prism.M1.M2.wrapper used to get df.prism.M1_M2 first groups
At.indices and A.indices by line segment and then nests them so
each row of the nested data frame contains all basis function
evaluation data each line segment. Once data are nested, information
on the index of the basis functions and its associated value is stored
in new column variables named as data.x for the temporal basis
functions, and as data.y for the spatio.directional basis functions.
Information on the times, head directions, and coordinates is also
appended to the nested data frame, that is, for every line segment
(indexed by variable tk). Information on lags and leads is also
included because these are required to compute the integration weights
based on the trapezoidal rule (details will be added in an Appendix of
the statistical version of the paper). For the computation of the
weights, the lengths of the line segments (dGamma) together with
their lags and leads are also appended. Finally, in the last column
a variable named val. Fix a line segment, say the first one tk=1.
Then, for example, the first elements of the column variables data.x
and data.y (these are lists due to the nest operation) are:
> data.y[[1]]
# A tibble: 6 x 2
i psi.ot
<dbl> <dbl>
1 7660 0.0475
2 7726 0.405
3 8037 0.246
4 8932 0.0205
5 8998 0.175
6 9309 0.106
> data.x[[1]]
l psi.t
<dbl> <dbl>
1 1 1
2 2 0The code creates the Cartesian product {1,2} X {7660, 7726, 8037,
8932, 8998, 9309}, that is, the set of all ordered pairs (a,b) where
a is in {1,2} and b in {7660, 7726, 8037, 8932, 8998, 9309} with
expand.grid, and calculates, for each pair, the product of the basis
functions \psi_{T} * \psi_{Omega x Theta}. Lastly, the function
returns a data frame that has the index of the temporal basis function
l, the index of the spatio-directional basis function i, and the
product of the basis functions val. This data framed is stored in a
column variable named val. The final commands discard data.x and
data.y which are no longer used and unnests the data frame to bring
it back in standard form
df.prism.M0 <- df.prism.M0.wrapper(Aosc.indices = Aosc.indices, dGamma=dGamma, T.data=T.data, HD.data=HD.data,
coords.trap=coords.trap) %>% unnest(cols=c(val.M0))
df.prism.M1_M2 <- df.prism.M1.M2.wrapper(At.indices= At.indices, A.indices=A.indices, dGamma=dGamma, T.data=T.data, HD.data=HD.data, coords.trap=coords.trap)
df.prism.M1 <- df.prism.M1_M2 %>% dplyr::select(-val.M2) %>% unnest(cols=c(val.M1))
df.prism.M2 <- df.prism.M1_M2 %>% dplyr::select(-val.M1) %>% unnest(cols=c(val.M2))df.W.M0 <- rbind(df.prism.M0 %>% mutate(group=tk, dGamma.lag=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, i, val.M0),
df.prism.M0 %>%
filter(tk!=1) %>%
mutate(time=time.lag, direction=direction.lag, coords=coords.lag,
group=tk-1,
dGamma=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, i, val.M0)) %>%
arrange(group) %>%
mutate(dGamma.trap = dGamma + dGamma.lag)
tol <- 0
df.dGamma.sum.k.kplus1.M0 <- df.W.M0 %>% group_by(group, i) %>%
summarize(val = sum(max(dGamma.trap*val.M0, tol))/2,
time = unique(time),
direction=unique(direction),
coords=unique(coords)) %>%
ungroup %>% group_by(i) %>%
summarize(val = sum(val))
W.M0 <- sparseVector(i=df.dGamma.sum.k.kplus1.M0$i,
x=df.dGamma.sum.k.kplus1.M0$val,
length=mesh$n)
W.M0.vector <- sparseVector(i=df.dGamma.sum.k.kplus1.M0$i,
x=df.dGamma.sum.k.kplus1.M0$val,
length=mesh$n)
W.ipoints.M0 <- as(W.M0, "sparseMatrix")
W.ipoints.M0 <- data.frame(coords.x1 = mesh$loc[W.ipoints.M0@i+1,1],
coords.x2 = mesh$loc[W.ipoints.M0@i+1,2],
weight=W.ipoints.M0@x) df.W.M1 <- rbind(df.prism.M1 %>% mutate(group=tk, dGamma.lag=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, i, val.M1),
df.prism.M1 %>%
filter(tk!=1) %>%
mutate(time=time.lag, direction=direction.lag, coords=coords.lag,
group=tk-1,
dGamma=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, i, val.M1)) %>%
arrange(group) %>%
mutate(dGamma.trap = dGamma + dGamma.lag)
tol <- 0
df.dGamma.sum.k.kplus1.M1 <- df.W.M1 %>% group_by(group, i) %>%
summarize(val = sum(max(dGamma.trap*val.M1, tol))/2,
time = unique(time),
direction=unique(direction),
coords=unique(coords)) %>%
ungroup %>% group_by(i) %>%
summarize(val = sum(val))
W.M1 <- sparseVector(i=df.dGamma.sum.k.kplus1.M1$i,
x=df.dGamma.sum.k.kplus1.M1$val,
length=mesh$n * mesh.hd$n)
W.M1.vector <- sparseVector(i=df.dGamma.sum.k.kplus1.M1$i,
x=df.dGamma.sum.k.kplus1.M1$val,
length=mesh$n * mesh.hd$n)
df.indices <- data.frame(dir = sort(rep(1:mesh.hd$n, mesh$n)), space = rep(1:mesh$n, mesh.hd$n), cross = 1:(mesh$n*mesh.hd$n))
mapindex2space.direction_index <- function(index){
f<-function(index.single){
as.numeric(df.indices[which(df.indices$cross==index.single),c("dir","space")])
}
t((Vectorize(f, vectorize.args="index.single"))(index))
}
mapindex2space.direction_basis <- function(index){
f<-function(index.single){
o <- as.numeric(df.indices[which(df.indices$cross==index.single),c("dir","space")])
return(c(mesh.hd$loc[o[1]], mesh$loc[o[2],-3]))
}
t((Vectorize(f, vectorize.args="index.single"))(index))
}
W.ipoints.M1 <- as(W.M1, "sparseMatrix")
W.ipoints.M1 <- data.frame(hd=mapindex2space.direction_basis(W.ipoints.M1@i+1)[,1],
coords.x1 =mapindex2space.direction_basis(W.ipoints.M1@i+1)[,2],
coords.x2 =mapindex2space.direction_basis(W.ipoints.M1@i+1)[,3],
weight=W.ipoints.M1@x) df.W.M2 <- rbind(df.prism.M2 %>% mutate(group=tk, dGamma.lag=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, l, i, val.M2),
df.prism.M2 %>%
filter(tk!=1) %>%
mutate(time=time.lag, direction=direction.lag, coords=coords.lag,
group=tk-1,
dGamma=0) %>%
dplyr::select(group, time, direction, coords, dGamma, dGamma.lag, l, i, val.M2)) %>%
arrange(group) %>%
mutate(dGamma.trap = dGamma + dGamma.lag)
tol <- 0
df.dGamma.sum.k.kplus1.M2 <- df.W.M2 %>% group_by(group, l, i) %>%
summarize(val = sum(max(dGamma.trap*val.M2, tol)),
time = unique(time),
direction=unique(direction),
coords=unique(coords))
W <- sparseMatrix(i=df.dGamma.sum.k.kplus1.M2$l,
j=df.dGamma.sum.k.kplus1.M2$i,
x=df.dGamma.sum.k.kplus1.M2$val/2)
W <- W %>% cbind(Matrix(0, nrow=nrow(W), ncol=ncol(A)-ncol(W)))
W.ipoints.M2 <- as(W, "dgTMatrix")
W.ipoints.M2 <- data.frame(firing_times=mesh1d$loc[W.ipoints.M2@i+1], hd=mapindex2space.direction_basis(W.ipoints.M2@j+1)[,1],
coords.x1 =mapindex2space.direction_basis(W.ipoints.M2@j+1)[,2],
coords.x2 =mapindex2space.direction_basis(W.ipoints.M2@j+1)[,3],
weight=W.ipoints.M2@x) %>% arrange(firing_times)The following B matrices are intended to be used with
inla.spde2.generic and see spde2_implementation.pdf. there are two
possibilities for defining models: one with inla.spde2.generic and
the other with inla.rgeneric.define. The function
inla.spde2.generic provides support for Matern models (this includes
oscillating models too) whilst inla.rgeneric.define allows user to
build the model from scratch. The latter permits priors of
hyperparameters to be defined by the user.
B.phi0.matern = matrix(c(0,1,0), nrow=1)
B.phi1.matern = matrix(c(0,0,1), nrow=1)
B.phi0.oscillating = matrix(c(0,1,0,0), nrow=1)
B.phi1.oscillating = matrix(c(0,0,1,0), nrow=1)
B.phi2.oscillating = matrix(c(0,0,0,1), nrow=1)The following commands implement the finite element method and are
used to obtain the M matrices (defined in
spde2_implementation.pdf). These are used both in
inla.spde2.generic and inla.rgeneric.define
fem.mesh <- inla.mesh.fem(mesh, order = 2)
fem.mesh.hd <- inla.mesh.fem(mesh.hd, order = 2)
fem.mesh.temporal <- inla.mesh.fem(mesh1d, order = 2)
## M matrices for spatial oscillating model
M0 = fem.mesh$c0
M1 = fem.mesh$g1
M2 = fem.mesh$g2
## M matrices for temporal model
M0.temporal = fem.mesh.temporal$c0
M1.temporal = fem.mesh.temporal$g1
M2.temporal = fem.mesh.temporal$g2
## M matrices for circular/directional model
M0.hd = fem.mesh.hd$c0
M1.hd = fem.mesh.hd$g1
M2.hd = fem.mesh.hd$g2Below we assign a prior distribution to each hyperparameter of models
M0, M1 and M2. The priors are defined in the
R/rgeneric_models.R but the specification of the hyperparameters is
given externally below.
## ------------------------------------------------------
## specification of prior distribution of hyperparameters
## ------------------------------------------------------
sigma.range.spatial.oscillating <- .4
mu.range.spatial.oscillating <- 20
sigma.spatial.oscillating <- 1/2
a.par.phi.prior.spatial.oscillating <- 2
b.par.phi.prior.spatial.oscillating <- 20
## directional model
rho.directional <- 1/(2*pi)
sigma.directional <- 1
##
rho.temporal <- 1/100
sigma.temporal <- 1/3
l = -0.98
u = 1
## initial values for optimization
initial.space <- list(theta1=2.6,theta2=0.5, theta3=-1.4)
initial.direction <- list(theta4=log(pi), theta5=0)and we all custom-made built models and use them in
inla.rgeneric.define to define our models.
source("R/rgeneric_models.R")
## define models
## oscillating.rgeneric is used for M0
## space.direction.rgeneric is used for M1 and M2
## temporal.rgeneric is used for M1 and M2
space.rgeneric <- inla.rgeneric.define(oscillating.model,
M = list(M0=M0, M1=M1, M2=M2),
theta.functions = list(theta.2.phi = theta.2.phi,
theta.2.sigma = theta.2.sigma,
theta.2.rho = theta.2.rho,
l=l, u=u),
hyperpar = list(
mu.range.spatial.oscillating = mu.range.spatial.oscillating,
sigma.range.spatial.oscillating = sigma.range.spatial.oscillating,
sigma.spatial.oscillating = sigma.spatial.oscillating,
a.par.phi.prior.spatial.oscillating = a.par.phi.prior.spatial.oscillating,
b.par.phi.prior.spatial.oscillating = b.par.phi.prior.spatial.oscillating),
prior.functions = list(prior.phi_osc = prior.phi_osc),
initial.space=initial.space)
##
space.direction.rgeneric <- inla.rgeneric.define(space.direction.model,
M=list(M0.space=M0, M1.space=M1, M2.space=M2,
M0.direction=M0.hd, M1.direction=M1.hd, M2.direction=M2.hd),
theta.functions = list(theta.2.rho = theta.2.rho,
theta.2.sigma = theta.2.sigma,
theta.2.phi = theta.2.phi,
theta.2.rho.direction = theta.2.rho.direction,
theta.2.sigma.direction = theta.2.sigma.direction,
l=l, u=u),
hyperpar = list(
mu.range.spatial.oscillating = mu.range.spatial.oscillating,
sigma.range.spatial.oscillating = sigma.range.spatial.oscillating,
sigma.spatial.oscillating = sigma.spatial.oscillating,
a.par.phi.prior.spatial.oscillating = a.par.phi.prior.spatial.oscillating,
b.par.phi.prior.spatial.oscillating = b.par.phi.prior.spatial.oscillating,
rho.directional = rho.directional,
sigma.directional = sigma.directional),
prior.functions = list(prior.phi_osc = prior.phi_osc),
initial.space=initial.space,
initial.direction = initial.direction)
##
time.rgeneric <- inla.rgeneric.define(temporal.model,
M=list(M0.temporal=M0.temporal, M1.temporal=M1.temporal, M2.temporal=M2.temporal),
hyperpar = list(
rho.temporal = rho.temporal,
sigma.temporal = sigma.temporal
))Below Y.spdf is the Y data frame except that coords are encoded
as SpatialPoints. Similarly, Ypos.sldf is the Ypos data frame
except that the segments between coords and coords.lead are
encoded as SpatialLines. For inlabru we will only use Y.spdf as
we manually compute the integration points.
Y.spdf <- SpatialPointsDataFrame(coords = SpatialPoints(cbind(Y$position_x, Y$position_y)),
data = as.data.frame(Y%>%dplyr::select(-c(position_x, position_y))))
Ypos.sldf <- SpatialLinesDataFrame(sl = SpatialLines(lapply(as.list(1:nrow(Ypos)),
function(k) Lines(list(Line(cbind(c(Ypos$coords[k,1],
Ypos$coords.lead[k,1]),
c(Ypos$coords[k,2],
Ypos$coords.lead[k,2])))), ID=k))),
data = Ypos %>% dplyr::select(-c(coords, coords.lead)))
data <- list(Ypos=Ypos, Y=Y, Yspdf=Y.spdf, Ypos.sldf = Ypos.sldf)Next, we define a rectangular area using SpatialPolygons (object
SPls.Omega) which is needed for constraining the the spatial field
to integrate to zero (see next section). This rectangle corresponds to
the arena within which the animal explores space and is defined by the
coordinate wise minima and maxima of the visited positions.
Pl.Omega <- Polygon(expand.grid(c(min(X$position_x),max(X$position_x)), c(min(X$position_y),max(X$position_y)))[c(1,2,4,3),])
ID.Omega <- "Omega"
Pls.Omega <- Polygons(list(Pl.Omega), ID=ID.Omega)
SPls.Omega <- SpatialPolygons(list(Pls.Omega))
weights.domain <- ipoints(domain=mesh, samplers=SPls.Omega)
locs <- weights.domain@coords
rownames(locs) <- NULLTo ensure identifiability between the intercept and the oscillating field, we impose integral-to-zero constraints for the spatial field, that is, we enforce the spatial oscillating field in the arena integrates to 0. This is implemented via the following matrix of contraints.
A.spatial.field_constr <- inla.spde.make.A(mesh=mesh, loc=locs,
weights=weights.domain@data[,1],
block=rep(1, nrow(weights.domain@coords)))With this in place, we can fit M0 as follows
cmp.space <- firing_times ~
spde2(cbind(coords.x1, coords.x2), model=space.rgeneric, mapper=bru_mapper(mesh,indexed=TRUE),
extraconstr=list(A=as.matrix(A.spatial.field_constr,nrow=1), e=0)) + Intercept
fit.space <- lgcp(cmp.space,
data = Y.spdf,
ips = W.ipoints.M0,
domain = list(firing_times = mesh1d),
options = list( num.threads=8,verbose = TRUE, bru_max_iter=1) )Run on a computer with at least 32GB of shared
memory. computationally demanding. Specification of code below
assumes computer has 8 cores. Switch num.threads accordingly.
For this model, we impose similar constraints for the spatial field as
in M0 but for each direction knot in the circular mesh. To do so, we
construct the matrix of constraints as follows
A.spatial.field_constr_along.directions <- as.matrix(kronecker(Diagonal(mesh.hd$n),
as.matrix(inla.spde.make.A(mesh=mesh, loc=locs,
weights=weights.domain@data[,1],
block=rep(1, nrow(weights.domain@coords))), nrow=1)))With this in place, we can fit model M1 as follows:
cmp.space.direction <- firing_times ~
spde2(list(spatial=cbind(coords.x1, coords.x2), direction=hd), model=space.direction.rgeneric,
mapper=bru_mapper_multi(list(spatial=bru_mapper(mesh,indexed=TRUE), direction=bru_mapper(mesh.hd, indexed=TRUE))),
extraconstr=list(A=as.matrix(A.spatial.field_constr_along.directions,nrow=19), e=rep(0,19))) +
Intercept
fit.space.direction <- lgcp(cmp.space.direction, data = Y.spdf,
ips = W.ipoints.M1,
domain = list(firing_times = mesh1d),
options = list( num.threads=8, verbose = TRUE, bru_max_iter=1) )Run on a computer with at least 64GB of shared
memory. computationally demanding. Switch num.threads accordingly
cmp.space.direction.time <- firing_times ~
spde2(list(spatial=cbind(coords.x1, coords.x2), direction=hd), model=space.direction.rgeneric,
mapper=bru_mapper_multi(list(spatial=bru_mapper(mesh,indexed=TRUE), direction=bru_mapper(mesh.hd, indexed=TRUE)))) +
time(firing_times, model=time.rgeneric, mapper=bru_mapper(mesh1d, indexed=TRUE)) + Intercept
fit.space.direction.time <- lgcp(cmp.space.direction.time, data = as.data.frame(Y.spdf),
ips=W.ipoints.M2,
domain = list(firing_times = mesh1d),
options=list(
num.threads=8,
verbose = TRUE, bru_max_iter=1))
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent