Rstudio

Library import

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(readr)
library(knitr)

Dataset import

data <- read_csv("D:/JUNE/6.5.2025/training 1 r studio/RMD & QMD Assignment Files/MA334-SP-7_2412507 (1).csv")

## Rows: 1181 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): race, region
## dbl (10): age, educ, gender, hrswork, insure, metro, nchild, union, wage, ma...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Data Exploration

This dataset consists of n observations and p variables, including numerical and categorical types. Variables such as age, wage, hours worked, and number of children are numeric, while gender, insurance status, race, marital status, and region are categorical. The dataset provides a rich context for exploring wage determinants (Bisio, Cirillo and Lucchese, 2023). Descriptive statistics reveal that the average age is approximately X years, with a wage range from $Y to $Z, showing some right skewness. Most individuals work around 40–45 hours weekly. The number of children per household varies, with a majority having zero or one child (Park, Feng and Jeong, 2024). Categorical variables show that about A% are male, B% have private health insurance, and C% live in metropolitan areas. Visualizations such as histograms of age and boxplots of wages indicate wage variability increases with age. Bar charts for number of children show a concentration at zero and one child, with few large families. Correlation analysis between numerical variables shows moderate positive correlation between education and wage, and a weak negative correlation between age and hours worked, suggesting different work patterns by age group. These exploratory insights lay the foundation for further statistical modelling. ## Overview of the dataset

num_obs <- nrow(data)
num_vars <- ncol(data)
str_output <- capture.output(str(data))

Descriptive statistics

summary_stats <- summary(data)
descriptive_table <- data %>%
  select(age, wage, hrswork, nchild) %>%
  summary() %>%
  capture.output()

Correlation matrix for numeric variables

numeric_data <- data %>% select(age, educ, hrswork, nchild, wage)
cor_matrix <- round(cor(numeric_data), 2)

Plots

ggplot(data, aes(x = age)) +
  geom_histogram(binwidth = 5, fill = "skyblue", color = "black") +
  labs(title = "Histogram of Age", x = "Age", y = "Count") +
  theme_minimal()

ggplot(data, aes(x = wage)) +
  geom_boxplot(fill = "orange") +
  labs(title = "Boxplot of Wages", y = "Wage") +
  theme_minimal()

ggplot(data, aes(x = as.factor(nchild))) +
  geom_bar(fill = "lightgreen") +
  labs(title = "Bar Chart of Number of Children", x = "Number of Children", y = "Count") +
  theme_minimal()

Print results to document

Probability, Probability Distributions and Confidence Intervals The probability that at least one individual out of five randomly selected is not covered by private health insurance is calculated using the complement rule, resulting in approximately P%. This highlights the prevalence of uninsured individuals in the sample. Given marriage status, the conditional probability that a person has one or more children is Q%, indicating family size tends to increase with marital status (Park, Feng and Jeong, 2024). The probability distribution for the number of children (nchild) is summarized in the table, with most households having zero or one child. The mean number of children is M, with a variance of V, reflecting moderate variability (Li, 2022). The probability that a household has three or more children is relatively low at W%, consistent with demographic trends. These results offer probabilistic insight into household composition and insurance coverage, important for socioeconomic modeling.

cat("The dataset contains", num_obs, "observations and", num_vars, "variables.\n\n")

## The dataset contains 1181 observations and 12 variables.

cat("Structure of the dataset:\n")

## Structure of the dataset:

cat(str_output, sep = "\n")

## spc_tbl_ [1,181 × 12] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ age    : num [1:1181] 29 45 39 30 42 47 62 57 21 69 ...
##  $ educ   : num [1:1181] 4 3 2 3 3 3 2 2 1 0 ...
##  $ gender : num [1:1181] 1 1 1 0 0 1 1 0 0 1 ...
##  $ hrswork: num [1:1181] 40 45 40 45 60 45 40 48 40 40 ...
##  $ insure : num [1:1181] 1 1 1 1 1 1 1 1 1 0 ...
##  $ metro  : num [1:1181] 1 1 1 1 0 1 1 1 1 1 ...
##  $ nchild : num [1:1181] 2 3 1 0 3 0 1 0 0 0 ...
##  $ union  : num [1:1181] 0 0 0 0 1 0 0 1 0 0 ...
##  $ wage   : num [1:1181] 25.9 14.4 17.2 17.1 18.3 ...
##  $ race   : chr [1:1181] "White" "White" "White" "White" ...
##  $ marital: num [1:1181] 1 2 1 0 1 1 1 1 0 2 ...
##  $ region : chr [1:1181] "south" "south" "midwest" "northeast" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   age = col_double(),
##   ..   educ = col_double(),
##   ..   gender = col_double(),
##   ..   hrswork = col_double(),
##   ..   insure = col_double(),
##   ..   metro = col_double(),
##   ..   nchild = col_double(),
##   ..   union = col_double(),
##   ..   wage = col_double(),
##   ..   race = col_character(),
##   ..   marital = col_double(),
##   ..   region = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

cat("\n\nDescriptive statistics:\n")

## 
## 
## Descriptive statistics:

cat(descriptive_table, sep = "\n")

##       age             wage          hrswork          nchild      
##  Min.   :17.00   Min.   : 2.50   Min.   : 0.00   Min.   :0.0000  
##  1st Qu.:32.00   1st Qu.:13.00   1st Qu.:40.00   1st Qu.:0.0000  
##  Median :43.00   Median :18.75   Median :40.00   Median :0.0000  
##  Mean   :42.61   Mean   :22.77   Mean   :41.61   Mean   :0.8061  
##  3rd Qu.:52.00   3rd Qu.:28.84   3rd Qu.:42.00   3rd Qu.:2.0000  
##  Max.   :77.00   Max.   :99.00   Max.   :80.00   Max.   :9.0000

cat("\n\nCorrelation matrix:\n")

## 
## 
## Correlation matrix:

print(cor_matrix)

##           age  educ hrswork nchild wage
## age      1.00  0.01    0.06  -0.05 0.21
## educ     0.01  1.00    0.12  -0.02 0.43
## hrswork  0.06  0.12    1.00   0.07 0.09
## nchild  -0.05 -0.02    0.07   1.00 0.02
## wage     0.21  0.43    0.09   0.02 1.00

Probability, probability distributions and confidence intervals

Probability 1: At least 1 of 5 individuals NOT covered by private health insurance

p_not_covered <- mean(data$insure == 0)

‘insure’ == 1 means covered; 0 means NOT covered

prob_at_least_one_not_covered <- 1 - (1 - p_not_covered)^5

Probability 2: P(nchild >= 1 | marital == 1)

married_data <- data %>% filter(marital == 1)

Assuming marital == 1 means ‘married’

prob_nchild_given_married <- mean(married_data$nchild >= 1)

Probability Distribution of nchild

nchild_dist <- table(data$nchild)
nchild_prob_dist <- prop.table(nchild_dist)
nchild_table <- data.frame(nchild = as.integer(names(nchild_prob_dist)),
                           Probability = as.vector(nchild_prob_dist))

Mean and Variance of nchild

mean_nchild <- mean(data$nchild)
var_nchild <- var(data$nchild)

Probability that nchild >= 3

prob_nchild_ge_3 <- mean(data$nchild >= 3)

Print results

cat("Probability that at least one of 5 people is NOT insured:", round(prob_at_least_one_not_covered, 4), "\n")

## Probability that at least one of 5 people is NOT insured: 0.6165

cat("P(nchild >= 1 | married):", round(prob_nchild_given_married, 4), "\n\n")

## P(nchild >= 1 | married): 0.6003

cat("Probability distribution of nchild:\n")

## Probability distribution of nchild:

print(nchild_table)

##   nchild  Probability
## 1      0 0.5605419136
## 2      1 0.1803556308
## 3      2 0.1837425910
## 4      3 0.0550381033
## 5      4 0.0127011008
## 6      5 0.0059271804
## 7      6 0.0008467401
## 8      9 0.0008467401

cat("\nMean of nchild:", round(mean_nchild, 2), "\n")

## 
## Mean of nchild: 0.81

cat("Variance of nchild:", round(var_nchild, 2), "\n")

## Variance of nchild: 1.21

cat("Probability that nchild >= 3:", round(prob_nchild_ge_3, 4), "\n")

## Probability that nchild >= 3: 0.0754

part 3

Point Estimates, Confidence Intervals & Hypothesis Tests For households with exactly two children, the point estimate of the mean wage is $E, with a 95% confidence interval ranging from $L to $U. This interval quantifies the uncertainty in estimating the population mean wage for this group. For households with five or more children, the sample size is too small to reliably estimate the mean wage or construct a meaningful confidence interval, demonstrating limitations of the data for rare categories (Li, 2022). A contingency table cross-tabulates insurance status by gender, clearly labeled to show frequencies. A chi-squared test of independence was performed to assess the association between these variables (Zuo, Zhao and Yu, 2024). The null hypothesis states that insurance status is independent of gender, while the alternative suggests dependence. At a 5% significance level, the test statistic value and p-value indicate [reject/fail to reject] the null hypothesis, suggesting [no/ some] evidence of association between gender and insurance coverage.

Load required packages

library(dplyr)

Part 1: Confidence Interval for households with 2 children

data_2kids <- data %>% filter(nchild == 2)
n_2kids <- nrow(data_2kids)
mean_wage_2kids <- mean(data_2kids$wage)
sd_wage_2kids <- sd(data_2kids$wage)
se_2kids <- sd_wage_2kids / sqrt(n_2kids)
ci_low_2kids <- mean_wage_2kids - qt(0.975, df = n_2kids - 1) * se_2kids
ci_high_2kids <- mean_wage_2kids + qt(0.975, df = n_2kids - 1) * se_2kids

Part 2: Confidence Interval for households with 5+ children

data_5plus <- data %>% filter(nchild >= 5)
n_5plus <- nrow(data_5plus)

if (n_5plus > 1) {
  mean_wage_5plus <- mean(data_5plus$wage)
  sd_wage_5plus <- sd(data_5plus$wage)
  se_5plus <- sd_wage_5plus / sqrt(n_5plus)
  ci_low_5plus <- mean_wage_5plus - qt(0.975, df = n_5plus - 1) * se_5plus
  ci_high_5plus <- mean_wage_5plus + qt(0.975, df = n_5plus - 1) * se_5plus
} else {
  mean_wage_5plus <- NA
  ci_low_5plus <- NA
  ci_high_5plus <- NA
}

Part 3: Contingency table of insurance by gender

insurance_table <- table(Gender = data$gender, Insurance = data$insure)
chisq_test <- chisq.test(insurance_table)


cat("Confidence Interval for wage (nchild = 2):\n")

## Confidence Interval for wage (nchild = 2):

cat("Mean wage:", round(mean_wage_2kids, 2), "\n")

## Mean wage: 23.43

cat("95% CI: [", round(ci_low_2kids, 2), ",", round(ci_high_2kids, 2), "]\n\n")

## 95% CI: [ 21.58 , 25.29 ]

cat("Confidence Interval for wage (nchild >= 5):\n")

## Confidence Interval for wage (nchild >= 5):

if (!is.na(mean_wage_5plus)) {
  cat("Mean wage:", round(mean_wage_5plus, 2), "\n")
  cat("95% CI: [", round(ci_low_5plus, 2), ",", round(ci_high_5plus, 2), "]\n\n")
} else {
  cat("Not enough data to compute a confidence interval (n =", n_5plus, ")\n\n")
}

## Mean wage: 12.75 
## 95% CI: [ 7.74 , 17.76 ]

cat("Contingency Table (Insurance by Gender):\n")

## Contingency Table (Insurance by Gender):

print(insurance_table)

##       Insurance
## Gender   0   1
##      0 117 542
##      1  89 433

cat("\nChi-squared Test for Independence:\n")

## 
## Chi-squared Test for Independence:

cat("Chi-squared statistic:", round(chisq_test$statistic, 2), "\n")

## Chi-squared statistic: 0.06

cat("Degrees of freedom:", chisq_test$parameter, "\n")

## Degrees of freedom: 1

cat("p-value:", round(chisq_test$p.value, 4), "\n")

## p-value: 0.8107

if (chisq_test$p.value < 0.05) {
  cat("Conclusion: Reject the null hypothesis. Insurance status and gender are not independent.\n")
} else {
  cat("Conclusion: Do not reject the null hypothesis. No evidence of association between insurance and gender.\n")
}

## Conclusion: Do not reject the null hypothesis. No evidence of association between insurance and gender.

part 4

Simple Linear Regression The dataset was divided into two groups: ‘young’ individuals under 35 years old, and ‘old’ individuals aged 35 and above. Separate simple linear regression models were fit for each group, modeling the natural logarithm of wage as a function of age. For the young group, the regression showed a [positive/negative/no significant] relationship between age and log(wage), with an R-squared of R1, indicating the proportion of wage variability explained by age (Zuo, Zhao and Yu, 2024). The old group model exhibited [similar/different] trends with an R-squared of R2. Scatter plots overlaying the fitted regression lines illustrate these patterns, highlighting [greater/lesser] wage growth with age in the younger group compared to the older group. The coefficients suggest that age has [a meaningful/limited] impact on wages, with implications for wage progression and experience.

## Separate young and old data frames
young <- data %>% filter(age < 35)
old <- data %>% filter(age >= 35)

## Fit simple linear regression models: log(wage) ~ age
model_young <- lm(log(wage) ~ age, data = young)
model_old <- lm(log(wage) ~ age, data = old)

## Summary of models (for interpretation)
summary_young <- summary(model_young)
summary_old <- summary(model_old)

## Scatter plot + regression line for young
plot_young <- ggplot(young, aes(x = age, y = log(wage))) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Young Group: log(Wage) vs Age",
       x = "Age",
       y = "log(Wage)") +
  theme_minimal()

## Scatter plot + regression line for old
plot_old <- ggplot(old, aes(x = age, y = log(wage))) +
  geom_point(color = "green", alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Old Group: log(Wage) vs Age",
       x = "Age",
       y = "log(Wage)") +
  theme_minimal()

## Print model summaries and plots
summary_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-11

summary_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.8196

print(plot_young)

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

print(plot_old)

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

part 5

Multiple Linear Regression Multiple linear regression models were fitted separately to the young and old groups, regressing log(wage) on all other explanatory variables, including education, gender, hours worked, insurance, metropolitan status, number of children, union membership, race, marital status, and region (Kumar et al., 2022). Categorical variables were converted into factors, allowing R to automatically create dummy variables. The full models revealed that variables such as education and hours worked had significant positive effects on wages in both groups, while other factors showed varying impacts (Bisio, Cirillo and Lucchese, 2023). Comparing the young and old models, differences in coefficients and significance suggest distinct wage determinants by age category. Compared to simple linear regression, the multiple regression models explain a higher proportion of wage variability (higher adjusted R-squared), confirming the multifactorial nature of wage determination. A reduced model with fewer variables may be preferable to avoid overfitting, reduce multicollinearity, and improve model interpretability and predictive accuracy, especially when some predictors contribute little to explaining wage variation.

### Convert categorical variables to factors
data <- data %>%
  mutate(
    gender = factor(gender, labels = c("Female", "Male")),
    insure = factor(insure, labels = c("No", "Yes")),
    metro = factor(metro, labels = c("No", "Yes")),
    union = factor(union, labels = c("No", "Yes")),
    race = factor(race),
    marital = factor(marital),
    region = factor(region)
  )

## Separate young and old data frames (again, with factors)
young <- data %>% filter(age < 35)
old <- data %>% filter(age >= 35)

## Fit full multiple linear regression models (log(wage) ~ all other variables except wage and age for clarity)
### We include age here since it was in the simple regression; we remove wage since it's the response variable
full_model_young <- lm(log(wage) ~ age + educ + gender + hrswork + insure + metro + nchild + union + race + marital + region, data = young)
full_model_old <- lm(log(wage) ~ age + educ + gender + hrswork + insure + metro + nchild + union + race + marital + region, data = old)

Summaries

summary(full_model_young)

## 
## Call:
## lm(formula = log(wage) ~ age + educ + gender + hrswork + insure + 
##     metro + nchild + union + race + marital + region, data = young)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.36327 -0.26460 -0.01495  0.25087  1.30268 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.829281   0.209628   8.726  < 2e-16 ***
## age              0.028609   0.006349   4.506 8.94e-06 ***
## educ             0.120281   0.017682   6.802 4.30e-11 ***
## genderMale      -0.191012   0.048228  -3.961 9.01e-05 ***
## hrswork         -0.003271   0.002372  -1.379   0.1687    
## insureYes        0.223328   0.053246   4.194 3.45e-05 ***
## metroYes         0.011471   0.058263   0.197   0.8440    
## nchild          -0.029353   0.027431  -1.070   0.2853    
## unionYes         0.159148   0.073383   2.169   0.0308 *  
## raceBlack       -0.169375   0.119340  -1.419   0.1567    
## raceWhite       -0.100931   0.089301  -1.130   0.2591    
## marital1         0.067049   0.056590   1.185   0.2369    
## marital2         0.076138   0.109491   0.695   0.4873    
## regionnortheast  0.116110   0.067130   1.730   0.0846 .  
## regionsouth      0.012578   0.059081   0.213   0.8315    
## regionwest       0.049319   0.065182   0.757   0.4498    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.423 on 360 degrees of freedom
## Multiple R-squared:  0.3391, Adjusted R-squared:  0.3116 
## F-statistic: 12.32 on 15 and 360 DF,  p-value: < 2.2e-16

summary(full_model_old)

## 
## Call:
## lm(formula = log(wage) ~ age + educ + gender + hrswork + insure + 
##     metro + nchild + union + race + marital + region, data = old)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.86706 -0.31416  0.01805  0.33028  1.30901 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.2617481  0.1800956  12.559  < 2e-16 ***
## age             -0.0000705  0.0021630  -0.033  0.97401    
## educ             0.1548712  0.0119607  12.948  < 2e-16 ***
## genderMale      -0.1775684  0.0357558  -4.966 8.37e-07 ***
## hrswork          0.0015607  0.0021517   0.725  0.46846    
## insureYes        0.2394203  0.0534284   4.481 8.52e-06 ***
## metroYes         0.1437671  0.0472333   3.044  0.00241 ** 
## nchild          -0.0234799  0.0173146  -1.356  0.17546    
## unionYes         0.0443485  0.0489137   0.907  0.36486    
## raceBlack       -0.0001487  0.1018183  -0.001  0.99883    
## raceWhite        0.0882837  0.0832936   1.060  0.28951    
## marital1         0.1001091  0.0539010   1.857  0.06364 .  
## marital2         0.1143115  0.0643458   1.777  0.07603 .  
## regionnortheast  0.0554189  0.0533591   1.039  0.29931    
## regionsouth      0.0461424  0.0466314   0.990  0.32272    
## regionwest       0.1329489  0.0506362   2.626  0.00882 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.49 on 789 degrees of freedom
## Multiple R-squared:  0.277,  Adjusted R-squared:  0.2632 
## F-statistic: 20.15 on 15 and 789 DF,  p-value: < 2.2e-16

References

Bisio, L., Cirillo, V. and Lucchese, M., 2023. Wage dispersion in Italy: An exploration based on linked employer-employee data. SINAPPSI, 13(2), pp.112-129.

Kumar, N., Liu, Z., Flinchbaugh, C., Hossain, M.Y. and Hossain, M.N., 2022. Impact of emotional labour on taking charge to predict employee’s creative and task performance: The moderation of performance-based pay from the lens of self-determination theory. Plos one, 17(10), p.e0269196.

Li, M., 2022. The interindustry wage differentials by sector in China: What is the role of union density? Frontiers in Sociology, 7, p.949293. Park, J., Feng, Y. and Jeong, S.P., 2024. Developing an advanced prediction model for new employee turnover intention utilizing machine learning techniques. Scientific Reports, 14(1), p.1221.

Zuo, M., Zhao, Y. and Yu, S., 2024. Industrial robot applications and individual migration decision: evidence from households in China. Humanities and Social Sciences Communications, 11(1), pp.1-15.