MA 568 Statistical Analysis of Point Process Data
Problem set 3 solution
Long Tao and Uri Eden
Department of Mathematics and Statistics, Boston university
1
Question 1
We download the hipp data.mat file, which consists of the spiking activity of two neurons recorded
simultaneously from the CA1 regions of a rat’s hippocampus as it performed a free-foraging task.
We plot the animal’s movement trajectory with its position at the spiking times for each neuron
overlaid.
From the plot above, we find that these two neurons are both tuned to the spatial location of the
animal, and both have similar tuning properties in that their place fields nearly completely overlap.
The center of both place fields occur when both the x and y coordinates are slightly negative.
1 sp ik e data <�read . csv ( ”hipp data . csv ” , header=T) 2 attach ( sp ike data ) 3 #1. a Plot movements t r a j e c t o r y with po s i t i o n at each sp i k ing t imes 4 p lo t (xN,yN, cex =.2 , c o l=’ darkgray ’ , type=’ l ’ ) 5 idx 1 = sp ikes >0; idx 2 = sp ikes2>0
6 po in t s (xN [ idx 1 ] ,yN [ idx 1 ] , c o l=’ red ’ , cex=.5)
7 po in t s (xN [ idx 2 ] ,yN [ idx 2 ] , c o l=’ blue ’ , cex=.5)
8 t i t l e ( ’ Tra j ec tory p l o t with sp i k i ng p o s i t i o n s ’ )
Question 2
Use glm command to fit the exponential linear model of equation
�(t) = exp{�
0
+ �
1
xN (t) + �2yN (t)} (1)
2
From the summary, we see all of the p values are highly significant (p < 2e� 16). We plot out the
likelihood model intensity as a function of (x, y) position. The exponential linear model captures
the fact that the firing rate increases as we move from the origin to negative values for x and y.
However, the firing rate continues to increase in this direction even though the number of spikes
falls o↵ as x and y become more negative. In order to capture this phenomenon, we will need to
include quadratic terms of x and y.
1 model1 <� glm ( sp i k e s ˜xN+yN, poisson , data )
2 summary(model1 )
3 Cal l :
4 glm ( formula = sp i k e s ˜ xN + yN, fami ly = poisson , data = data )
5 Co e f f i c i e n t s :
6 Estimate Std . Error z va lue Pr(>| z | )
7 ( I n t e r c ep t ) �3.24918 0.02760 �117.74 <2e�16 ⇤⇤⇤
8 xN �0.50623 0.04724 �10.72 <2e�16 ⇤⇤⇤
9 yN �1.18010 0.04477 �26.36 <2e�16 ⇤⇤⇤
1 # cons t ruc t a g r id o f p o s i t i o n s to p l o t the model aga in s t . . .
2 x new <� seq ( �1 ,1 , . 1 ) ; y new <� seq (�1 ,1 , .1)
3
4 # compute lambda f o r each po int on t h i s g r i d us ing the GLM model
5 lambda <� f unc t i on (x , y ) { exp ( p r ed i c t (model1 , data . frame (xN =x , yN=y ) ) ) }
6 z<� outer ( x new , y new , lambda )
7
8 # plo t lambda as a func t i on po s i t i o n over t h i s g r i d
9 persp (x new , y new , z , theta = 30 , phi = 30 , expand = 0 . 5 , c o l = ” l i g h t b l u e ” ) �> r e s
10 phi <� seq (0 , 2⇤pi , l en = 201)
11 xr <� cos ( phi )
12 yr <� s i n ( phi )
13 l i n e s ( trans3d ( xr , yr , lambda ( xr , yr ) , r e s ) , c o l = ”pink” , lwd = 2)
3
Question 3
One way to improve this exponential linear model is to add quadratic functions of xN and yN .
That is
�(t) = exp{�
0
+ �
1
xN (t) + �2yN (t) + �3xN (t)
2 + �
4
yN (t)
2 + �
5
xN : yN} (2)
all of the parameters except the one corresponding to the xN : yN term are significant. This
suggests that a more parsimonious model could be constructed by removing this term from the
model.
The AIC-values for exponential linear and quadratic model are 11394.2 and 9153.2 respectively.
Actually here we are fitting the Gaussian place field model, and it is able to more accurately
describe the firing properties of each of these neurons. Each place field has a center where firing
rate is maximal, and falls o↵ evenly when in any direction away. Even with the extraneous xN : yN
term, the Gaussian model has significantly lower AIC value.
1 # f i t a GLM model to the x and y po s i t i o n s . Adjust t h i s l i n e to improve model f i t .
2 model2 <� glm ( sp i k e s ˜xN⇤yN + I (xNˆ2) + I (yNˆ2) , po isson , data )
3 summary(model2 )
4 # v i s u a l i z e your model
5 # cons t ruc t a g r id o f p o s i t i o n s to p l o t the model aga in s t . . .
6 x new <� seq (�1 ,1 , .1)
7 y new <� seq (�1 ,1 , .1)
8
9 # compute lambda f o r each po int on t h i s g r i d us ing the GLM model
4
10 lambda <� f unc t i on (x , y ) { exp ( p r ed i c t (model2 , data . frame (xN=x ,yN=y ) ) ) }
11 z<� outer ( x new , y new , lambda )
12 # plo t lambda as a func t i on po s i t i o n over t h i s g r i d
13 persp (x new , y new , z , theta = 30 , phi = 30 , expand = 0 . 5 , c o l = ” l i g h t b l u e ” ) �> r e s
14 phi <� seq (0 , 2⇤pi , l en = 201)
15 xr <� cos ( phi )
16 yr <� s i n ( phi )
17 l i n e s ( trans3d ( xr , yr , lambda ( xr , yr ) , r e s ) , c o l = ”pink” , lwd = 2)
18 # Compute the AIC
19 model2$ dev iance +2⇤6 # 9153.217
20 model1$ dev iance +2⇤3 # 11394.16
Question 4
The AIC is minimized by a history component that goes back 8 ms. The optimal model from this
model class would include the constant term, the 4 quadratic terms with significant parameters
and 8 history terms.
1 l i b r a r y ( quantmod )
2 sp ik e h i s t = Lag ( sp ikes , k=1:20)
3 n=length (T)
4 #sp ike h i s t [ 1 : 2 0 , 1 : 2 0 ] = data [ 1 : 2 0 , 5 : 2 4 ]
5 idx=21:n
6 f o r (p in 1 : 20 ) {
7 model <� glm ( sp i k e s [ idx ] ˜xN [ idx ] ⇤yN[ idx ] + I (xN [ idx ] ˆ 2 ) + I (yN [ idx ] ˆ 2 ) + sp ike h i s t [ idx , 1 : p ] ,
8 f ami ly = poisson , data=data )
9 model a i c [ p ] = model$ dev iance + 2⇤ l ength (model$ c o e f f i c i e n t s )
10 }
11 p lo t ( 1 : 2 0 , model a ic , type=’ l ’ , c o l=’ blue ’ )
12 ab l i n e (h=min (model a i c ) , l t y =2, c o l=’ red ’ )
13 which . min (model a i c ) # 8
Question 5
Augment the GLM model from question 4 to include the past spiking history of the other recorded
neuron. And we find that none of the parameter estimates related to the history of the second
neuron is significant. This suggests that there is no functional interaction from neuron 2 to neuron
5
1, once position e↵ects have been modeled. The optiomal model remains the one with a constant
term, the 4 quadratic terms with significant parameters, 8 history terms and no interaction terms.
1 sp ik e h i s t = Lag ( sp ikes , k=1:20)
2 n=length (T)
3 sp ik e h i s t 2 = Lag ( sp ikes2 , k = 1 : 20 )
4 model <� glm ( sp i k e s [ idx ] ˜xN [ idx ] ⇤yN[ idx ] + I (xN [ idx ] ˆ 2 ) + I (yN [ idx ] ˆ 2 ) +
5 sp ik e h i s t [ idx , 1 : 8 ] + sp ike h i s t 2 [ idx , ] , f ami ly = poisson , data=data )
6 summary(model )
Question 6
When we eliminate the position related terms from the model, all of the interaction parameters
become significant. This is not surprising, since any spiking neuron 2 is likely to correspond to
periods when the animal is near the center of place field 2, and therefore, near the center of the
palce field of neuron 1. This model uses information about position that is important in the firing
of neuron 2 to predict the spiking of neuron 1. When we add the position term, this information is
already present and no significant interaction terms are observed.
1 idx=21:n
2 model <� glm ( sp i k e s [ idx ] ˜ I (xN [ idx ] ˆ 2 ) + I (yN [ idx ] ˆ 2 ) +
3 sp ik e h i s t [ idx , 1 : 8 ] + sp ike h i s t 2 [ idx , ] , f ami ly = poisson , data=data )
4 summary(model )
Question 7
The most parsimonious model does not completely pass the KS test. This suggests that there still
may be other statistical structure in the spiking activity that cannot be captured within this model
class. However, the KS statistic is relatively small suggesting that this model is able to explain
most of the observed structure in the data.
1 model <� glm ( sp i k e s ˜xN+yN + I (xNˆ2) + I (yNˆ2) +
2 s p i k e s h i s t 1+sp i k e s h i s t 2+sp i k e s h i s t 3+sp i k e s h i s t 4+sp i k e s h i s t 5+
3 s p i k e s h i s t 6+sp i k e s h i s t 7+sp i k e s h i s t8 , f ami ly = poisson , data=data )
4 summary(model )
5 lambda hat = exp ( p r ed i c t (model ) ) ; summary( lambda hat )
6 lambda hat = model$ f i t t e d . va lue s ; summary( lambda hat )
7 lambda in t = cumsum( lambda hat )
8 Z = d i f f ( c (0 , lambda i n t [ sp ikes >0]))
9 U = pexp (Z)
10 m=length (U)
11 p lo t ( s o r t (U) , ( 1 :m�1)/m, type=’ l ’ , lwd=1, c o l=’ blue ’ , x lab=’ Empir ica l CDF’ ,
12 ylab=’Model CDF’ )
13 l i n e s ( ( 1 :m�1)/m, ( 1 :m�1)/m+1.36/ sq r t (m) , l t y =2)
14 l i n e s ( ( 1 :m�1)/m, ( 1 :m�1)/m�1.36/ sq r t (m) , l t y =2)
15 t i t l e ( ’KS�p lo t ’ )
6
c
7