Packages

install.packages("lmtest")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
library(MASS)      
library(lmtest)  
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

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. Visualise

plot(df_firms$Total_Assets, df_firms$RD_Expenditure,
     xlab = "Total Assets (mn)", ylab = "R&D Expenditure (mn)",
     main = "Total Assets vs R&D Expenditure",
     pch = 19, col = rgb(0, 0, 1, 0.4))
abline(lm(RD_Expenditure ~ Total_Assets, data = df_firms), col = "red", lwd = 2)

#2. Diagnose

model_raw <- lm(RD_Expenditure ~ Total_Assets, data = df_firms)
summary(model_raw)
## 
## 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(mar = c(4, 4, 2, 1))
par(mfrow = c(2, 2))
plot(model_raw)

par(mfrow = c(1, 1))

shapiro.test(residuals(model_raw))
bptest(model_raw)
#3. Box-Cox Transform and refine

bc <- boxcox(model_raw, lambda = seq(-2, 2, 0.05))

lambda_opt <- bc$x[which.max(bc$y)]
cat("Optimal lambda:", lambda_opt, "\n")
## Optimal lambda: 0.1818182
df_firms$RD_transformed <- (df_firms$RD_Expenditure ^ lambda_opt - 1) / lambda_opt

model_bc <- lm(RD_transformed ~ Total_Assets, data = df_firms)
summary(model_bc)
## 
## 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_bc)

par(mfrow = c(1, 1))

shapiro.test(residuals(model_bc))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_bc)
## W = 0.98986, p-value = 0.1707
bptest(model_bc)               
## 
##  studentized Breusch-Pagan test
## 
## data:  model_bc
## BP = 0.0051706, df = 1, p-value = 0.9427