---
title: "Bold or Bashful? How Fur Color Predicts Human Approach Behavior in NYC's Central Park Squirrels"
subtitle: "SOCI 252: Final Project"
author: "Akshara Kolipaka"
date: today
format:
html:
theme: cosmo
toc: true
toc-depth: 2
toc-title: "Contents"
number-sections: true
code-fold: true
code-summary: "Show the code"
code-tools: true
embed-resMyces: true
execute:
echo: true
warning: false
message: false
---
```{r setup, echo=FALSE, warning=FALSE}
#| include: false
knitr::opts_chunk$set(
collapse = TRUE,
warning = FALSE,
message = FALSE,
error = TRUE,
fig.height = 4,
fig.width = 6,
fig.env = 'figure',
fig.pos = 'h',
fig.align = 'center')
library(patchwork)
library(tidyverse)
library(broom)
library(table1)
nyc_squirrels <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2019/2019-10-29/nyc_squirrels.csv")
squirrels <- nyc_squirrels %>%
filter(!is.na(primary_fur_color)) %>%
mutate(
boldness = as.integer(approaches) - as.integer(runs_from),
primary_fur_color = factor(primary_fur_color, levels = c("Gray", "Cinnamon", "Black")),
age = factor(age),
shift = factor(shift)
)
squirrels <- squirrels %>%
filter(age %in% c("Adult", "Juvenile")) %>%
mutate(age = factor(age, levels = c("Adult", "Juvenile")))
hectare_density <- squirrels %>%
count(hectare, name = "squirrel_density")
squirrels <- squirrels %>%
left_join(hectare_density, by = "hectare")
```
## Introduction
Central Park is home to thousands of cute and amazing squirrels with different personalities. Some squirrels will eat out of your hand while others scatter the moment you get close. In this project I investigate if difference in squirrel personality has anything to do with fur color. Specifically, I want to know if Gray, Cinnamon, and Black squirrels differ in how they respond to humans.
This question matters because urban squirrels live in unusually close contact with people in Central Park. Researching about which squirrels are bolder or more avoidant can tell us something meaningful about how animals adapt to city life. In wildlife biology, color differences within a species are sometimes linked to differences in personality and stress response. If that is true for Central Park squirrels, fur color could serve as a simple, observable indicator of behavioral tendencies.
Most people would probably assume all Central Park squirrels behave the same way, since they all share the same park and have grown up around humans. But, what I think is less obvious is whether the rarer color squirrels like Cinnamon and Black, have distinct behavioral profiles compared to the much more common Gray squirrel.
The data for this project comes from the 2018 NYC Squirrel Census, a citizen science project in which trained volunteers observed and recorded squirrels across Central Park over two weeks in October. The dataset is publicly available through the NYC Open Data portal and was shared through the TidyTuesday project.
## Data Collection
The 2018 NYC Squirrel Census was conducted in October 2018 in Central Park, Manhattan where volunteer sighters organized by the Squirrel Census team observed squirrels across the park. Each volunteer was assigned a section of the park and recorded detailed information about every squirrel they spotted during two-week morning and afternoon sessions. The goal was to create a repeatable count of the Central Park squirrel population.
To clean my dataset I loaded the raw NYC Squirrel Census data directly from TidyTuesday and removed any rows missing a fur color, since fur color is my main explanatory variable of interest. Next, I created a new outcome variable called `boldness`. The raw dataset records two separate TRUE/FALSE columns:
- `approaches` marks whether a squirrel moved toward the human observer
- `runs_from` marks whether it fled.
Since neither column alone captures the full picture, I combined these two columns into a single numeric score by converting both to integers (TRUE = 1, FALSE = 0) and subtracting `runs_from` from `approaches`. This gave me a three-point scale: +1 if the squirrel approached, -1 if it ran away, and 0 if it ignored the observer.
I also converted three variables into factors so R treats them as categorical in the models. These factors are:
- `primary_fur_color` (the dominant coat color: Gray, Cinnamon, or Black). Gray is set as the reference level for fur color since it is the most common color type in Central Park.
- `age` (Adult or Juvenile)
- `shift` (AM or PM sighting session).
Finally, I created `squirrel_density`. The park was divided into a grid of hectares for the census, and this variable counts how many squirrels were observed in each hectare. That count is then joined back onto the main dataset so every sighting has a measure of how crowded its local patch of the park was. I included this variable as a control in the full model.
There are 2,968 cases in the cleaned dataset. Each row represents one unique squirrel sighting recorded by a volunteer observer. The final squirrels dataset contains 38 variables in total after I added `boldness` and `squirrel_density`. The first six rows of the cleaned dataset are shown below.
```{r}
head(squirrels)
dim(squirrels)
```
My primary explanatory variable is `primary_fur_color`. I chose it because fur color is the most visually fun physical trait recorded in this dataset, and wildlife biology research suggests fur color within a species can reflect differences in stress response and behavioral boldness. An example of this research can be found here: https://pmc.ncbi.nlm.nih.gov/articles/PMC11617328/
My secondary quantitative explanatory variable is `squirrel_density`, which I created from the data. I tought that squirrels in more crowded hectares may face more competition for food, which could push them toward bolder behavior around humans regardless of their fur color. Therefore, including this variable helps me isolate the effect of fur color on boldness.
The key variables used in this analysis are:
- `boldness`: A composite score I created. It equals +1 if the squirrel approached the observer, -1 if it ran away, and 0 if it was indifferent. This is the outcome variable.
- `primary_fur_color`: The dominant coat color of the squirrel. Values are Gray, Cinnamon, or Black.
- `squirrel_density`: The total number of squirrels observed in the same hectare across the census. Higher values mean a more crowded local environment.
- `age`: Whether the squirrel appeared to be an Adult or Juvenile at the time of sighting.
- `shift`: Whether the sighting took place during the AM or PM observation session.
- `foraging`: A TRUE/FALSE indicator for whether the squirrel was observed searching for food at the time of the sighting.
- `hectare`: The section of Central Park where the squirrel was spotted. The park was divided into a grid of hectares and each hectare can be identified by a number and letter code.
## Data Exploration
### Response Variable
My response variable is `boldness`. This is a composite score I created to measure how each squirrel behaved toward the human observer. A score of +1 means the squirrel approached the human observer., 0 means it was indifferent, and -1 means it fled from the human observer.. I use this variable instead of the raw `approaches` and `runs_from` variable columns because it captures the full behavioral picture in a single numeric value, which really helps me for linear regression.
To understand how boldness is distributed across the sample, I created a bar chart and summary statistics below.
```{r}
p_bold <- squirrels %>%
ggplot(aes(x = boldness)) +
geom_bar(color = "black", fill = "#F6B8D0", alpha = 0.6) +
scale_x_continuous(breaks = c(-1, 0, 1),
labels = c("-1 (Fled)", "0 (Indifferent)", "1 (Approached)")) +
labs(
title = "Distribution of Squirrel Boldness Scores",
x = "Boldness Score",
y = "Count"
) +
theme_minimal()
p_bold
squirrels %>%
summarise(
mean_boldness = mean(boldness, na.rm = TRUE),
sd_boldness = sd(boldness, na.rm = TRUE),
min_boldness = min(boldness, na.rm = TRUE),
max_boldness = max(boldness, na.rm = TRUE),
n = n()
)
```
The distribution of boldness is heavily skewed toward 0. Out of 2,968 squirrels, the vast majority were indifferent to the observer, around 650 squirrels fled, and fewer than 200 squirrels approached the human observers. The mean boldness score is -0.17. This means that the average squirrel leans to be slightly avoidant rather than bold. The standard deviation is 0.49, which reflects the fact that most values cluster tightly at 0 with relatively few squirrels at the extremes. Boldness ranges from -1 to +1 as expected. This distribution tells us that approaching behavior is rare in Central Park squirrels, and that most variation in the data comes from the difference between fleeing and being indifferent rather than between fleeing and approaching.
### Key Explanatory Variables
#### Squirrel Density Per Hectare
I first looked at how local squirrel density relates to boldness. The logic behind this is that in more crowded hectares, squirrels face more competition for natural food, which may push them toward bolder behavior around humans as an alternative food supply.
```{r}
p_density_scatter <- squirrels %>%
ggplot(aes(x = squirrel_density, y = boldness)) +
geom_jitter(alpha = 0.2, height = 0.1, color = "#5C8374") +
geom_smooth(method = "lm", se = TRUE, color = "#183D3D") +
labs(
title = "Squirrel Density vs. Boldness",
x = "Squirrel Density (per hectare)",
y = "Boldness Score"
) +
theme_minimal()
p_density_scatter
squirrels %>%
summarise(
mean_density = mean(squirrel_density, na.rm = TRUE),
sd_density = sd(squirrel_density, na.rm = TRUE),
min_density = min(squirrel_density, na.rm = TRUE),
max_density = max(squirrel_density, na.rm = TRUE)
)
```
The boldness variable can only take three discrete values: -1, 0, and 1. Thefore if regular geom_point scatterplot was used, thousands of dots would just stack directly on top of each other at those three horizontal bands. This makes it impossible to see how many squirrels fall at each point.Thefore I used geom_jitter to add a small amount of random vertical noise to each point so the density of observations at each value becomes visible. The horizontal spread of points across density values is real data. The vertical spread within each band is just the jitter.
The Hectare density ranges from 1 to 32 squirrels, with a mean of 12.51 and standard deviation of 6.49. This means that most hectares had somewhere between 6 and 19 squirrels observed. The scatter plot shows that the regression line is nearly flat. This suggestes that there is a very weak relationship between local crowding and boldness. The points are spread across all three boldness values at every density level, and the confidence band is narrow, which indicates that this flat trend is estimated with reasonable certainty. Overall, squirrel density does not appear to be a strong predictor of boldness on its own.
#### Primary Fur Color
My primary categorical explanatory variable is `primary_fur_color`. Since boldness only takes three values, I used a stacked proportional bar graph. This shows what share of each fur color group fled, was indifferent, or approached the observer.
```{r}
squirrels %>%
mutate(boldness_label = case_when(
boldness == -1 ~ "Fled",
boldness == 0 ~ "Indifferent",
boldness == 1 ~ "Approached"
),
boldness_label = factor(boldness_label, levels = c("Fled", "Indifferent", "Approached"))) %>%
group_by(primary_fur_color, boldness_label) %>%
summarise(n = n()) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = primary_fur_color, y = prop, fill = boldness_label)) +
geom_bar(stat = "identity") +
scale_fill_manual(values = c("Fled" = "#FF7F7F", "Indifferent" = "lightyellow", "Approached" = "lightgreen")) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Boldness Behavior by Primary Fur Color",
x = "Primary Fur Color",
y = "Proportion of Squirrels",
fill = "Behavior"
) +
theme_minimal()
squirrels %>%
group_by(primary_fur_color) %>%
summarise(
n = n(),
mean_bold = mean(boldness, na.rm = TRUE),
sd_bold = sd(boldness, na.rm = TRUE)
)
```
#### Full Summary Table
```{r}
summary_table <- table1(~ boldness + squirrel_density + age + shift + foraging | primary_fur_color,data = squirrels)
summary_table
```
Here I created a summary table breaks down all key variables by fur color group.
For boldness, Cinnamon squirrels have the highest mean at -0.110, followed by Gray squirrels at -0.172, and Black squirrels at -0.252. All three squirrelsgroups have a median of 0. This confirms that indifference is the most common behavior regardless of fur color.
Squirrel density is very similar across all three groups. Gray squirrels have a mean of 12.5, Cinnamon squirrels have a mean of 12.6, and Black squirrels have a mean of 11.9 This is helpful because it suggests that the boldness differences between fur color groups are not simply a result of one group living in more crowded parts of the park.
The age breakdown is also similar across groups. Adult squirrels make up roughly 83 to 89 percent of each group. There are 101 sightings with missing age values and 2 recorded as "?". This will be excluded by R when I fit the models.
The shift variable is nearly evenly split across all three groups, with slightly more PM sightings than AM in every group. Foraging rates are also comparable, ranging from 41.7% for Black squirrels to 51.8% for Cinnamon squirrels.
Overall, the three fur color groups look similar across all control variables. This strengthens the case that any boldness differences found in the models are more likely related to fur color itself rather than to differences in age, time of day, or local density.
## Modeling
### One Quantitative Explanatory Variable: Squirrel Density
Here I want to find out if how crowded a hectare is (measured by squirrel_density) tells me anything about how bold a squirrel acts toward humans. The idea is that squirrels living in more crowded areas might have to compete harder for natural food, so they might be more willing to approach a human for food as a backup option. To test this, I used a simple linear regression model.
My expected model equation is:
$$E[\text{boldness} \mid \text{squirrel\_density}] = \beta_0 + \beta_1 \cdot \text{squirrel\_density}$$
Here, Y is boldness and X is squirrel density. The equation says that my model predicts the average boldness score using an intercept plus a slope times density.
```{r}
model_density <- lm(boldness ~ squirrel_density, data = squirrels)
squirrels %>%
ggplot(aes(x = squirrel_density, y = boldness)) +
geom_jitter(alpha = 0.15, height = 0.1, color = "#5C8374") +
geom_smooth(method = "lm", se = TRUE, color = "#183D3D") +
labs(
title = "Model: Squirrel Density Predicting Boldness",
x = "Squirrel Density (per hectare)",
y = "Boldness Score"
) +
theme_minimal()
summary(model_density)
tidy(model_density, conf.int = TRUE)
```
**Plot Interpretation:** The dark line in the middle of the plot is the regression line, and the gray shaded area around it is the confidence band that shows the range where the true line probably falls. This band is narrow here, so the model is fairly confident about where the line sits, even though the line itself is nearly flat. The dots above and below are not errors. They are jittered, meaning small random noise was added just so the thousands of overlapping points at -1, 0, and 1 become visible instead of forming three solid horizontal bars.
**Model Summary Interpretation:**
The estimated model equation given by my model is:
$$\hat{\text{boldness}} = -0.202 + 0.003 \cdot \text{squirrel\_density}$$
- The **intercept** (-0.202) is the predicted boldness score for a squirrel in a hectare with zero other squirrels. This means that a squirrel living alone in a hectare with zero other squirrels observed is predicted to have a boldness score of -0.202, so it leans to be slightly avoidant.
- The **slope** (0.003) tells me that for every one-unit increase in squirrel density, boldness is predicted to go up by 0.003 points on average.
- My **confidence interval** is (0.0001, 0.0056). Therefore, I am 95% confident that the true effect of adding one more squirrel to a hectare is somewhere between almost zero and 0.0056 points of boldness. Because this range does not include zero, I can say the relationship is statistically significant. This means that it is unlikely this small positive trend happened purely by chance.
- The **R-squared** is 0.0014. This means squirrel density explains only about 0.14% of the differences in boldness between squirrels. In plain terms, density barely matters at all for predicting boldness, even though there is a real pattern.
- The **Slope** is 0.003, p = 0.041. The p-value tells me the probability of seeing a slope this large (or larger) purely by random chance if there were truly no relationship at all. Since 0.041 is below the common cutoff of 0.05, I treat this as a real, though very small, effect.
This analysis tells me that squirrel density does have a statistically real relationship with boldness, but the effect is so small that it has almost no practical impact. A squirrel in the most crowded hectare (32 squirrels) is only predicted to be about 0.09 points bolder than one in the least crowded hectare (1 squirrel), which is a tiny shift on a scale that only ranges from -1 to 1.
**Important note:** Boldness only takes three possible values: -1, 0, and 1. Therefore it is not a smooth, continuous measurement like height or temperature, where the difference between 1.1 and 1.2 means the same thing as the difference between 5.1 and 5.2. I created boldness to be closer to an ordered category like fled, indifferent, or approached.
But the problem is that Linear regression assumes the outcome is continuous and that a one-unit change means the same thing no matter where it happens on the scale. Here, that assumption does not perfectly hold. Going from "indifferent" (0) to "approached" (1) might represent a very different kind of behavioral jump than going from "fled" (-1) to "indifferent" (0). The model treats both jumps as equally sized, which is a big simplification.
This is why the plot looks a little unusual, with all the dots clumped into three horizontal bands instead of being spread smoothly across the whole range. A normal scatterplot assumes points could fall anywhere along the y-axis. But mine cannot because boldness was never built to be a smooth number and was actually built by combining two yes/no behaviors into a score.
Because of this, the slope, intercept, R-squared, and p-value reported above should be read as a simplified summary of an average linear trend rather than a precise description of how boldness actually works. After doing a bit more research, I think more statistically correct approach for this kind of three-level outcome would be ordinal logistic regression, which is built specifically for ranked categories like these. However, I used linear regression becuase I wanted to showcase specifically what I learned in this course!
### One Categorical Explanatory Variable: Primary Fur Color
Next, I want to find out if a squirrel's fur color (by itself), predicts how bold it is. To do this with a categorical variable like fur color, R picks one category to be the "reference group" that everything else gets compared against. Since Gray is the most common fur color, R automatically uses it as the reference. This means the model will tell me how much bolder or less bold Cinnamon and Black squirrels are compared to Gray squirrels, instead of giving each color its own separate number.
My expected model equation is:
$$E[\text{boldness} \mid \text{fur color}] = \beta_0 + \beta_1 \cdot \mathbb{1}[\text{Cinnamon}] + \beta_2 \cdot \mathbb{1}[\text{Black}]$$
The notation 1[Cinnamon] is an indicator variable. It equals 1 if the squirrel is Cinnamon and 0 if it is not. The same goes for 1[Black]. Because Gray is the reference group, there is no indicator for Gray. It is already represented by the intercept.
```{r}
squirrels %>%
group_by(primary_fur_color) %>%
summarise(
mean_boldness = mean(boldness, na.rm = TRUE),
sd_boldness = sd(boldness, na.rm = TRUE),
n = n()
)
```
This table shows the actual average boldness for each fur color group before any modeling. Gray squirrels average -0.172, Cinnamon squirrels average -0.110, and Black squirrels average -0.252. These raw averages should match up with what the regression model produces, and this was a way I double checked if the model is working correctly.
```{r}
model_color <- lm(boldness ~ primary_fur_color, data = squirrels)
summary(model_color)
tidy(model_color, conf.int = TRUE)
```
The fitted model equation is:
$$\hat{\text{boldness}} = -0.172 + 0.062 \cdot \mathbb{1}[\text{Cinnamon}] + (-0.081) \cdot \mathbb{1}[\text{Black}]$$
This means that the predicted boldness for a Gray squirrel is -0.172.
- To get the predicted boldness for a Cinnamon squirrel, I add 0.062 to that number, making it -0.110.
- To get the predicted boldness for a Black squirrel, I add -0.081, which gives us -0.252.
These numbers match the raw averages in the table above and I took this as a good sign that the model is doing what it should !
**Intercept:** The intercept of -0.172 is the predicted average boldness for a Gray squirrel, since Gray is the reference group. Every other prediction in this model is built by starting at this number and adjusting up or down.
**Slope:**
- Gray squirrels do not have a slope of their own since they are the reference group. Their predicted boldness comes directly from the intercept, -0.172.
- The Cinnamon slope is 0.062 with p = 0.021. This means Cinnamon squirrels are predicted to be 0.062 points bolder than Gray squirrels on average. The p-value tells me there is only about a 2.1% chance of seeing a difference this large between Cinnamon and Gray squirrels if there were truly no difference at all in the real world. Since this is below the usual 0.05 cutoff, this result is considered statistically significant.
- The Black slope is -0.081 with p = 0.104. This means Black squirrels are predicted to be 0.081 points more avoidant than Gray squirrels on average. The p-value tells me there is about a 10.4% chance of seeing a difference this large just by random chance, even if Black and Gray squirrels behave identically. Since that is above 0.05, this result is not statistically significant, meaning there is not enough evidence to say Black squirrels are truly different from Gray squirrels.
**Confidence intervals:**
- Gray squirrels do not get their own confidence interval here since they are the baseline that Cinnamon and Black are both being compared against.
- For Cinnamon squirrels, the 95% confidence interval is (0.0096, 0.115). Therefore, I am 95% confident that the true boldness advantage of Cinnamon squirrels over Gray squirrels is somewhere between about 0.01 and 0.11 points. Since this entire range is above zero, it supports the conclusion that Cinnamon squirrels are genuinely bolder and not just bolder by chance in this particular sample.
- For Black squirrels, the confidence interval is (-0.178, 0.017). This range includes zero, which means I cannot confidently say Black squirrels are different from Gray squirrels at all. The true difference could be negative (more avoidant), positive (bolder), or essentially zero.
**R-squared:** The R-squared is 0.0029. This means fur color alone, across all three groups, explains only about 0.29% of all the variation in boldness among squirrels. So fur color is a very weak predictor on its own. Thus, most of what makes a squirrel bold or avoidant comes from something other than its fur color.
### Full Model: Multiple Linear Regression
The two models above each looked at one variable at a time, but they cannot tell me whether fur color truly creates differences in squirrels boldness, or if those differences are actually explained by something else. For example, if Cinnamon squirrels happen to be observed foraging more often than Gray squirrels, their higher boldness might really be a foraging effect rather than something about their fur color specifically. To separate out these overlapping explanations, I built a multiple linear regression model that includes `fur color` and also `squirrel_density`, `age`, `shift`, and `foraging`. This helps account for possible confounding variables that might be tangled up with both fur color and boldness.
```{r}
model_full <- lm(
boldness ~ squirrel_density + primary_fur_color + age + shift + foraging,
data = squirrels
)
summary(model_full)
tidy(model_full, conf.int = TRUE)
```
In this multiple regression, I interpreted every slope "holding all else constant," by describing the effect of one variable while imagining that every other variable in the model stays fixed at some value. In this context, "all else" refers to squirrel density, fur color, age, shift, and foraging status all together.
The **intercept** of -0.220 represents the predicted boldness of an Adult, Gray squirrel observed in the AM time, not foraging, in a hectare with zero other squirrels.
- Holding all else constant, each additional squirrel observed in a hectare is associated with a 0.002 point increase in boldness. The standard error here (0.0014) is similar in size to the estimate itself, and the p-value of 0.081 is above the 0.05 cutoff, so this effect is not statistically significant in the full model.
- Holding all else constant, Cinnamon squirrels are predicted to be 0.057 points bolder than Gray squirrels. The standard error (0.027) is small relative to the estimate, and the p-value of 0.035 is below 0.05, so this effect remains statistically significant even after adding all the controls. This is a meaningful finding because it shows the Cinnamon boldness advantage is not simply explained away by age, time of day, foraging, or local crowding.
- Holding all else constant, Black squirrels are predicted to be 0.065 points more avoidant than Gray squirrels. The standard error (0.050) is large relative to the estimate, and the p-value of 0.197 is well above 0.05, so this effect is not statistically significant. This matches the earlier categorical model, where the Black coefficient also failed to reach significance, likely because there are only 103 Black squirrels in the dataset.
- Holding all else constant, Juvenile squirrels are predicted to be 0.020 points more avoidant than Adult squirrels. The standard error (0.029) is larger than the estimate itself, and the p-value of 0.489 is far above 0.05, so this effect is not statistically significant. There is not enough evidence to say age affects boldness once the other variables are accounted for.
- Holding all else constant, foraging is associated with a 0.072 point increase in boldness. The standard error (0.019) is small relative to the estimate, and the p-value is extremely small (p < 0.001), making this the strongest and most statistically significant effect in the entire model. This makes intuitive sense: a squirrel that is already searching for food is naturally more likely to approach a human as a potential food source.
- Holding all else constant, PM sightings are associated with a 0.026 point decrease in boldness compared to AM sightings, but the p-value of 0.166 means this is not statistically significant.
Comparing the Cinnamon coefficient across models is informative was really helpful for me to gain some helpful points.
- In the simple categorical model, Cinnamon squirrels were 0.062 points bolder than Gray.
- In the full model, that estimate barely moved to 0.057, and it remained statistically significant in both.
I think this stability is a good sign becuase it suggests that the Cinnamon effect is not just a byproduct of Cinnamon squirrels happening to be more often juvenile, more often foraging, or living in different density hectares. It looks like a real and independent relationship between fur color and boldness.
The R-squared of the full model is 0.0101. This means that all six variables together explain only about 1% of the variation in boldness. This is still a small number, but it is larger than the R-squared from the fur-color-only model (0.0029). This shows that adding foraging, density, age, and shift does meaningfully improve the model, even though boldness is still mostly unexplained by the variables collected in this dataset.
## Model Evaluation
It is important to check whether a model is correct before trusting its results becuase a model can show a statistically significant coefficient and still be built on assumptions that do not hold for the data. The LINE conditions, Linearity, Independence, Normality, and Equal variance, are the standard checks for whether a linear regression model is appropriate. If these conditions are badly violated, the p-values and confidence intervals from the model cannot be fully trusted.
```{r}
model_full %>%
ggplot(aes(x = .fitted, y = .resid)) +
geom_point(alpha = 0.2, color = "#5C8374") +
geom_hline(yintercept = 0, color = "black") +
geom_smooth(se = FALSE, color = "tomato") +
labs(
title = "Residuals vs. Fitted Values",
x = "Fitted Values",
y = "Residuals"
) +
theme_minimal()
```
**Linearity:** The residual plot shows three distinct diagonal bands instead of a random, formless scatter around zero. This pattern happens because boldness only takes three values (-1, 0, 1). Since a residual is just the observed value minus the predicted value, and the predicted value changes continuously while the observed value is locked to one of three numbers, each band traces out a straight diagonal line as the fitted value shifts. This is a clear violation of linearity, but it is a direct result of how the outcome variable was built, not a mistake in the model itself.
**Independence:** Each row represents a different squirrel sighting recorded by a different observer moving through the park's hectare grid, so observations are largely independent of each other. However, squirrels sighted in the same hectare may share some behavioral similarity due to a shared environment, which is a mild violation I think is worth noting.
**Normality:** Because boldness is between -1 and 1 and can only take three values, the residuals are split into three separate diagonal lines rather than forming one smooth bell-shaped spread. This is a clear violation of normality. With nearly 3,000 observations, the Central Limit Theorem offers some protection for the validity of the p-values, but normality is still technically violated here.
**Equal Variance:** The three bands in the plot are roughly parallel to each other and evenly spaced, which actually suggests the spread within each band stays fairly consistent across fitted values. This is a less severe violation than expected, though the three-banded structure itself still makes it hard to call this true equal variance in the traditional sense, since residuals are not scattered randomly above and below zero at any given fitted value.
**Addressing these limitations in future work:** All of these patterns trace back to the same root cause: boldness is an ordinal categorical variable with three levels, not a continuous measurement, even though it was coded as numbers for this project. A more appropriate model for future work would be ordinal logistic regression, which is built specifically for outcomes with ordered categories like fled, indifferent, and approached. This method is beyond the scope of this course, so linear regression was used here as a reasonable simplification, but the LINE violations above should be read as a result of that simplification rather than evidence that the data or model were handled incorrectly.
## Predictions and Conclusions
Using the full model, I predict the boldness of an Adult, Gray squirrel observed in the AM, not foraging, in a hectare with a squirrel density of 25, which represents a fairly crowded part of the park.
```{r}
new_quant <- data.frame(
squirrel_density = 25,
primary_fur_color = factor("Gray", levels = c("Gray", "Cinnamon", "Black")),
age = factor("Adult", levels = c("Adult", "Juvenile")),
shift = factor("AM"),
foraging = FALSE
)
predict(model_full, newdata = new_quant, interval = "prediction")
```
For an Adult, Gray squirrel in a hectare with 25 other squirrels, observed in the AM and not foraging, the model predicts a boldness score of -0.158. The 95% prediction interval is (-1.123, 0.807), meaning I am 95% confident that an individual squirrel matching this exact profile would have a boldness score somewhere in that wide range. This prediction is close to 0, which matches what I found earlier: most squirrels in this dataset are indifferent toward humans, and crowding by itself does very little to change that.
**Predicting using the categorical variable**
Next, I predict boldness for the exact same kind of squirrel, but now Cinnamon instead of Gray, to see how much fur color alone shifts the prediction.
```{r}
new_cat <- data.frame(
squirrel_density = 25,
primary_fur_color = factor("Cinnamon", levels = c("Gray", "Cinnamon", "Black")),
age = factor("Adult", levels = c("Adult", "Juvenile")),
shift = factor("AM"),
foraging = FALSE
)
predict(model_full, newdata = new_cat, interval = "prediction")
```
Comparing these two predictions isolates the effect of fur color directly. The Cinnamon prediction is -0.101, while the Gray prediction is -0.158. The difference between them is 0.057, the same number as the Cinnamon coefficient in the full model. In plain terms, this means that if I take an identical squirrel and only change its fur color from Gray to Cinnamon, holding density, age, shift, and foraging fixed, the model predicts a boldness score 0.057 points higher. 0.057 points higher.
One thing worth noting: both prediction intervals are very wide, stretching from about -1.1 to 0.85, almost the entire possible range of boldness. This is consistent with the model's low R-squared. The model can confidently estimate the average difference between groups, like the 0.057 Cinnamon effect, but it cannot confidently predict any individual squirrel's exact behavior, since so much of an individual squirrel's boldness is driven by factors outside this model. This would be worth a one-sentence mention in your key takeaways if you want to be thorough.
**Key takeaways**
I think most Central Park squirrels are indifferent to humans because the boldness distribution sharply piled up at 0 rather than spread evenly across approaching or fleeing.
- Fur color does matter, but only for one group: Cinnamon squirrels are reliably bolder than Gray squirrels, and this difference held up even after controlling for age, time of day, foraging, and local density.
- Black squirrels could not be reliably distinguished from Gray squirrels, most likely because there were only 103 Black squirrels in the dataset compared to over 2,000 Gray squirrels, making it harder to detect a real difference even if one exists.
- Foraging behavior turned out to be the single strongest predictor of boldness in the entire project, even stronger than fur color, which makes sense since a squirrel already searching for food has an obvious reason to approach a person.
Finally, the prediction intervals from the full model were very wide, stretching across almost the entire possible range of boldness from about -1.1 to 0.85. This tells me the model can confidently estimate average differences between groups, like the 0.057 point Cinnamon effect, but it cannot confidently predict how any single squirrel will behave, since most of what drives an individual squirrel's boldness comes from something outside this model.
**Answering the research question**
My research question asked whether a squirrel's fur color predicts its behavioral personality profile toward humans. Based on this analysis, I say the the answer is that fur color partially predicts squirrel behavior. Cinnamon squirrels show a small but statistically significant boldness advantage over Gray squirrels, and this advantage survives even after controlling for several other explanations, holding steady at about 0.057 points whether predicting between groups or comparing two individual squirrels that differ only in fur color. However, the full model's R-squared was only about 0.01, meaning all the variables in this project together explain just 1% of the variation in boldness. This tells me that while fur color is a real, measurable piece of the puzzle, it is a very small piece. Most of what determines whether an individual squirrel is bold or skittish around humans comes from something this dataset does not capture, possibly individual personality, past experiences with people, or specific locations within the park that were not measured here.
