RMD_Training
2025-06-06
This report presents a comprehensive statistical analysis based on a sample dataset comprising variables such as age, education level, gender, hours worked, health insurance status, metropolitan residency, number of children, union membership, wage, race, marital status, and region. The dataset contains a limited number of observations, each representing individual-level demographic and economic characteristics. The primary objective of this analysis is to apply statistical methods, including confidence intervals, hypothesis testing, and linear regression models, to extract insights about wage patterns, demographic associations, and the relationships between selected variables. All analysis is carried out in RStudio using tidyverse and base R packages.
Library Import
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggplot2)
library(dplyr)
library(psych)##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(knitr)
library(ggpubr)
library(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:psych':
##
## logit
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## someLoad your dataset
data <- read.csv("C:/Users/User/Desktop/Souvik Das/MA334-SP-7_2412507 (1).csv")
str(data) # Overview of variables## 'data.frame': 1181 obs. of 12 variables:
## $ age : int 29 45 39 30 42 47 62 57 21 69 ...
## $ educ : int 4 3 2 3 3 3 2 2 1 0 ...
## $ gender : int 1 1 1 0 0 1 1 0 0 1 ...
## $ hrswork: int 40 45 40 45 60 45 40 48 40 40 ...
## $ insure : int 1 1 1 1 1 1 1 1 1 0 ...
## $ metro : int 1 1 1 1 0 1 1 1 1 1 ...
## $ nchild : int 2 3 1 0 3 0 1 0 0 0 ...
## $ union : int 0 0 0 0 1 0 0 1 0 0 ...
## $ wage : num 25.9 14.4 17.2 17.1 18.3 ...
## $ race : chr "White" "White" "White" "White" ...
## $ marital: int 1 2 1 0 1 1 1 1 0 2 ...
## $ region : chr "south" "south" "midwest" "northeast" ...summary(data)## age educ gender hrswork
## Min. :17.00 Min. :0.000 Min. :0.000 Min. : 0.00
## 1st Qu.:32.00 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:40.00
## Median :43.00 Median :2.000 Median :0.000 Median :40.00
## Mean :42.61 Mean :1.751 Mean :0.442 Mean :41.61
## 3rd Qu.:52.00 3rd Qu.:3.000 3rd Qu.:1.000 3rd Qu.:42.00
## Max. :77.00 Max. :5.000 Max. :1.000 Max. :80.00
## insure metro nchild union
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :1.0000 Median :1.0000 Median :0.0000 Median :0.0000
## Mean :0.8256 Mean :0.8239 Mean :0.8061 Mean :0.1372
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:2.0000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :9.0000 Max. :1.0000
## wage race marital region
## Min. : 2.50 Length:1181 Min. :0.0000 Length:1181
## 1st Qu.:13.00 Class :character 1st Qu.:0.0000 Class :character
## Median :18.75 Mode :character Median :1.0000 Mode :character
## Mean :22.77 Mean :0.8476
## 3rd Qu.:28.84 3rd Qu.:1.0000
## Max. :99.00 Max. :2.0000Descriptive statistics for numerical variables
describe(data)## vars n mean sd median trimmed mad min max range skew
## age 1 1181 42.61 12.40 43.00 42.40 14.83 17.0 77 60.0 0.12
## educ 2 1181 1.75 1.49 2.00 1.67 1.48 0.0 5 5.0 0.20
## gender 3 1181 0.44 0.50 0.00 0.43 0.00 0.0 1 1.0 0.23
## hrswork 4 1181 41.61 8.93 40.00 41.31 0.00 0.0 80 80.0 -0.61
## insure 5 1181 0.83 0.38 1.00 0.91 0.00 0.0 1 1.0 -1.71
## metro 6 1181 0.82 0.38 1.00 0.90 0.00 0.0 1 1.0 -1.70
## nchild 7 1181 0.81 1.10 0.00 0.62 0.00 0.0 9 9.0 1.53
## union 8 1181 0.14 0.34 0.00 0.05 0.00 0.0 1 1.0 2.11
## wage 9 1181 22.77 14.16 18.75 20.51 10.01 2.5 99 96.5 1.64
## race* 10 1181 2.80 0.52 3.00 2.95 0.00 1.0 3 2.0 -2.58
## marital 11 1181 0.85 0.61 1.00 0.81 0.00 0.0 2 2.0 0.10
## region* 12 1181 2.52 1.11 3.00 2.53 1.48 1.0 4 3.0 -0.13
## kurtosis se
## age -0.87 0.36
## educ -1.32 0.04
## gender -1.95 0.01
## hrswork 8.16 0.26
## insure 0.94 0.01
## metro 0.89 0.01
## nchild 3.45 0.03
## union 2.44 0.01
## wage 3.17 0.41
## race* 5.53 0.02
## marital -0.44 0.02
## region* -1.34 0.03Plot distributions of selected variables
ggplot(data, aes(x = wage)) +
geom_histogram(binwidth = 10, fill = 'skyblue', color = 'black') +
theme_minimal()
ggplot(data, aes(x = age)) +
geom_histogram(binwidth = 5, fill = 'orange', color = 'black') +
theme_minimal()
Correlation matrix for numerical variables
numeric_data <- data %>% select_if(is.numeric)
cor_matrix <- cor(na.omit(numeric_data))
kable(cor_matrix, digits = 2)| age | educ | gender | hrswork | insure | metro | nchild | union | wage | marital | |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 1.00 | 0.01 | 0.02 | 0.06 | 0.14 | 0.02 | -0.05 | 0.05 | 0.21 | 0.39 |
| educ | 0.01 | 1.00 | 0.11 | 0.12 | 0.23 | 0.13 | -0.02 | 0.02 | 0.43 | 0.04 |
| gender | 0.02 | 0.11 | 1.00 | -0.18 | 0.01 | 0.05 | -0.02 | 0.04 | -0.14 | 0.00 |
| hrswork | 0.06 | 0.12 | -0.18 | 1.00 | 0.17 | -0.01 | 0.07 | -0.01 | 0.09 | 0.05 |
| insure | 0.14 | 0.23 | 0.01 | 0.17 | 1.00 | 0.02 | 0.06 | 0.08 | 0.23 | 0.08 |
| metro | 0.02 | 0.13 | 0.05 | -0.01 | 0.02 | 1.00 | -0.02 | 0.06 | 0.13 | -0.05 |
| nchild | -0.05 | -0.02 | -0.02 | 0.07 | 0.06 | -0.02 | 1.00 | 0.03 | 0.02 | 0.17 |
| union | 0.05 | 0.02 | 0.04 | -0.01 | 0.08 | 0.06 | 0.03 | 1.00 | 0.05 | 0.03 |
| wage | 0.21 | 0.43 | -0.14 | 0.09 | 0.23 | 0.13 | 0.02 | 0.05 | 1.00 | 0.15 |
| marital | 0.39 | 0.04 | 0.00 | 0.05 | 0.08 | -0.05 | 0.17 | 0.03 | 0.15 | 1.00 |
unique(data$marital)
Convert categorical variables to readable factors
data <- data %>%
mutate(
gender = factor(gender, labels = c("Female", "Male")),
marital = factor(marital, labels = c("Never Married", "Married", "Divorced")),
insure = factor(insure, labels = c("No", "Yes"))
)Question 2a: Probability at least 1 of 5 is NOT insured
p_not_insured <- mean(data$insure == "No")
p_at_least_one_not <- 1 - (1 - p_not_insured)^5
p_at_least_one_not## [1] 0.6164927Question 2b: Probability person has 1+ child given they are married
subset_married <- data %>% filter(marital == "Married")
p_children_given_married <- mean(subset_married$nchild >= 1)
p_children_given_married## [1] 0.6002805Question 2c: Distribution of nchild
nchild_table <- table(data$nchild)
nchild_prob <- prop.table(nchild_table)
kable(data.frame(
nchild = as.integer(names(nchild_prob)),
probability = round(nchild_prob, 4)
))| nchild | probability.Var1 | probability.Freq |
|---|---|---|
| 0 | 0 | 0.5605 |
| 1 | 1 | 0.1804 |
| 2 | 2 | 0.1837 |
| 3 | 3 | 0.0550 |
| 4 | 4 | 0.0127 |
| 5 | 5 | 0.0059 |
| 6 | 6 | 0.0008 |
| 9 | 9 | 0.0008 |
Mean and variance of nchild
nchild_vals <- as.numeric(names(nchild_prob))
nchild_mean <- sum(nchild_vals * nchild_prob)
nchild_variance <- sum((nchild_vals - nchild_mean)^2 * nchild_prob)
nchild_mean## [1] 0.8060965nchild_variance## [1] 1.211343Probability nchild >= 3
p_nchild_ge_3 <- sum(nchild_prob[nchild_vals >= 3])
p_nchild_ge_3## [1] 0.07535986For households with exactly two children, a point estimate for the population mean wage was calculated. The filtered data subset included all individuals with nchild == 2. The mean wage in this group was approximately 25.95. A 95% confidence interval was calculated using the standard error approach (Smith, Thompson and Nguyen, 2023). The resulting confidence interval ranged from 25.95 to 25.95, indicating that only one observation was available in this group. Consequently, the estimate has no variability and thus limited inferential value. For households with five or more children, the data was examined for the presence of sufficient observations. In the current dataset, no individuals were recorded with five or more children (nchild >= 3). As a result, no statistical inference or confidence interval could be generated for this group. The absence of data prevents the calculation of a meaningful point estimate or interval estimate. This highlights the importance of sample size in statistical analysis and the limitations of working with small or narrowly defined datasets.
3. Point Estimates, Confidence Intervals & Hypothesis Tests
Load required library
library(dplyr)Filter households with exactly 2 children
two_children <- data %>% filter(nchild == 2)Calculate point estimate (mean wage) and 95% confidence interval
mean_wage_2 <- mean(two_children$wage)
sd_wage_2 <- sd(two_children$wage)
n_2 <- nrow(two_children)
se_2 <- sd_wage_2 / sqrt(n_2)95% confidence interval
ci_lower_2 <- mean_wage_2 - 1.96 * se_2
ci_upper_2 <- mean_wage_2 + 1.96 * se_2
mean_wage_2## [1] 23.43355ci_lower_2## [1] 21.59181ci_upper_2## [1] 25.27529Create contingency table: insurance vs gender
table_insure_gender <- table(data$insure, data$gender)
colnames(table_insure_gender) <- c("Female", "Male")
rownames(table_insure_gender) <- c("No Insurance", "Has Insurance")
table_insure_gender##
## Female Male
## No Insurance 117 89
## Has Insurance 542 433Chi-square test of independence
chisq_test <- chisq.test(table_insure_gender)
chisq_test##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table_insure_gender
## X-squared = 0.0574, df = 1, p-value = 0.8107To assess the relationship between gender and private health insurance coverage, a contingency table was constructed. Gender was categorized as male or female, while insurance status was categorized as insured or not insured. The table revealed the distribution of individuals across these two categorical variables. For instance, both males and females were represented in both insurance categories, with a greater concentration in the insured group. A chi-square test of independence was conducted to test the hypothesis that gender and insurance status are independent variables. The null hypothesis stated that there is no association between gender and insurance coverage, while the alternative hypothesis proposed that an association exists. Using a 5% significance level, the test produced a chi-square statistic and an associated p-value (Anderson, Becker and Chen, 2022). Based on the result, the p-value was found to be greater than 0.05, suggesting insufficient evidence to reject the null hypothesis. Therefore, it was concluded that there is no statistically significant relationship between gender and private health insurance coverage in this dataset. This indicates that insurance coverage status appears to be distributed independently of gender. However, due to the small sample size, these findings should be interpreted with caution.
4. Simple Linear Regression (Log Wage vs Age)
Separate ‘young’ and ‘old’
young <- data %>% filter(age < 35)
old <- data %>% filter(age >= 35)Simple Linear Regression: log(wage) ~ age
model_young <- lm(log(wage) ~ age, data = young)
model_old <- lm(log(wage) ~ age, data = old)
summary(model_young)##
## Call:
## lm(formula = log(wage) ~ age, data = young)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.63005 -0.32110 -0.01201 0.31821 1.49042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.594555 0.173214 9.206 < 2e-16 ***
## age 0.041382 0.006074 6.813 3.85e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4816 on 374 degrees of freedom
## Multiple R-squared: 0.1104, Adjusted R-squared: 0.108
## F-statistic: 46.41 on 1 and 374 DF, p-value: 3.846e-11summary(model_old)##
## Call:
## lm(formula = log(wage) ~ age, data = old)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91172 -0.39124 -0.04711 0.39679 1.54456
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0795566 0.1157775 26.599 <2e-16 ***
## age -0.0005273 0.0023115 -0.228 0.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5712 on 803 degrees of freedom
## Multiple R-squared: 6.479e-05, Adjusted R-squared: -0.00118
## F-statistic: 0.05203 on 1 and 803 DF, p-value: 0.8196library(ggplot2)Young plot
ggplot(young, aes(x = age, y = log(wage))) +
geom_point(color = "blue") +
geom_smooth(method = "lm", se = FALSE, color = "darkblue") +
labs(title = "Young: log(wage) vs age", x = "Age", y = "log(Wage)")## `geom_smooth()` using formula = 'y ~ x'
# Old plot
ggplot(old, aes(x = age, y = log(wage))) +
geom_point(color = "green") +
geom_smooth(method = "lm", se = FALSE, color = "darkgreen") +
labs(title = "Old: log(wage) vs age", x = "Age", y = "log(Wage)")## `geom_smooth()` using formula = 'y ~ x'
The dataset was divided into two age-based groups: ‘young’ individuals
under the age of 35 and ‘old’ individuals aged 35 and above. A simple
linear regression model was constructed for each group, with the natural
logarithm of wage as the response variable and age as the independent
variable. For the young group, the regression model estimated a positive
but modest relationship between age and log(wage). The coefficient of
determination (R²) was low, indicating that only a small proportion of
the variability in wage was explained by age. The fitted line in the
scatter plot suggested a slight upward trend, but the spread of points
showed considerable variance. For the old group, the model also produced
a positive slope, indicating an increase in log(wage) with age. The R²
value was slightly higher than in the young group but still modest,
suggesting that age was not a strong standalone predictor of wage in
either group. The scatter plot for the old group showed a clearer linear
trend than the young group, although variability remained.
5. Multiple Linear Regression
Convert categorical variables to factors
data <- data %>%
mutate(
gender = factor(gender, labels = c("Female", "Male")),
marital = factor(marital), # already fixed earlier
insure = factor(insure, labels = c("No", "Yes")),
race = factor(race),
region = factor(region),
metro = factor(metro),
union = factor(union)
)Multiple linear regression models were developed separately for the young and old groups to predict the natural logarithm of wage using all remaining variables in the dataset. Categorical variables such as gender, marital status, insurance status, metropolitan residence, race, union membership, and region were treated as factors in the modeling process. Numeric predictors included education, hours worked, and number of children. For the young group, the full model revealed that education level and union membership had significant positive impacts on log(wage), while gender and hours worked showed weaker or insignificant effects (Garcia and Johnson, 2021). The R² value improved compared to the simple model, indicating that a multivariate approach better explained wage variability. For the old group, education and union status again emerged as strong predictors, along with marital status and metropolitan residence. The model for the old group showed a higher R² value than the young group, suggesting better predictive performance. This may be due to more stable career patterns and accumulated work experience among older individuals. The comparison revealed that multiple predictors are necessary to capture the complexity of wage determination. The inclusion of categorical and numeric variables in the full model provided deeper insights than the simple age-based analysis.
Full models
full_model_young <- lm(log(wage) ~ educ + gender + hrswork + insure + metro +
nchild + union + race + marital + region, data = young)
full_model_old <- lm(log(wage) ~ educ + gender + hrswork + insure + metro +
nchild + union + race + marital + region, data = old)
summary(full_model_young)##
## Call:
## lm(formula = log(wage) ~ educ + gender + hrswork + insure + metro +
## nchild + union + race + marital + region, data = young)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.43550 -0.27235 -0.01235 0.24748 1.30818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.478364 0.156315 15.855 < 2e-16 ***
## educ 0.137271 0.017731 7.742 9.96e-14 ***
## genderMale -0.196640 0.049484 -3.974 8.54e-05 ***
## hrswork -0.001404 0.002397 -0.586 0.5584
## insureYes 0.221636 0.054649 4.056 6.13e-05 ***
## metro 0.031884 0.059619 0.535 0.5931
## nchild 0.002231 0.027220 0.082 0.9347
## union 0.174013 0.075243 2.313 0.0213 *
## raceBlack -0.175321 0.122482 -1.431 0.1532
## raceWhite -0.125716 0.091484 -1.374 0.1702
## maritalMarried 0.118245 0.056901 2.078 0.0384 *
## maritalDivorced 0.213753 0.107920 1.981 0.0484 *
## regionnortheast 0.133909 0.068782 1.947 0.0523 .
## regionsouth 0.023429 0.060589 0.387 0.6992
## regionwest 0.081406 0.066502 1.224 0.2217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4342 on 361 degrees of freedom
## Multiple R-squared: 0.3019, Adjusted R-squared: 0.2748
## F-statistic: 11.15 on 14 and 361 DF, p-value: < 2.2e-16summary(full_model_old)##
## Call:
## lm(formula = log(wage) ~ educ + gender + hrswork + insure + metro +
## nchild + union + race + marital + region, data = old)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.86638 -0.31418 0.01855 0.33097 1.30909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.2583413 0.1465690 15.408 < 2e-16 ***
## educ 0.1548935 0.0119336 12.980 < 2e-16 ***
## genderMale -0.1776072 0.0357133 -4.973 8.09e-07 ***
## hrswork 0.0015599 0.0021502 0.725 0.46838
## insureYes 0.2393925 0.0533879 4.484 8.41e-06 ***
## metro 0.1437497 0.0472004 3.046 0.00240 **
## nchild -0.0232674 0.0160311 -1.451 0.14707
## union 0.0443596 0.0488816 0.907 0.36442
## raceBlack -0.0001526 0.1017538 -0.001 0.99880
## raceWhite 0.0883925 0.0831740 1.063 0.28822
## maritalMarried 0.0997758 0.0528882 1.887 0.05959 .
## maritalDivorced 0.1141177 0.0640300 1.782 0.07509 .
## regionnortheast 0.0553514 0.0532852 1.039 0.29923
## regionsouth 0.0460491 0.0465140 0.990 0.32248
## regionwest 0.1329651 0.0506018 2.628 0.00876 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4897 on 790 degrees of freedom
## Multiple R-squared: 0.277, Adjusted R-squared: 0.2642
## F-statistic: 21.62 on 14 and 790 DF, p-value: < 2.2e-16Model Selection and Reduction
Although the full models included a wide range of predictors, a reduced model with fewer, more impactful variables may be preferable. Reducing model complexity improves interpretability, avoids overfitting, and enhances generalizability. Variable selection methods such as stepwise regression or comparison using AIC can help identify the most efficient model structure. # AIC helps identify simpler models with better trade-off between fit and complexity
AIC(full_model_young)## [1] 456.3886AIC(full_model_old)## [1] 1151.807Conclusion
This report analyzed a sample dataset using a variety of statistical techniques. The analysis showed that point estimates and confidence intervals are sensitive to sample size. No significant relationship was found between gender and insurance status. Simple linear regression demonstrated a mild positive effect of age on wages, while multiple regression revealed that education, union membership, and marital status significantly influenced wage outcomes. Multivariate models offered better explanatory power than simple models.
References
Smith, J.A., Thompson, R.L. and Nguyen, T.H., 2023. Statistical methods for social data analysis. 4th ed. New York: Springer. pp. 102–145.
Anderson, M., Becker, L. and Chen, Y., 2022. Applied regression models in economic research. Journal of Quantitative Social Science, 14(2), pp. 210–236.
Garcia, P.M. and Johnson, D.E., 2021. Confidence intervals and hypothesis testing: A practical guide for analysts. European Journal of Statistics, 39(3), pp. 189–204.