Skip to content

Commit 00c3563

Browse files
camila-mauracamila-maura
committed
finished official first draft of intro text
1 parent e184b29 commit 00c3563

File tree

1 file changed

+115
-52
lines changed

1 file changed

+115
-52
lines changed

‎docs/higher-order/GLM_pynapple_nemos.ipynb‎

Lines changed: 115 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -8,52 +8,83 @@
88
"# Generalized Linear Models using Pynapple & NeMos\n",
99
"In this notebook, we will use NeMos and Pynapple packages (supported by the [Flatiron Institute](https://neurorse.flatironinstitute.org)), to model spiking neural data using [Generalized Linear Models (GLM)](https://en.wikipedia.org/wiki/Generalized_linear_model). We will explain what GLMs are and which are their components, then use Pynapple and NeMos python packages to preprocess real data from the Allen Institute and use a GLM model to predict spiking neural data as a function of stimuli.\n",
1010
"\n",
11-
"A GLM is a regression model which trains a filter to predict a value (output) as it relates to some other variable (or input). \n",
11+
"A GLM is a regression model which trains a filter to predict a value (output) as it relates to some other variable (or input). It is called \"generalized\" because it constitutes a generalization of [ordinary linear regression](https://en.wikipedia.org/wiki/Linear_regression). Ordinary linear regression assumes that a constant change in a predictor leads to a constant change in the response variable, but this assumption is no longer useful for some types of response variables. In particular, when interested in modeling spikes, these are never expected to be negative. Moreover, spikes do not vary in a constant manner (usually, the variance of spike counts changes with the mean firing rate: neurons with average higher firing rates tend to have a higher variability than neurons with average lower firing rates). Thus, if interested in modeling this kind of phenomenon, the model of choice must account for these restrictions (:\n",
1212
"\n",
13-
":::{admonition} Why are GLMs called GLMs?\n",
14-
":class: info\n",
15-
":class: dropdown\n",
16-
"\n",
17-
"It is called \"generalized\" because it constitutes a generalization of [ordinary linear regression](https://en.wikipedia.org/wiki/Linear_regression). Ordinary linear regression assumes that a constant change in a predictor leads to a constant change in the response variable, but this assumption is no longer useful for some types of response variables. In particular, when interested in modeling spikes, these are never expected to be negative. Moreover, spikes do not vary in a constant manner (usually, the variance of spike counts changes with the mean firing rate: neurons with average higher firing rates tend to have a higher variability than neurons with average lower firing rates). Thus, if interested in modeling this kind of phenomena, the model of choice must account for these restrictions (:\n",
18-
":::\n",
19-
"\n",
20-
"\n",
21-
"\n",
22-
"In the neuroscience context, we can use a particular type of GLM to predict spikes: linear-nonlinear-Poisson (LNP) model.\n",
13+
"In the neuroscience context, we can use a particular type of GLM to predict spikes: the linear-nonlinear-Poisson (LNP) model.\n",
2314
"<figure>\n",
24-
"<img src=\"lnp_model.svg\" style=\"width:130%\"/>\n",
15+
"<img src=\"lnp_model.svg\" style=\"width:50%\"/>\n",
2516
"<figcaption align = \"center\"> LNP model schematic. Modified from <a href=\"https://www.nature.com/articles/nature07140\">Pillow et al., 2008</a></figcaption>\n",
2617
"</figure>\n",
27-
" That is, the model receives one or more inputs and then:\n",
28-
"\n",
29-
"1. Sends them through a linear filter or transformation\n",
18+
" This type of model receives one or more inputs and then sends them through a linear \"filter\" or transformation, passes said transformation through a nonlinearity to get the firing rate and uses that firing rate as the mean of a Poisson process to generate spikes. We will go through each of these steps one by one:\n",
3019
"\n",
31-
" The input (s) (also known as \"predictor(s)\") are first passed through a linear transformation: \n",
20+
"1. Sends them through a linear \"filter\" or transformation\n",
21+
" \n",
22+
" The inputs (also known as \"predictors\") are first passed through a linear transformation:\n",
3223
" \n",
33-
" $WX + c$\n",
24+
" $$\n",
25+
" \\begin{aligned}\n",
26+
" L(X) = WX + c\n",
27+
" \\end{aligned}\n",
28+
" $$\n",
29+
"\n",
30+
" Where $L$ is the filter, $X$ is the input (in matrix form), $W$ is a matrix and $c$ is a vector.\n",
3431
"\n",
35-
" This scales (makes bigger or smaller) or shifts (up or down) the input. When there is zero input, this is equivalent to changing the baseline rate of the neuron, which is how the intercept should be interpreted. So far, this is the same treatment of an ordinary linear regression. \n",
32+
" $L$ scales (makes bigger or smaller) or shifts (up or down) the input. When there is zero input, this is equivalent to changing the baseline rate of the neuron, which is how the intercept should be interpreted. So far, this is the same treatment of an ordinary linear regression. \n",
3633
"\n",
37-
"2. Passes the transformation through a nonlinearity to obtain a firing rate.\n",
34+
"2. Passes the transformation through a nonlinearity to get the firing rate.\n",
3835
" \n",
39-
" The aim of a GLM is to predict spiking activity. In particular, to predict a neuron's firing rate, which must be non-negative. This, as mentioned, is what the non-linearity part of the model handles: by passing the linear transformation through an exponential function, it is assured that it will always be non-negative. \n",
36+
" The aim of a LNP model is to predict the firing rate of a neuron and use it to generate spikes, but if we were only to keep $L(X)$ as it is, we would quickly notice that we could obtain negative values for firing rates, which makes no sense! This is what the nonlinearity part of the model handles: by passing the linear transformation through an exponential function, it is assured that the resulting firing rate will always be non-negative. \n",
4037
"\n",
41-
" As such, the firing rate, according to GLMs, would be defined as:\n",
38+
" As such, the firing rate in a LNP model is defined:\n",
4239
"\n",
43-
" $\\lambda = exp(WX + c)$\n"
40+
" $$\n",
41+
" \\begin{aligned}\n",
42+
" \\lambda = exp(L(X))\n",
43+
" \\end{aligned}\n",
44+
" $$\n",
45+
"\n",
46+
" where $\\lambda$ is a vector containing the firing rates corresponding to each timepoint."
4447
]
4548
},
4649
{
4750
"cell_type": "markdown",
4851
"id": "2a3ef026",
4952
"metadata": {},
5053
"source": [
51-
":::{admonition} A note on non-linearity\n",
54+
":::{admonition} A note on nonlinearity\n",
55+
":class: info\n",
56+
":class: dropdown\n",
57+
"\n",
58+
"\n",
59+
"In NeMoS, the nonlinearity is kept fixed. We default to the exponential, but a small number of other choices, such as soft-plus, are allowed. The allowed choices guarantee both the non-negativity constraint described above, as well as convexity, i.e. a single optimal solution. In principle, one could choose a more complex nonlinearity, but convexity is not guaranteed in general.\n",
60+
":::"
61+
]
62+
},
63+
{
64+
"cell_type": "markdown",
65+
"id": "7d32d362",
66+
"metadata": {},
67+
"source": [
68+
":::{admonition} What is the difference between a \"link function\" and the \"nonlinearity\"?\n",
5269
":class: info\n",
5370
":class: dropdown\n",
5471
"\n",
72+
"The link function states the relationship between the linear predictor and the mean of the distribution function. If $g$ is a link function, $L(⋅)$ is the linear predictor and $\\lambda$ the mean of the distribution function:\n",
73+
"\n",
74+
"$$\n",
75+
"\\begin{aligned}\n",
76+
"g(\\lambda) = L(⋅)\n",
77+
"\\end{aligned}\n",
78+
"$$\n",
79+
"\n",
80+
"$$\n",
81+
"\\begin{aligned}\n",
82+
"\\lambda = g^{-1}(L(⋅))\n",
83+
"\\end{aligned}\n",
84+
"$$\n",
85+
"\n",
86+
"the \"nonlinearity\" is the name for the inverse of the link function $g^{-1}(⋅)$.\n",
5587
"\n",
56-
"In NeMoS, the non-linearity is kept fixed. We default to the exponential, but a small number of other choices, such as soft-plus, are allowed. The allowed choices guarantee both the non-negativity constraint described above, as well as convexity, i.e. a single optimal solution. In principle, one could choose a more complex non-linearity, but convexity is not guaranteed in general.\n",
5788
":::"
5889
]
5990
},
@@ -66,23 +97,69 @@
6697
"\n",
6798
" A Poisson process is a special type of point process, in which the events are statistically independent. With these type of GLMs, each spike train is a sample from a Poisson process with the mean equal to the firing rate, i.e., output of the linear-nonlinear parts of the model. \n",
6899
"\n",
69-
" Remember, spiking is a stochastic process. That means that a given firing rate can give rise to a variety of different spike trains. Given that this is a stochastic process that could produce an infinite number of possible spike trains, how do we compare our model against the single observed spike train we have? We use the log-likelihood. This quantifies how likely it is to observe the given spike train for the computed firing rate: if $y(t)$ is the spike counts and $\\lambda(t)$ the firing rate, the equation for the log-likelihood is\n",
100+
" Remember, spiking is a stochastic process. That means that a given firing rate can give rise to a variety of different spike trains. Given that this is a stochastic process that could produce an infinite number of possible spike trains, how do we compare our model against the single observed spike train we have? We use the log-likelihood. This quantifies how likely it is to observe the given spike train for the computed firing rate: if $y(t)$ is the spike counts and $\\lambda(t)$ the firing rate at time $t$, the equation for the log-likelihood is\n",
101+
"\n",
102+
" $$\n",
103+
" \\begin{aligned}\n",
104+
"\n",
105+
" \\sum_{t}logP(y(t)|\\lambda(t)) = \\sum_{t}y(t)log(\\lambda(t)) - \\lambda(t)\n",
70106
"\n",
71-
" $\\sum_{t}logP(y(t)|\\lambda(t) = \\sum_{t}y(t)log(\\lambda(t))) - \\lambda(t) - log(y(t)!)$"
107+
" \\end{aligned}\n",
108+
" $$\n",
109+
"\n",
110+
" This is the objective function of the GLM model: we are trying to find the firing rate that maximizes the likelihood of the observed spike train."
111+
]
112+
},
113+
{
114+
"cell_type": "markdown",
115+
"id": "02063283",
116+
"metadata": {},
117+
"source": [
118+
":::{admonition} A note on log-likelihood of a Poisson distribution\n",
119+
":class: info\n",
120+
":class: dropdown\n",
121+
"\n",
122+
"To be precise, the log-likelihood of a Poisson distribution is:\n",
123+
"\n",
124+
"$$\n",
125+
"\\begin{aligned}\n",
126+
"\n",
127+
"\\sum_{t}logP(y(t)|\\lambda(t)) = \\sum_{t}y(t)log(\\lambda(t)) - \\lambda(t) - \\sum_{t}log(x(t)!)\n",
128+
"\n",
129+
"\\end{aligned}\n",
130+
"$$\n",
131+
"\n",
132+
"However, one can see that the term $- \\sum_{t}log(x(t)!)$ does not depend on $\\lambda$, and thus is independent of the model. Because of that, it is normally ignored.\n",
133+
"\n",
134+
":::"
135+
]
136+
},
137+
{
138+
"cell_type": "markdown",
139+
"id": "b2d41fee",
140+
"metadata": {},
141+
"source": [
142+
":::{admonition} What is the \"link function\" in the case of a LNP model?\n",
143+
":class: info\n",
144+
":class: dropdown\n",
145+
"\n",
146+
"In the case of a LNP model, the distribution function is a Poisson process with a mean of $\\lambda$. The \"non-linearity\", as mentioned before, is an exponential. So the \"link function\" would be the inverse, the logarithm!\n",
147+
":::"
72148
]
73149
},
74150
{
75151
"cell_type": "markdown",
76152
"id": "55cf8020",
77153
"metadata": {},
78154
"source": [
79-
":::{admonition} More resources\n",
155+
":::{admonition} More resources on GLMs\n",
80156
":class: seealso\n",
81157
":class: dropdown\n",
82158
"\n",
83-
"If you would like to learn more in depth about GLMs, you can refer to:\n",
159+
"If you would like to learn more about GLMs, you can refer to:\n",
84160
"\n",
85161
"- [Nemos GLM tutorial](https://nemos.readthedocs.io/en/latest/background/plot_00_conceptual_intro.html): for a bit more detailed explanation of all the components of a GLM within the NEMOS framework, as well as some nice visualizations of all the steps of the input transformation!\n",
162+
"- [Neuromatch Academy GLM tutorial](https://compneuro.neuromatch.io/tutorials/W1D3_GeneralizedLinearModels/student/W1D3_Tutorial1.html): for a bit more detailed explanation of the components of a GLM, slides and some coding exercises to ensure comprehension.\n",
86163
"- [Youtube video on LNP](https://www.youtube.com/watch?v=i62gffPrZYA): Although outside of the context of Neuroscience, it does go step by step explaining what LNPs are with visualizations and notes on limitations of the model (:\n",
87164
":::"
88165
]
@@ -92,8 +169,7 @@
92169
"id": "c5945cfd",
93170
"metadata": {},
94171
"source": [
95-
"\n",
96-
"We will be analyzing data from the [Visual Coding - Neuropixels dataset](https://portal.brain-map.org/circuits-behavior/visual-coding-neuropixels), published by the Allen Institute. This dataset uses [extracellular electrophysiology probes](https://www.nature.com/articles/nature24636) to record spikes from multiple regions in the brain during passive visual stimulation. For simplicity, we will focus on the activity of neurons in the visual cortex (VISp) during passive visual stimulation: full-field flashes, of color either black or white. \n",
172+
"For this tutorial, we will be analyzing data from the [Visual Coding - Neuropixels dataset](https://portal.brain-map.org/circuits-behavior/visual-coding-neuropixels), published by the Allen Institute. This dataset uses [extracellular electrophysiology probes](https://www.nature.com/articles/nature24636) to record spikes from multiple regions in the brain during passive visual stimulation. For simplicity, we will focus on the activity of neurons in the visual cortex (VISp) during passive visual stimulation: full-field flashes, of color either black or white. \n",
97173
"\n",
98174
"Our aim is to model spiking activity from neurons (spiking rates and spike trains) from the Primary Visual cortex (VISp) as a function of the presented visual stimuli. Moreover, we will also compare this model including an extra predictor: spike history, to see whether adding history as a predictor improves the performance of the model."
99175
]
@@ -1959,27 +2035,6 @@
19592035
"### References"
19602036
]
19612037
},
1962-
{
1963-
"cell_type": "markdown",
1964-
"id": "228bc865",
1965-
"metadata": {},
1966-
"source": [
1967-
"### Delete later\n",
1968-
"\n",
1969-
"<div class=\"admonition info\">\n",
1970-
"<p class=\"admonition-title\">Resources</p>\n",
1971-
"<p>\n",
1972-
"If you would like to learn more in depth about GLMs, you can refer to:\n",
1973-
"\n",
1974-
"<p>\n",
1975-
"\n",
1976-
"- [Nemos GLM tutorial](https://nemos.readthedocs.io/en/latest/background/plot_00_conceptual_intro.html): for a bit more detailed explanation of all the components of a GLM within the nemos framework, as well as some nice visualizations of all the steps of the input transformation!</p>\n",
1977-
"\n",
1978-
"\n",
1979-
"</p>\n",
1980-
"</div>"
1981-
]
1982-
},
19832038
{
19842039
"cell_type": "markdown",
19852040
"id": "79dec3fc",
@@ -2019,8 +2074,16 @@
20192074
"- improve responsiveness explanation/notation\n",
20202075
"- Should I z-score firing rates so its easier to compare? in single unit raster\n",
20212076
"- Add smoothing to plotting added to hist!!!\n",
2022-
"- Add single cell prediction perievent"
2077+
"- Add single cell prediction perievent\n",
2078+
"- Read all of the texts\n",
2079+
"- center image"
20232080
]
2081+
},
2082+
{
2083+
"cell_type": "markdown",
2084+
"id": "338acb72",
2085+
"metadata": {},
2086+
"source": []
20242087
}
20252088
],
20262089
"metadata": {

0 commit comments

Comments
 (0)