Data Understanding
- for language enrollment rates prediction model
Doris Chen
2024-04-19
# Load packages
library(cowplot)
library(tidyverse)
library(corrplot)
library(leaps)Data processing please check “1_Data_Processing.Rmd”
# Load the data
all_data <- read_csv("Data/all_data.csv")
totl_lan <- read_csv("Data/All_Lang_Enrol.csv")
college_enrl <- read_csv("Data/UsCollege_Enrol.csv")
summary(all_data)## Year Language Lang_Enrol Lang_Enrol_Pct
## Min. :1958 Length:328 Min. : 0.0 Min. :0.00000
## 1st Qu.:1981 Class :character 1st Qu.: 727.8 1st Qu.:0.00070
## Median :1994 Mode :character Median : 2389.5 Median :0.00170
## Mean :1993 Mean :13147.6 Mean :0.01004
## 3rd Qu.:2008 3rd Qu.:14369.5 3rd Qu.:0.01323
## Max. :2021 Max. :72379.0 Max. :0.05550
##
## Pop Pop_UsPct Lang_Pop Lang_Pop_UsPct
## Min. : 74264 Min. :0.000400 Min. : 31614 Min. :0.00020
## 1st Qu.: 606186 1st Qu.:0.002500 1st Qu.: 343175 1st Qu.:0.00150
## Median : 895283 Median :0.003550 Median : 736797 Median :0.00290
## Mean :1316853 Mean :0.004671 Mean : 973738 Mean :0.00337
## 3rd Qu.:1723252 3rd Qu.:0.006025 3rd Qu.:1382340 3rd Qu.:0.00450
## Max. :4404678 Max. :0.013400 Max. :3348155 Max. :0.01020
## NA's :36 NA's :36
## Eng_LesVw_Pop Eng_LesVw_Pop_UsPct HEdu_Pop HEdu_UsPct
## Min. : 44146 Min. :0.00020 Min. : 4343 Min. :0.000600
## 1st Qu.: 204825 1st Qu.:0.00080 1st Qu.: 90257 1st Qu.:0.004600
## Median : 377164 Median :0.00130 Median : 258221 Median :0.006300
## Mean : 448984 Mean :0.00156 Mean : 423966 Mean :0.008906
## 3rd Qu.: 543873 3rd Qu.:0.00202 3rd Qu.: 552831 3rd Qu.:0.012025
## Max. :1492502 Max. :0.00460 Max. :2336589 Max. :0.030100
## NA's :76 NA's :76
## HEdu_Pct HEdu_Ratio Med_Inc Inc_Ratio
## Min. :0.0550 Min. :0.760 Min. : 2817 Min. :0.760
## 1st Qu.:0.2567 1st Qu.:1.620 1st Qu.: 19928 1st Qu.:1.067
## Median :0.3925 Median :1.800 Median : 42842 Median :1.235
## Mean :0.3904 Mean :1.840 Mean : 44165 Mean :1.243
## 3rd Qu.:0.5050 3rd Qu.:2.212 3rd Qu.: 63466 3rd Qu.:1.383
## Max. :0.7680 Max. :3.300 Max. :141906 Max. :2.040
##
## Pvt_Pct Pvt_Ratio Intl_Stud Intl_Stud_Pct
## Min. :0.0260 Min. :0.2500 Min. : 0 Min. :0.00000
## 1st Qu.:0.0660 1st Qu.:0.5700 1st Qu.: 3224 1st Qu.:0.01100
## Median :0.0980 Median :0.8500 Median : 13474 Median :0.03550
## Mean :0.1084 Mean :0.9129 Mean : 35901 Mean :0.05625
## 3rd Qu.:0.1340 3rd Qu.:1.1200 3rd Qu.: 42587 3rd Qu.:0.08125
## Max. :0.3550 Max. :3.0600 Max. :372532 Max. :0.34700
##
## Total_Trade Trade_Pct Total_Stud Presidency
## Min. : 0 Min. :0.00000 Min. : 3242000 Length:328
## 1st Qu.: 3305 1st Qu.:0.00640 1st Qu.:12205679 Class :character
## Median : 17709 Median :0.01575 Median :14358953 Mode :character
## Mean : 68730 Mean :0.03786 Mean :14352656
## 3rd Qu.: 78101 3rd Qu.:0.05158 3rd Qu.:18658756
## Max. :658795 Max. :0.18380 Max. :21019438
##
## Party_Gov Lan_UsPct
## Length:328 Min. :0.0000000
## Class :character 1st Qu.:0.0000575
## Mode :character Median :0.0001400
## Mean :0.0008159
## 3rd Qu.:0.0010300
## Max. :0.0037600
## It is a time series data set, spanning from 1958 to 2021. We will take a look at the trend of language enrollment across the entire data frame first.
# Check the trend of each language course enrollments
LanEnr1 <- ggplot(data = totl_lan, aes(x = factor(Year), y = Totl_Lan, group = 1))+
geom_point()+
geom_line()+
theme_classic()+
theme(axis.text.x = element_text(angle = 90))+
labs(title = "Total Foreign Language Courses Enrollments in US Institution", x = "Year", y = "Enrollment")
LanEnr2 <- ggplot(data = all_data, aes(x = as.factor(Year), y = Lang_Enrol, group = Language, col = Language))+
geom_point()+
geom_line()+
theme_classic()+
theme(axis.text.x = element_text(angle = 90))+
theme(legend.position = 'none')+
labs(title = "Asian Language Courses Enrollments", x = "Year", y = "Enrollment")
LanEnr3 <- ggplot(data = all_data, aes(x = as.factor(Year), y = Lang_Enrol_Pct, group = Language, col = Language))+
geom_point()+
geom_line()+
theme_classic()+
theme(axis.text.x = element_text(angle = 90))+
scale_y_continuous(labels = scales::percent)+
theme(legend.position = "bottom")+
theme(legend.title = element_blank())+
guides(col = guide_legend(nrow = 1))+
labs(title = "Asian Language Courses Enrollments Percentage", x = "Year", y = "Enrollment Percentage")
plot_grid(LanEnr1, LanEnr2, LanEnr3, align = "v", nrow = 3, rel_heights = c(0.3, 0.3, 0.4))
As can be seen from the plots above, the enrollments in Asian language courses roughly follow the general trend of total language enrollments with some deviations: Japanese enrollments increased rapidly in the 1980s and remain the highest among the six languages; Korean enrollments have risen slowly and steadily, even in the last decade when total language enrollments are declining; Filipino, Vietnamese, and Hindi have lower enrollments compared to the other three languages, and not showing the trend clearly enough in the plot. Let’s check the three languages alone.
ggplot(data = all_data[all_data$Language %in% c("Vietnamese", "Filipino", "Hindi"), ],
aes(x = as.factor(Year), y = Lang_Enrol_Pct, group = Language, col = Language))+
geom_point()+
geom_line()+
theme_classic()+
theme(axis.text.x = element_text(angle = 90))+
scale_y_continuous(labels = scales::percent)+
theme(legend.position = "bottom")+
theme(legend.title = element_blank())+
guides(col = guide_legend(nrow = 1))+
labs(x = "Year", y = "Enrollment Percentage")
The percentages of these three language course enrollments over the
total language courses enrollments are going up, similar trend as the
other three Asian languages.
Now we will compare the language courses enrollments with the total college enrollments in US.
totl_lan <- left_join(totl_lan, college_enrl, by = "Year")
totl_lan <- totl_lan %>%
mutate(ToLan_UsPct = round(Totl_Lan/Total_Stud, 4))lanPct1 <- ggplot(data = totl_lan, aes(x = factor(Year), y = ToLan_UsPct, group = 1, label = scales::percent(ToLan_UsPct)))+
geom_point()+
geom_line()+
scale_y_continuous(labels = scales::percent)+
theme_classic()+
geom_text(vjust = -0.6)+
labs(title = "All foreign language courses enrollments over total college enrollments in US",
x = "Year", y = "Percentage")
lanPct2 <- ggplot(data = all_data, aes(x = factor(Year), y = Lan_UsPct, group = Language, col = Language))+
geom_point()+
geom_line()+
scale_y_continuous(labels = scales::percent)+
theme_classic()+
theme(axis.text.x = element_text(angle = 90))+
theme(legend.position = "bottom", legend.title = element_blank())+
guides(col = guide_legend(nrow = 1))+
labs(title = "Each language courses enrollments over total college enrollments in US",
x = "Year", y = "Percentage")
plot_grid(lanPct1, lanPct2, align = "v", nrow = 2)
The trends are almost reverse when comparing with the total students in
US colleges.The percentage of total foreign language course enrollment
over total students in US colleges is going down since 1965 while the
percentage of each Asian language enrollments are going up, especially
Japanese in 1980s and Chinese in early 2000. Korean enrollment
percentage is going up steadily since 1980s.
# Check the distribution of "Lan_UsPct"
ggplot(data = all_data, aes(x = Lan_UsPct))+
geom_histogram(aes(y = after_stat(density)), bins = 100, col = "white", fill = "lightblue")+
geom_density()+
theme_classic()
The variable “Lan_UsPct” is percentage data with the range of (0, 1),
which should follow a beta distribution.
The histogram above looks like a beta distribution with very small \(\alpha\) and very large \(\beta\). The imperfect shape might be due to the small size of the data we have.
# Subset potential features for model building
all_ratios <- all_data[, c(1, 2, 24, 25, grep("Pct|Ratio", colnames(all_data)))]
# Multiply all percentage data with 100
for (i in grep("Pct", colnames(all_ratios))){
all_ratios[, i] = all_ratios[, i] * 100
}
# All ratio data subtract with 1
for (i in grep("Ratio", colnames(all_ratios))){
all_ratios[, i] = all_ratios[, i] - 1
}# Randomly extract 80% of the data for further data exploration and feature selection.
set.seed(123)
exploring <- all_ratios[sample(floor(nrow(all_ratios)*0.8)), ]Now, let’s check the correlation among the variables. First, Chi-square test of independence between the two categorical political variables under this hypothesis: H_0: The two variables are independent; H_1: the two variables are dependent; alpha = 0.05
chisq.test(exploring$Presidency, exploring$Party_Gov)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: exploring$Presidency and exploring$Party_Gov
## X-squared = 42.524, df = 1, p-value = 6.98e-11The p-value is very small, we reject the H_0, meaning the two variables are related to each other.
g1 <- ggplot(data = exploring, aes(x = Presidency, y = Lan_UsPct))+
geom_boxplot()
g2 <- ggplot(data = exploring, aes(x = Party_Gov, y = Lan_UsPct))+
geom_boxplot()
plot_grid(g1, g2, nrow = 1)
The plots do not look like there are much difference between the groups with respect to the percentage of language course enrollment. We will run anova test to double confirm: H_0: There are no differences across groups; H_1: There are some difference across groups; alpha = 0.05.
anova(aov(Lan_UsPct ~ Presidency, data = exploring))## Analysis of Variance Table
##
## Response: Lan_UsPct
## Df Sum Sq Mean Sq F value Pr(>F)
## Presidency 1 0.0091 0.0091155 0.6204 0.4316
## Residuals 260 3.8201 0.0146929anova(aov(Lan_UsPct ~ Party_Gov, data = exploring))## Analysis of Variance Table
##
## Response: Lan_UsPct
## Df Sum Sq Mean Sq F value Pr(>F)
## Party_Gov 1 0.0001 0.0001457 0.0099 0.9208
## Residuals 260 3.8291 0.0147274All the p-values are much greater than 0.05. We fail to reject H_0, meaning there are no differences across groups. We will exclude the two political variables from the model design.
We’ll move on to check the correlation coefficient among all possible pairs of the numeric data related to percentage or ratio.
corrplot(cor(exploring[, grep("Pct|Ratio", colnames(exploring))], use = "pairwise.complete.obs"), method = "number")
From the corrplot above,we can see that a high correlation coefficient(0.99) between the two potential response variable Lan_UsPct and Lang_Enrol_Pct. Five variables(Eng_LesVw_Pop_UsPct, HEdu_Ratio, Pvt_Pct, Intl_Stud_Pct, and Trade-Ratio) have correlation over 0.2 with both Lan_UsPct and Lang_Enrol_Pct.
The correlation between most variables are below 0.8, except for between Pop_UsPct, HEdu_UsPct, Lang_Pop_UsPct and Eng_LesVw_Pop_UsPct.
Note that, this correlation is between data in the same year, and the prediction model we are going to build will be base on historical data.
We will run linear models to select the best feature combinations. Because the availability of three population data (Pop_UsPct, Lang_Pop_UsPct and Eng_LesVw_Pop_UsPct) are different (Pop_UsPct 1958-2021, Lang_Pop_UsPct 1970-2021 and Eng_LesVw_Pop_UsPct 1980-2021), and since these three variables are highly correlated, we will combine them with other predictors separately.
rs58 <- summary(regsubsets(Lan_UsPct ~ Pop_UsPct + HEdu_Pct + Inc_Ratio + Pvt_Pct
+ Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring))
AIC58 = 2 * (2:8) + nrow(exploring) * log(rs58$rss/nrow(exploring))
BIC58 = log(nrow(exploring)) * (2:8) + nrow(exploring) * log(rs58$rss/nrow(exploring))
par(mfrow = c(1, 3))
plot(AIC58 ~ I(1:7), xlab = "Number of Predictors", ylab = "AIC")
plot(BIC58 ~ I(1:7), xlab = "Number of Predictors", ylab = "BIC")
plot(1:7, rs58$adjr2, xlab = "Number of Predictors", ylab = "Adjusted R-Squared")
Based on the plots above, the lowest value of AIC and BIC, and the highest value of adjusted R-square value are all landed in five predictors.
rs58$which## (Intercept) Pop_UsPct HEdu_Pct Inc_Ratio Pvt_Pct Pvt_Ratio Intl_Stud_Pct
## 1 TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## 3 TRUE FALSE TRUE TRUE FALSE FALSE FALSE
## 4 TRUE FALSE TRUE TRUE FALSE TRUE FALSE
## 5 TRUE FALSE TRUE TRUE FALSE TRUE TRUE
## 6 TRUE FALSE TRUE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## Trade_Pct
## 1 TRUE
## 2 TRUE
## 3 TRUE
## 4 TRUE
## 5 TRUE
## 6 TRUE
## 7 TRUErs58$adjr2## [1] 0.4804787 0.5162760 0.6056830 0.6830008 0.7115034 0.7110906 0.7107911summary(lm(Lan_UsPct ~ HEdu_Pct + Inc_Ratio + Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring))##
## Call:
## lm(formula = Lan_UsPct ~ HEdu_Pct + Inc_Ratio + Pvt_Ratio + Intl_Stud_Pct +
## Trade_Pct, data = exploring)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.153145 -0.041465 -0.004607 0.052665 0.128720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1565168 0.0159077 -9.839 < 2e-16 ***
## HEdu_Pct 0.0073184 0.0005344 13.693 < 2e-16 ***
## Inc_Ratio -0.4529443 0.0363954 -12.445 < 2e-16 ***
## Pvt_Ratio -0.1197629 0.0206021 -5.813 1.82e-08 ***
## Intl_Stud_Pct -0.0052086 0.0010139 -5.137 5.53e-07 ***
## Trade_Pct 0.0220685 0.0012299 17.943 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06506 on 256 degrees of freedom
## Multiple R-squared: 0.717, Adjusted R-squared: 0.7115
## F-statistic: 129.7 on 5 and 256 DF, p-value: < 2.2e-16The five predictors are significant in the summary above, the Adjusted R-squared is 0.7115.
So, here is our first feature group: Lan_UsPct ~ HEdu_Pct + Inc_Ratio + Pvt_Ratio + Intl_Stud_Pct + Trade_Pct
Let’s do the same for another two data combinations.
rs70 <- summary(regsubsets(Lan_UsPct ~ Lang_Pop_UsPct + HEdu_Pct + Inc_Ratio + Pvt_Pct
+ Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring[exploring$Year >= 1970, ]))
n <- nrow(exploring[exploring$Year >= 1970, ])
AIC70 = 2 * (2:8) + n * log(rs70$rss / n)
BIC70 = log(n) * (2:8) + n * log(rs70$rss / n)
par(mfrow = c(1, 3))
plot(AIC70 ~ I(1:7), xlab = "Number of Predictors", ylab = "AIC")
plot(BIC70 ~ I(1:7), xlab = "Number of Predictors", ylab = "BIC")
plot(1:7, rs70$adjr2, xlab = "Number of Predictors", ylab = "Adjusted R-Squared")
The lowest value of AIC and BIC landed in different number of predictors(AIC in 6, BIC in 5), while the highest value of Adjusted R-squared landed at the number of 6 predictors.
rs70$which## (Intercept) Lang_Pop_UsPct HEdu_Pct Inc_Ratio Pvt_Pct Pvt_Ratio Intl_Stud_Pct
## 1 TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## 3 TRUE FALSE TRUE TRUE FALSE FALSE FALSE
## 4 TRUE FALSE TRUE TRUE FALSE TRUE FALSE
## 5 TRUE FALSE TRUE TRUE FALSE TRUE TRUE
## 6 TRUE TRUE TRUE TRUE FALSE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## Trade_Pct
## 1 TRUE
## 2 TRUE
## 3 TRUE
## 4 TRUE
## 5 TRUE
## 6 TRUE
## 7 TRUEThe matrix above shows that “Lang_Pop_UsPct” was not included in the first five models until the last two models, which indicates the insignificance of the variable.
summary(lm(Lan_UsPct ~ Lang_Pop_UsPct + HEdu_Pct + Inc_Ratio
+ Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = all_ratios[all_ratios$Year >= 1970, ]))##
## Call:
## lm(formula = Lan_UsPct ~ Lang_Pop_UsPct + HEdu_Pct + Inc_Ratio +
## Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = all_ratios[all_ratios$Year >=
## 1970, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.176067 -0.049300 -0.009753 0.039747 0.213710
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0993491 0.0250057 -3.973 9.00e-05 ***
## Lang_Pop_UsPct 0.0337997 0.0256254 1.319 0.188231
## HEdu_Pct 0.0040450 0.0006374 6.346 8.65e-10 ***
## Inc_Ratio -0.1761841 0.0314080 -5.610 4.80e-08 ***
## Pvt_Ratio -0.0108634 0.0172572 -0.630 0.529525
## Intl_Stud_Pct -0.0058926 0.0015750 -3.741 0.000221 ***
## Trade_Pct 0.0228524 0.0016464 13.880 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07248 on 285 degrees of freedom
## Multiple R-squared: 0.6285, Adjusted R-squared: 0.6206
## F-statistic: 80.35 on 6 and 285 DF, p-value: < 2.2e-16The summary above confirms our conclusion. We will exclude this feature combination.
rs80 <- summary(regsubsets(Lan_UsPct ~ Eng_LesVw_Pop_UsPct + HEdu_Pct + Inc_Ratio + Pvt_Pct
+ Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring[exploring$Year >= 1980, ]))
n <- nrow(exploring[exploring$Year >= 1980, ])
AIC80 <- 2 * (2:8) + n * log(rs80$rss / n)
BIC80 <- log(n) * (2:8) + n * log(rs80$rss / n)
par(mfrow = c(1, 3))
plot(AIC80 ~ I(1:7), xlab = "Number of Predictors", ylab = "AIC")
plot(BIC80 ~ I(1:7), xlab = "Number of Predictors", ylab = "BIC")
plot(1:7, rs80$adjr2, xlab = "Number of Predictors", ylab = "Adjusted R-Squared")
All three plots indicates that the model with six predictors is the
best.
rs80$which## (Intercept) Eng_LesVw_Pop_UsPct HEdu_Pct Inc_Ratio Pvt_Pct Pvt_Ratio
## 1 TRUE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE FALSE TRUE FALSE FALSE
## 3 TRUE FALSE FALSE TRUE TRUE FALSE
## 4 TRUE FALSE TRUE TRUE FALSE TRUE
## 5 TRUE FALSE TRUE TRUE FALSE TRUE
## 6 TRUE TRUE TRUE TRUE FALSE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE
## Intl_Stud_Pct Trade_Pct
## 1 FALSE TRUE
## 2 FALSE TRUE
## 3 FALSE TRUE
## 4 FALSE TRUE
## 5 TRUE TRUE
## 6 TRUE TRUE
## 7 TRUE TRUEsummary(lm(Lan_UsPct ~ Eng_LesVw_Pop_UsPct + HEdu_Pct + Inc_Ratio
+ Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring[exploring$Year >= 1980, ]))##
## Call:
## lm(formula = Lan_UsPct ~ Eng_LesVw_Pop_UsPct + HEdu_Pct + Inc_Ratio +
## Pvt_Ratio + Intl_Stud_Pct + Trade_Pct, data = exploring[exploring$Year >=
## 1980, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.157781 -0.047181 -0.006834 0.054719 0.130483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2865073 0.0382096 -7.498 2.89e-12 ***
## Eng_LesVw_Pop_UsPct 0.1859675 0.0564396 3.295 0.00119 **
## HEdu_Pct 0.0098011 0.0009771 10.031 < 2e-16 ***
## Inc_Ratio -0.5390323 0.0476730 -11.307 < 2e-16 ***
## Pvt_Ratio -0.1902271 0.0332520 -5.721 4.38e-08 ***
## Intl_Stud_Pct -0.0079982 0.0016678 -4.796 3.40e-06 ***
## Trade_Pct 0.0248680 0.0018076 13.758 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06482 on 179 degrees of freedom
## Multiple R-squared: 0.7784, Adjusted R-squared: 0.771
## F-statistic: 104.8 on 6 and 179 DF, p-value: < 2.2e-16All six predictors are significant in the summary above, and the Adjusted R-squared in this model is 0.771, higher than the first linear model. We will choose the feature group with six predictors.
In summary, we got two potential feature groups:
Model_58 (for all data since 1958): Lan_UsPct ~ HEdu_Pct + Inc_Ratio + Pvt_Ratio + Intl_Stud_Pct + Trade_Pct
Model_80 (for data since 1980): Lan_UsPct ~ Eng_LesVw_Pop_UsPct + HEdu_Pct + Inc_Ratio + Pvt_Ratio + Intl_Stud_Pct + Trade_Pct
Feature extraction
We will add the rolling average value of prior 3-5 years of all predictors;
# Extract required features
all_features <- all_ratios %>%
select(Year, Language, Lan_UsPct, Lang_Enrol_Pct, Eng_LesVw_Pop_UsPct, HEdu_Pct,
Inc_Ratio, Pvt_Ratio, Intl_Stud_Pct, Trade_Pct)# Add features of rolling average of 3-5 years into the data frame
col_list <- list()
l = 1
for (i in c(3:10)){
feature_list = list()
k = 1
for (lan in unique(all_features$Language)){
LanGr <- all_features[all_features$Language == lan, ]
for (j in c(3:5)){
LanGr[seq(j+1, nrow(LanGr), by = 1), ncol(LanGr)+1] = sapply(1:(nrow(LanGr)-j),
function(x){round(mean(LanGr[[i]][x:(x+j-1)]),3)})
colnames(LanGr)[ncol(LanGr)] = paste0(colnames(LanGr[i]), "_avg", j)
}
feature_list[[k]] = LanGr
k = k+1
}
col_list[[l]] = dplyr::bind_rows(feature_list)
l = l+1
}all_features <- col_list[[1]][, c(1:3, 11:13)]
for (i in c(2:length(col_list))){
all_features <- dplyr::left_join(all_features, col_list[[i]][, -c(3:10)], by = c("Year", "Language"))
}# Save the dataset
write_csv(all_features, "Data/all_features.csv")