Easy Problem: Large N And Half The Variables Are Unrelated
library(GlmSimulatoR)
library(ggplot2)
library(MASS)
#Creating data to work with
set.seed(1)
simdata <- simulate_inverse_gaussian(N = 100000, link = "1/mu^2",
weights = c(1, 2, 3), unrelated = 3)
#Y looks like an inverse gaussian distribution.
ggplot(simdata, aes(x=Y)) +
geom_histogram(bins = 30)
#Setting the simplest model and the most complex model.
scopeArg <- list(
lower = Y ~ 1,
upper = Y ~ X1 + X2 + X3 + Unrelated1 + Unrelated2 + Unrelated3
)
#Run search
startingModel <- glm(Y ~ 1, data = simdata, family = inverse.gaussian(link = "1/mu^2"))
glmSearch <- stepAIC(startingModel, scopeArg)
#> Start: AIC=-209841.1
#> Y ~ 1
#>
#> Df Deviance AIC
#> + X3 1 33538 -211199
#> + X2 1 33789 -210467
#> + X1 1 33952 -209991
#> <none> 34005 -209841
#> + Unrelated3 1 34005 -209839
#> + Unrelated1 1 34005 -209839
#> + Unrelated2 1 34005 -209839
#>
#> Step: AIC=-211221
#> Y ~ X3
#>
#> Df Deviance AIC
#> + X2 1 33324 -211853
#> + X1 1 33485 -211375
#> <none> 33538 -211221
#> + Unrelated3 1 33538 -211219
#> + Unrelated1 1 33538 -211219
#> + Unrelated2 1 33538 -211219
#> - X3 1 34005 -209839
#>
#> Step: AIC=-211858.7
#> Y ~ X3 + X2
#>
#> Df Deviance AIC
#> + X1 1 33270 -212018
#> <none> 33324 -211859
#> + Unrelated3 1 33324 -211857
#> + Unrelated1 1 33324 -211857
#> + Unrelated2 1 33324 -211857
#> - X2 1 33538 -211221
#> - X3 1 33789 -210469
#>
#> Step: AIC=-212018.7
#> Y ~ X3 + X2 + X1
#>
#> Df Deviance AIC
#> <none> 33270 -212019
#> + Unrelated3 1 33270 -212017
#> + Unrelated1 1 33270 -212017
#> + Unrelated2 1 33270 -212017
#> - X1 1 33324 -211859
#> - X2 1 33485 -211376
#> - X3 1 33736 -210625
summary(glmSearch)
#>
#> Call:
#> glm(formula = Y ~ X3 + X2 + X1, family = inverse.gaussian(link = "1/mu^2"),
#> data = simdata)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.68546 -0.47420 -0.08869 0.29790 2.36148
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.81196 0.20842 13.49 <2e-16 ***
#> X3 3.03192 0.08116 37.36 <2e-16 ***
#> X2 2.05731 0.08100 25.40 <2e-16 ***
#> X1 1.02594 0.08066 12.72 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for inverse.gaussian family taken to be 0.3335593)
#>
#> Null deviance: 34005 on 99999 degrees of freedom
#> Residual deviance: 33270 on 99996 degrees of freedom
#> AIC: -212019
#>
#> Number of Fisher Scoring iterations: 5
rm(simdata, scopeArg, glmSearch, startingModel)Looking at the summary, the correct model was found. Stepwise search worked perfectly!
