I'm having more conceptual difficulty with how to formulate my model and which statistical method to choose to analyze it.
In this case, I need to analyze the resilience of the cultural sector in different neighborhoods of a city in the face of the impact of Covid-19. My interest is to compare neighborhoods that had a specific policy (treated group) with those that did not (control group).
I was thinking of applying a Diff-in-Diff model to analyze differences in variation in the number of cultural centers and in their income from treated versus untreated neighborhoods (communes). The problem is that the policy was not created as a response to the pandemic, it has been in existence for 10 years. I have data from 2018, 2019, and 2020 (I can request the data until 2014, but there's no newer data after 2020).
The idea is to understand whether the presence of the program correlated with the increase in resilience or not of these neighborhoods - how the number of cultural centers and their income after the first pandemic hit was affected.
I have the data from the following variables:
id - identification variable grouping entries for commune, type, income_range
treatment - a dummy variable for the communes treated or not
Commune - the number of the commune (from 1 to 15)
type - the type of the cultural center (library, Music Hall, etc.)
income_range - the income range of the centers (low, mid, high)
Number - number of the cultural centers of a type and income range in that commune in that year
Income - income of the cultural centers of a type and income range in that commune in that
Change_N - the variation of cultural centers compared to the last year’s (year x - year x-1)
Change_Inc - the variation of income compared to the last year’s (year x - year x-1)
And thought of the following model:
model <- lm(Change_N ~ Treatment + Number + Treatment * Number + income_range)
for the visualization I have
plot_data <- table %>%
mutate(Treatment = factor(Treatment)) %>%
mutate(Year = factor(Year)) %>%
group_by(Year, Treatment) %>%
summarize(mean_change = mean(Change_N),
se_change = sd(Change_N) / sqrt(n()),
upper = mean_change + (-1.96 * se_change),
lower = mean_change + (1.96 * se_change))
plot_1 <- ggplot(plot_data, aes(x = Year, y = mean_change, color = Treatment)) +
geom_pointrange(aes(ymin = lower, ymax = upper)) +
geom_line(aes(group = Treatment)) +
labs(x = "Year", y = "Mean Change")
I would like to know if the model is valid or if I know my stats better and search for another one lol. If so, any tips on how I can answer this question with the data I have? Also, if this type of visualization could be useful for what I'm thinking.