程序代写代做代考 MA 568 Statistical Analysis of Point Process Data

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