Neurons fire electric signal in order to transmit sensory information throughout the body. Can the intensity of a stimulus predict how frequently a neuron fires?

• Predictor: Stimulus intensity
• Response: Firing rate in Hz
• Method: Simple linear regression

For this lab, I manually modeled equation of firing rate as simple linear regression, where Y is firing rate, X is stimulus intensity, and ϵ is random error.

\[Y_i = 5 + 3X_i + \epsilon_i, \quad \epsilon \sim N(0, 5^2)\]

set.seed(123)
stimulus_intensity <- runif(100, 0, 10)

firing_rate <- 5 + 3 * stimulus_intensity + rnorm(100, mean = 0, sd = 5)

df <- data.frame(
  stimulus_intensity, firing_rate)

summary(df)

##  stimulus_intensity  firing_rate   
##  Min.   :0.006248   Min.   :-2.18  
##  1st Qu.:2.454707   1st Qu.:12.95  
##  Median :4.663708   Median :19.49  
##  Mean   :4.985590   Mean   :19.69  
##  3rd Qu.:7.554712   3rd Qu.:27.26  
##  Max.   :9.942698   Max.   :43.89

• The parameters of a simple linear regression model are the intercept \(\beta_0\) and the slope \(\beta_1\). The intercept is the expected value of the response when the predictor (explanatory variable) is zero, and the slope represents the expected change in the response for a one-unit increase in the predictor. Because the true parameters are unknown, we estimate them from the data using \(\hat{\beta}_0\) and \(\hat{\beta}_1\).

model <- lm(firing_rate ~ stimulus_intensity,
            data = df)
summary(model)

## 
## Call:
## lm(formula = firing_rate ~ stimulus_intensity, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1899  -3.0661  -0.0987   2.9817  11.0861 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.9552     0.9803   5.055    2e-06 ***
## stimulus_intensity   2.9551     0.1709  17.290   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.846 on 98 degrees of freedom
## Multiple R-squared:  0.7531, Adjusted R-squared:  0.7506 
## F-statistic:   299 on 1 and 98 DF,  p-value: < 2.2e-16

• p-value: < 2.2e-16, smaller than 0.001

• parameters

coef(model)

##        (Intercept) stimulus_intensity 
##           4.955201           2.955084

• confident interval at 95% confidence interval

confint(model, level=0.95)

##                       2.5 %   97.5 %
## (Intercept)        3.009877 6.900525
## stimulus_intensity 2.615918 3.294251

\[\begin{aligned} \hat{\beta}_0 &= 4.9552, \quad \hat{\beta}_1 = 2.9551 \\ \hat{Y}_i &= 4.9552 + 2.9551X_i \end{aligned}\] The parameters of linear model are close to the original parameters of original model.

plt <- ggplot(df,
       aes(x=stimulus_intensity,
           y=firing_rate)) +
  geom_point() +
  geom_smooth(method='lm', se=TRUE) +
  labs(x ='Stimulus intensity',
       y = 'Firing rate')

plt

## `geom_smooth()` using formula = 'y ~ x'

The scatterplot is observed data, and the blue line is fitted regression line which estimates mean firing rate at each level of stimulus intensity. The shaded area represents 95% confidence interval around the estimated regression line.

• Residual Errors: TSS measures the total variation of the model, ESS measures explained variation of model, and RSS measures unexplained error of model. The best fitting line is where the RSS is minimized.

• $$TSS = ESS + RSS$$

  $$TSS = \sum_{i=1}^{n}(y_i-\bar{y})^2$$

  $$ESS = \sum_{i=1}^{n}(\hat{y}_i-\bar{y})^2$$

  $$RSS = \sum_{i=1}^{n}(y_i-\hat{y}_i)^2$$

df$fitted <- fitted(model)
# column fitted values, yhat
y_bar <- mean(df$firing_rate)

ggplotly(tss_plot)

ggplotly(ess_plot)

## `geom_smooth()` using formula = 'y ~ x'

ggplotly(rss_plot)

## `geom_smooth()` using formula = 'y ~ x'

This plot shows the residuals across different levels of stimulus intensity. Each point represents the difference between the observed and fitted firing rate. The red dashed line represents zero residual.

residual >0 means model is overpredicting and residual <0 means model is underpredicting.

Because the errors were generated randomly, the residuals are randomly scattered around zero without a clear pattern, which suggests that the linear model fits the data.

The residuals are approximately bell-shaped and centered around zero. Since we assumed that the errors are normally distributed, the histogram is consistent with the initial equation.