GENERALIZED LINEAR
MODELHideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM Linear modelLinear modelParameters

Measurement noise

Y = Xβ + ε

ε ~ N ( 0, σ 2 )Linear predictor
Dependent variable

Explanatory variables

Probabilistic representation: A normal distributionf ( y | µ, σ 2 ) =12πσ2e( y−µ )−σ22,where the mean is modeled as µ = X β.
Linear predictor
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Generalized linear modelWhat if the measurement noise follows non-Gaussian distribution?
GLM-analysis extends the regression analysis when the residual isgiven by an exponential family of distributions.

The exponential family of distributions include many popular distributions:normal, log-normal, exponential, gamma, chi-squared, beta, Dirichlet, Bernoulli,
categorical, Poisson, inverse Gaussian, Wishart, Inverse Wishart, von Mises, etc.

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Exponential family of distributions
Parameter(s) of the family(Sufficient statistics)

f ( y | θ ) = h ( y) exp "#η (θ ) T ( y) − A (θ )$%Natural parameter

log-partition function

By re-parameterizing the distribution by η , an exponential family issimplified to a ‘natural form’ (canonical form) asf ( y | η ) = h ( y) exp "#ηT ( y) − A (η )$%Relation between a natural parameter and log-partition function:∂A (η )= EY "#T (Y )$%∂ηHideaki Shimazaki, Ph.D. http://goo.gl/viSNG

∂2 A (η )2#%= EY $T (Y ) − EY T (Y )& = var #$T (Y )%&2∂η
Point process-GLM Examples
Bernoulli probability mass function.
y1−yf ( y; µ ) = µ (1− µ )"%µ= exp $ y log+ log (1− µ )'1− µ#&y = {0,1}

Poisson probability mass function.µ y −µ 1f ( y; µ ) = e = exp [ y log µ − µ ]y!
y!

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

y = {0,1, 2,}Point process-GLM 

Generalized linear model
Modeling the mean by a linear predictor.

f ( y;θ ) = h ( y) exp "#η (θ ) T ( y) − A (θ )$%Meanµ ≡ EY [Y ]Mean is boundedLinear predictor X β0 < µ < 1 (Bernoulli)
0 < µ (Poisson)

A link function relates the mean anda linear predictor:µ = Xβ1.0µ=g−1(Xβ )0.80.6µ = g−1 ( X β )X β = g (µ )g−1 ( x ) =11+ e− x0.40.2so that µ is bounded.
-4-224Xβ

A link function can be selected arbitrary. However, there is a functionthat arises ‘naturally’, depending on the type of distribution.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM

Canonical link function: Bernoulli
Poisson probability mass function.
y1−yf ( y; µ ) = µ (1− µ )"%µ= exp $ y log+ log (1− µ )'1− µ#&y = {0,1}Linear predictor:  η = X βModel the mean using a link function:  µ = g
Canonical link function:(logit function)

1−1µ=gη=,( )If we use −η1+ eµη = g (µ ) = log,namely 
1− µthen the distribution is written in a naturalform as

−1(η )
sigmoid function

µ1.00.80.6µ = g−1 (η )0.4

0.2-4-22η

4logit functionη420.2

0.4

0.6

0.8

1.0-2f ( y;η ) = exp #$ yη − A" (η )%&Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

-4-6

η = g (µ )µPoint process-GLM

Canonical link function: Poisson
Bernoulli probability mass function.µ y −µ 1f ( y; µ ) = e = exp [ y log µ − µ ]y!
y!

y = {0,1, 2,}Linear predictor:  η = X βModel the mean using a link function:  µ = g
Canonical link function:(log function)
If we use µ=g−1−1(η )Exp. function

µ6050

η

(η ) = e ,40µ = g−1 (η )

30
20
10

namely η = g (µ ) = log µ,then the distribution is written in a naturalform as

1f ( y;η ) = exp "# yη − eη $%y!

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

-4-22η4log function

η0.2

0.4

0.6

0.8

1.0-1-2-3-4

η = g (µ )µPoint process-GLM

Relating a stimulus input/motor output with neural activity

POINT PROCESS-GLM

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM

Relating stimulus and neural activity
Stimulus signal x (t )(Conditional intensity λ t | x (t )Point process

Nt)

dN t

t t+Δ
P ( a spike in [t, t + Δ) | x (t )) = λ (t | x (t )) Δ + o ( Δ )Stimulus signal (Covariate)
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Discretizing time into bins
Stimulus signal x (t )kth-binxk (≡ x ( kΔ ))(Conditional intensity λ t | x (t ))λk (≡ λ ( kΔ ))Point process N t

1

100110001010100100

ΔN kΔ
Likelihood (Joint probability mass function)
N(P ( ΔN Δ , ΔN 2Δ ,, ΔN NΔ ) = ∏ P ΔN kΔ | ΔN Δ ,, ΔN (k−1)Δ
k=1

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG)

Conditional probability mass function
Point process-GLM Bernoulli-GLM
Conditional Bernoulli probability:()

P ΔN kΔ | ΔN Δ ,, ΔN (k−1)Δ = ( λk Δ )ΔN kΔ

1−ΔN kΔλk ≡ λ ( kΔ | x ( kΔ ))(1− λk Δ)#&λk Δ= exp %ΔN kΔ log+ log {1− λk Δ}(1−λΔ$'kMean of the Bernoulli probability at k-th bin:µ = E ( ΔN kΔ ) = λk Δ

Modeling the mean response via a canonical link function:

λk Δ =

11+ e ( 0− β + β xk )log

sigmoid function

logit function

1.00.8λk Δ= β 0 + β xk1− λk Δ

420.6

0.2

0.4

0.4

0.6

0.8

1.0-2

0.2-4-4-2

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

24-6misprint!
Point process-GLM Poisson-GLM
Approximation by Poisson distribution:
ΔN# λk Δ & kΔ
P ΔN kΔ | ΔN Δ ,, ΔN (k−1)Δ = %( (1− λk Δ )$ 1− λk Δ '≈ exp *+ΔN kΔ log {λk Δ} − λk Δ,-()

Mean of the Poisson distribution at k-th bin:µ = E ( ΔN kΔ ) = λk Δ

Modeling the mean response via a canonical link function:

λk Δ = exp [ β0 + β xk ]log [ λk Δ ] = β0 + β xk

Exp. function

log function
0.260
50

0.4

0.6

0.8

1.0-140-230
20-310-4-2

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

24-4

Point process-GLM 

Relation between Bernoulli and Poisson-GLM
Inverse link function

Logit (Bernoulli)

2.5Exp (Poisson)

λΔ = exp [µ ]Link function
4log

2.02λΔ=µ1− λΔ

1.5
0.2

0.4

0.60.81.0

1.0-20.5-10-8-6-4-2Sigmoid (Bernoulli)
2λΔ =11+ e− µ4-4Log (Poisson)

log [ λΔ ] = µmisprint

If the mean λΔ is small, an exponential and sigmoid functions are nearly
identical. Similarly, log and sigmoid functions are nearly identical for small λΔ. Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Is the point process a Poisson?
We now recall that the number of spike count in a certain interval follows
the Poisson distribution. The spike count in a small bin is then approximated as()P NΔ = n =(λΔ)

nnλΔ ) #&2(1− λΔe =%1− λΔ + ( λΔ ) +(n!
n! $2'In particular,!
$
21P N Δ = 0 = 1#1− λΔ + λΔ +& = 1− λΔ + o Δ
2
"%!
$
21
P N Δ = 1 = λΔ #1− λΔ + λΔ +& = λΔ + o Δ
2"%$2!
21P N Δ = 2 = λΔ #1− λΔ + λΔ +& = o Δ
2"%

()

( )

( )

()

( )

( )

() ( )( )( )The probability of having a spike/no spike is an approximation of thePoisson count distribution for a small time bin.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM From discrete to continuous timeLikelihood of discrete-time conditional Poisson distributions: KP ( ΔN Δ , ΔN 2Δ ,, ΔN NΔ ) = ∏ exp #$ΔN kΔ log {λk Δ} − λk Δ%& + o ( Δ n )k=1#K%= Δ exp )∑ ΔN kΔ log λk − λk Δ* + o ( Δ n )$ k=1&nLikelihood of a point process by continuous-time limit:

P ( ΔN Δ , ΔN 2Δ ,, ΔN NΔ )p (t1, t2 ,…, tn ∩ N T = n ) = lim
Δ→0ΔnK&K) nexp (∑ ΔN kΔ log {λk } − ∑ λk Δ+ Δ + o ( Δ n )' k=1*k=1= lim
Δ→0ΔnT
T= exp &( ∫ log λu dN u − ∫ λu du)+0' 0*Truccolo et al. J Neurophysiol (2005)
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM

Point process
Covariate signal Xtxt

Point processNt

dN tt t+ΔP ( a spike in [t, t + Δ) | H t , Xt ) = λ (t | H t , Xt ) Δ + o ( Δ )History of spiking activity

H t = {t1 < t2 < < tn ∩ N t = n}
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

History of a covariate signal(Filtration)

X t = { xu : 0 < u < t }Point process-GLM

Point process-GLM
Likelihood function:nTp (t1, t2 ,…, tn ∩ N T = n ) = ∏ λ (ti | Xti , H ti ) exp %'− ∫ λ (u | Xu , H u ) du(*& 0)i=1nHere dN t =TT%
= exp ' ∫ log λ (u | Xu , H u ) dN u − ∫ λ (u | Xu , H u ) du(*0& 0)∑δ (t − t )dti

i=1Modeling the conditional mean rate via a canonical link function.

! dN t $E#= λ (t | H t , Xt ) = exp !"β0 + β1 f ( Xt ) + β2 g ( H t )$%&" dt %a log-link function
A canonical link function for a Poisson distribution.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLMCovariates for neural signalsThe covariates can be added in a (generalized) linear manner.

log λ (t | H t , Xt ) = β0 + β1 [ Input stimulus]
Log-link function+β1! [ Higher-order feature of stimulus]+β2 [Spike history]+β3 [ Ensemble spike history]+β 4 [ Local field potential]

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM

Maximum likelihood estimationLikelihood TTL ( β | N 0:T ) = exp #% ∫ log λ (u | Xti , H ti ) dN u − ∫ λ (u | Xu , H u ) du&(0$ 0'whereλ (t | H t , Xt ) = exp !" x (t ) β #$.

Optimize parameters under the principle of maximizing the likelihood:

β MLE = arg max log L ( β | N 0:T )βNewton-Raphson method:l ( β ) = log L ( β | N 0:T )β j+1 = β j − ∇∇l ( β j ) ⋅ ∇l ( β j )Hessian Score function

Matlab:
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

b = glmfit ( X, y, ' poisson')
Point process-GLM Applications of point process-GLM

2005

Here we reviewHideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM Modeling conditional intensity function

2005

Conditional intensity function:

Auto spike history Ensemble history Extrinsic covariate
Auto spike history

Expected to capture refractory effects, recovery periods, and oscillation.
Determined by AIC

Ensemble history

Expected to capture lagged spike synchrony between target and other cells.

Expected to capture ensemble effect on slower time scale

Extrinsic covariate (hand velocity model)

Spike count in a window of size W.
Angle of hand velocity
Magnitude of hand velocity

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

offset τ=+150ms
Point process-GLM 2005

Simulation (6 neurons)
Simulated conditional intensity ofa blue neuron
ISI

Coefficients ofauto-spikehistory effect

Excitatory

KS plot

History
Inhibitory

Excitatory
Inhibitory
Other neurons

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Application to MI spiking dataTwo models were compared. 
Velocity model
Velocity + Auto spike history model
Tuning curve of velocity modelHideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Auto spike history

Point process-GLM Test of goodness-of-fit: KS-plot
KS plot

is expected to be independent
uniformly distributed randomvariables in an unit interval.

Velocity model

Velocity + Auto-spikehistory model
95% confidencebound of KS test

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM Test of independence for rescaled ISIs

2005

If the conditional intensity model completely explains the data, is expected to be independently and uniformly distributed random variables
in an unit interval.

There is dependency
between zj and zj+1.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

There is no dependency
between zj and zj+1.
Point process-GLM

Model selection by AIC
Akaike information criterion(( ))AIC = −2log p DATA | β + 2 × dim βParameters of the model

Model dimension(Number of parameters)

The smaller the AIC, the better the predictive ability of the model is.The best model

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

AR: Auto-spike history
The models with the auto-spikehistory significantly decrease AIC. 

Point process-GLM

Point process residual analysis

2005

Point process residual

The difference between the actual spike count and spike count predicted bythe model.If the CIF completely characterizes response characteristics to covariates(stimulus/motor output), the residual should have no dependency with thecovariates.
Cross-correlation function between residual and covariates (hand movement).

direction

speed

Inclusion of velocity signal in a CIF
greatly reduced the correlation.
However, the auto-spike history did
not further reduce the correlation.
AR component does not carry additional
statistical information about hand velocity.

There is a component that was not
captured by the point process-GLMs.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLMWe will review
multi-neuron
GLM
Point process-GLM for the analysisof retinal
ganglion cells
exponential probabilistic
spiking
nonlinearity

stimulus filter+
post-spike filter
neuron 1

stimulus

coupling
filters+
neuron 2

Pillow et al (2008). Nature, 454(7207), 995–999. Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM

What we learned

1• Exponential family of distributions and generalizedliner model.

2• Discrete-time likelihood of conditional Bernoulli andPoisson distributions.

34• GLM framework for a continuous-time point process.
• Review of a paper by Truccolo et al. Focused on methods
for model validation (KS-plot, AIC, residual analysis).

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM 

Course overview
1• Introduction and a Poisson point process2• Renewal and non-Poisson processes

3• Point process-GLM: Stimulus and a point process (Encoding)4• Inference for a Poisson process: Spike-rate estimation

5
6• State-space model and a point process filter (Decoding)• Paper reviews by SCS 1st year students

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM From the feedback of the last lecture
Stimulus signal

x (t )

Instantaneous spike-rate λ (t )Point process
Estimation

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG

Point process-GLM