Libraries

library(ggplot2)
library(MASS)
library(lmtest)

Data Generation

set.seed(42)
n <- 200
firm_size <- runif(n, 10, 500)
error <- rnorm(n, mean = 0, sd = 0.5)
rd_expenditure <- exp(1.5 + 0.6 * log(firm_size) + error)

df_firms <- data.frame(
  Firm_ID        = 1:n,
  Total_Assets   = firm_size,
  RD_Expenditure = rd_expenditure
)

head(df_firms)
##   Firm_ID Total_Assets RD_Expenditure
## 1       1     458.2550      322.76939
## 2       2     469.1670      302.76313
## 3       3     150.2084       54.90529
## 4       4     416.9193      421.56611
## 5       5     324.4553      103.12089
## 6       6     264.3570      134.17397

1. Visualize: Total Assets vs R&D Expenditure

ggplot(df_firms, aes(x = Total_Assets, y = RD_Expenditure)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Relationship Between Total Assets and R&D Expenditure",
    x     = "Total Assets (Millions)",
    y     = "R&D Expenditure (Millions)"
  ) +
  theme_minimal()

2. Diagnose: OLS Regression and Residual Checks

model_ols <- lm(RD_Expenditure ~ Total_Assets, data = df_firms)
summary(model_ols)
## 
## Call:
## lm(formula = RD_Expenditure ~ Total_Assets, data = df_firms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -135.79  -42.06  -12.37   25.08  404.97 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  40.50788   11.25914   3.598 0.000405 ***
## Total_Assets  0.35091    0.03731   9.405  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 75.31 on 198 degrees of freedom
## Multiple R-squared:  0.3088, Adjusted R-squared:  0.3053 
## F-statistic: 88.46 on 1 and 198 DF,  p-value: < 2.2e-16
par(mfrow = c(2, 2))
plot(model_ols, main = "OLS Residual Diagnostics")

par(mfrow = c(1, 1))

Normality Test (Shapiro-Wilk)

shapiro.test(residuals(model_ols))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_ols)
## W = 0.87801, p-value = 1.231e-11

Homoscedasticity Test (Breusch-Pagan)

bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 13.298, df = 1, p-value = 0.0002657

Interpretation: A p-value < 0.05 in the Shapiro-Wilk test indicates non-normal residuals. A p-value < 0.05 in the Breusch-Pagan test indicates heteroscedasticity.

3. Transform: Box-Cox to Find Optimal Lambda

bc <- boxcox(model_ols, lambda = seq(-2, 2, by = 0.1))

optimal_lambda <- bc$x[which.max(bc$y)]
cat("Optimal Lambda:", optimal_lambda, "\n")
## Optimal Lambda: 0.1818182

Interpretation: Lambda ≈ 0 suggests a log transformation is most appropriate.

4. Refine: Re-run Model with Transformed Variable

if (abs(optimal_lambda) < 0.1) {
  df_firms$RD_transformed <- log(df_firms$RD_Expenditure)
  transform_label <- "log(RD_Expenditure)"
} else {
  df_firms$RD_transformed <- (df_firms$RD_Expenditure^optimal_lambda - 1) / optimal_lambda
  transform_label <- paste0("(RD_Expenditure^", round(optimal_lambda, 2), " - 1) / lambda")
}
model_refined <- lm(RD_transformed ~ Total_Assets, data = df_firms)
summary(model_refined)
## 
## Call:
## lm(formula = RD_transformed ~ Total_Assets, data = df_firms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6671 -0.9255 -0.0303  0.7700  3.2087 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.4627033  0.1839737   29.69   <2e-16 ***
## Total_Assets 0.0075331  0.0006096   12.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.231 on 198 degrees of freedom
## Multiple R-squared:  0.4354, Adjusted R-squared:  0.4325 
## F-statistic: 152.7 on 1 and 198 DF,  p-value: < 2.2e-16
par(mfrow = c(2, 2))
plot(model_refined, main = "Refined Model Residual Diagnostics")

par(mfrow = c(1, 1))

Normality Test (Shapiro-Wilk)

shapiro.test(residuals(model_refined))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_refined)
## W = 0.98986, p-value = 0.1707

Homoscedasticity Test (Breusch-Pagan)

bptest(model_refined)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_refined
## BP = 0.0051706, df = 1, p-value = 0.9427
ggplot(df_firms, aes(x = Total_Assets, y = RD_transformed)) +
  geom_point(alpha = 0.6, color = "darkorange") +
  geom_smooth(method = "lm", color = "darkred", se = TRUE) +
  labs(
    title = paste("Refined Model:", transform_label, "~ Total_Assets"),
    x     = "Total Assets (Millions)",
    y     = transform_label
  ) +
  theme_minimal()

Model Comparison Summary

comparison <- data.frame(
  Model         = c("OLS (Untransformed)", "Refined (Box-Cox)"),
  R_Squared     = c(round(summary(model_ols)$r.squared, 4),
                    round(summary(model_refined)$r.squared, 4)),
  Adj_R_Squared = c(round(summary(model_ols)$adj.r.squared, 4),
                    round(summary(model_refined)$adj.r.squared, 4)),
  Shapiro_pval  = c(round(shapiro.test(residuals(model_ols))$p.value, 4),
                    round(shapiro.test(residuals(model_refined))$p.value, 4)),
  BP_pval       = c(round(bptest(model_ols)$p.value, 4),
                    round(bptest(model_refined)$p.value, 4))
)

knitr::kable(comparison, caption = "Model Comparison: OLS vs Box-Cox Transformed")
Model Comparison: OLS vs Box-Cox Transformed
Model R_Squared Adj_R_Squared Shapiro_pval BP_pval
OLS (Untransformed) 0.3088 0.3053 0.0000 0.0003
Refined (Box-Cox) 0.4354 0.4325 0.1707 0.9427

Conclusion: The refined model shows improved normality and homoscedasticity in residuals, confirming that the Box-Cox transformation (log) better captures the true relationship between firm size and R&D expenditure.