Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-06-16
Modelo de Regressão Linear Múltipla - Multicolinearidade
Multicolinearidade
Na regressão linear múltipla, espera-se que as variáveis explicativas sejam, em alguma medida, independentes entre si. Quando existe uma alta correlação entre duas ou mais variáveis preditoras, estamos diante da multicolinearidade, o que pode levar a:
- Estimativas instáveis dos coeficientes;
- Aumento dos erros padrão;
- Dificuldade de interpretar o impacto de cada variável.
Correlação Amostral
Podemos identificar a multicolinearidade através da análise da matriz de correlação entre as variáveis explicativas. Uma correlação amostral elevada (valores próximos de 1 ou -1) entre duas variáveis indica que elas carregam informações semelhantes.
Diversos autores discutem os limiares de correlação que podem indicar a presença de multicolinearidade. É importante destacar que esses valores não são regras absolutas, mas orientações práticas baseadas na experiência aplicada com modelos de regressão. O ideal é que a avaliação da multicolinearidade considere tanto a análise da matriz de correlação quanto indicadores adicionais, como o Fator de Inflação da Variância (VIF).
Referências que consideram que há multicolinearidade entre as covariáveis \(X_{j}\) e \(X_{k}\) quando a correlação amostral entre elas for, em valor absoluto, maior que 0,7:
Gujarati, D. N. (2006). Econometria Básica. É uma referência clássica em cursos de econometria.
Wooldridge, J. M. (2006). Introdução à Econometria: Uma abordagem moderna.
Referências que consideram que há multicolinearidade entre as covariáveis \(x_{j}\) e \(x_{k}\) quando a correlação amostral entre elas for, em valor absoluto, maior que 0,8:
Kutner, Nachtsheim e Neter (2004). Applied Linear Regression Models.
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2012). Introduction to Linear Regression Analysis.
Fator de Inflação da Variância (VIF - Variance Inflation Factor)
O VIF (Variance Inflation Factor) mede o quanto a variância estimada de um coeficiente de regressão está sendo inflacionada devido à colinearidade com as demais variáveis explicativas do modelo.
Para cada variável explicativa \(X_j\), o VIF é calculado a partir do coeficiente de determinação \(R_j^2\), sendo
\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]
onde \(R_j^2\) representa o coeficiente de determinação da regressão da variável \(X_j\) sobre todas as outras variáveis independentes.
Exemplo com três variáveis explicativas
Suponha um modelo com três variáveis explicativas: \(X_1\), \(X_2\) e \(X_3\). Para calcular o VIF de cada uma delas:
- O \(R_1^2\) é obtido ajustando o seguinte modelo:
\[X_1 \stackrel{indep}{\sim} N(\gamma_0 + \gamma_1 X_2 + \gamma_2 X_3, V_1)\]
- O \(R_2^2\) é obtido a partir do modelo:
\[X_2 \stackrel{indep}{\sim} N(\gamma_0^{*} + \gamma_1^{*} X_1 + \gamma_2^{*} X_3, V_2)\]
- O \(R_3^2\) vem da regressão:
\[X_3 \stackrel{indep}{\sim} N(\gamma_0^{**} + \gamma_1^{**} X_1 + \gamma_2^{**} X_2, V_3).\]
Após calcular os \(R_j^2\), basta aplicar a fórmula para obter o VIF de cada variável.
Interpretação dos valores de VIF
- VIF ≈ 1: Não há multicolinearidade envolvendo a variável \(X_j\).
- 1 < VIF < 5: Indica uma correlação moderada, geralmente aceitável.
- VIF > 5: Indica forte multicolinearidade, sendo um sinal de alerta.
Importante: Esses valores de corte são orientações práticas, não limites fixos.
1 Simulando um conjunto de dados
- Simule um conjunto de dados do modelo de regressão linear considerando n=1000 (tamanho amostral), \(\beta_0=2\), \(\beta_1=3\), \(\beta_2=4\), \(\beta_3=5\), \(\sigma^2=1\) e que \(X_{i1} \stackrel{iid}{\sim} N(0,1)\), \(X_{i2} = X_{i1} + \nu_i\), \(\nu_i \stackrel{iid}{\sim} N(0; 10^{-5})\) e \(X_{i3} \stackrel{iid}{\sim} N(0,1)\).
# Configurações iniciais
n = 1000 # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = 1) # Covariável 1
x2 = x1 + rnorm(n, mean = 0, sd = sqrt(10^(-5))) # Covariável 2
x3 = rnorm(n, mean = 0, sd = 1) # Covariável 3
x = cbind(1, x1, x2, x3)
#definindo total de covariaveis
p = ncol(x)-1
# Definindo os coeficientes reais
beta_0 = 2 # Intercepto
beta_1 = 3 # Coeficiente para x1
beta_2 = 4 # Coeficiente para x2
beta_3 = 5 # Coeficiente para x3
beta = c(beta_0, beta_1, beta_2, beta_3)
rm(beta_0, beta_1, beta_2, beta_3)
# Gerando a variável resposta com erro aleatório
sigma2 = 1
e = rnorm(n, mean = 0, sd = sqrt(sigma2)) # Erro aleatório
y = x %*% beta + e
# Criando um data frame com os dados
dados = data.frame(y, x1=x1, x2=x2, x3=x3)2 Analisando os dados
- Analise o conjunto de dados simulado.
#analisando os dados gerados
par(mfrow=c(1,2))
hist(y, main="", freq=F, bty="n", ylab="densidade")
hist(e, main="", freq=F, bty="n", ylab="densidade")
par(mfrow=c(1,3))
plot(x1, y, bty="n")
plot(x2, y, bty="n")
plot(x3, y, bty="n")
par(mfrow=c(1,3))
plot(x1, x2, bty="n")
plot(x1, x3, bty="n")
plot(x2, x3, bty="n")
round(cor(dados),3)## y x1 x2 x3
## y 1.000 0.787 0.787 0.618
## x1 0.787 1.000 1.000 0.015
## x2 0.787 1.000 1.000 0.015
## x3 0.618 0.015 0.015 1.0003 Ajustando o MRL aos dados
- Ajuste um modelo de regressão linear múltipla e obtenha as estimativas pontual e intervalar para \(\beta_j\), \(j=0, \ldots, p\), considerando um nível de confiança de 90%.
#ajustando o modelo
ajuste = lm(y ~ x1 + x2 + x3, data = dados)
resultado = summary(ajuste)
#vendo os resultados do ajuste
#ajuste
resultado##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1892 -0.6461 -0.0422 0.7032 3.4832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.97868 0.03115 63.527 <2e-16 ***
## x1 -10.12726 9.86042 -1.027 0.3046
## x2 17.14320 9.85918 1.739 0.0824 .
## x3 5.03756 0.03070 164.089 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9832 on 996 degrees of freedom
## Multiple R-squared: 0.9864, Adjusted R-squared: 0.9864
## F-statistic: 2.41e+04 on 3 and 996 DF, p-value: < 2.2e-16cbind(ajuste$coefficients, confint(ajuste, level=0.9), beta)## 5 % 95 % beta
## (Intercept) 1.978684 1.9274042 2.029964 2
## x1 -10.127256 -26.3613105 6.106798 3
## x2 17.143202 0.9111889 33.375216 4
## x3 5.037555 4.9870113 5.088100 5#Notem o que acontece quando retiramos uma das covariaveis problematicas do ajuste
summary(lm(y ~ x1 + x3, data = dados) )##
## Call:
## lm(formula = y ~ x1 + x3, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2584 -0.6530 -0.0313 0.6816 3.5589
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.98073 0.03116 63.57 <2e-16 ***
## x1 7.01800 0.03338 210.25 <2e-16 ***
## x3 5.03639 0.03072 163.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9842 on 997 degrees of freedom
## Multiple R-squared: 0.9864, Adjusted R-squared: 0.9863
## F-statistic: 3.608e+04 on 2 and 997 DF, p-value: < 2.2e-16summary(lm(y ~ x2 + x3, data = dados) )##
## Call:
## lm(formula = y ~ x2 + x3, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2302 -0.6478 -0.0301 0.6969 3.5280
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.97990 0.03113 63.61 <2e-16 ***
## x2 7.01728 0.03334 210.46 <2e-16 ***
## x3 5.03687 0.03069 164.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9833 on 997 degrees of freedom
## Multiple R-squared: 0.9864, Adjusted R-squared: 0.9864
## F-statistic: 3.615e+04 on 2 and 997 DF, p-value: < 2.2e-16summary(lm(y ~ x1 + x2, data = dados) )##
## Call:
## lm(formula = y ~ x1 + x2, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.4703 -3.3235 -0.0517 3.3501 19.2324
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1736 0.1647 13.197 <2e-16 ***
## x1 25.2758 52.1689 0.485 0.628
## x2 -18.1712 52.1624 -0.348 0.728
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.203 on 997 degrees of freedom
## Multiple R-squared: 0.6191, Adjusted R-squared: 0.6183
## F-statistic: 810.2 on 2 and 997 DF, p-value: < 2.2e-164 Identificando multicolinearidade
- Identifique multicolinearidade pela correlação amostral.
#identificando multicolinearidade pela correlacao
rho_x = cor(dados[,-1])
rho_x## x1 x2 x3
## x1 1.0000000 0.99999428 0.01544570
## x2 0.9999943 1.00000000 0.01537179
## x3 0.0154457 0.01537179 1.00000000aux = rho_x
diag(aux) = NA
which(aux>0.7, arr.ind = T)## row col
## x2 2 1
## x1 1 2- Identifique multicolinearidade pelo fator de inflação da variância (VIF).
# Calculando o VIF para identificar multicolinearidade
# Função para calcular VIF manualmente
calcular_vif <- function(data) {
# Inicializando um vetor para armazenar os VIFs
vifs <- numeric(ncol(data))
for (j in seq_along(data)) {
# Variável dependente: a i-ésima coluna
var_resposta <- data[, j]
# Variáveis independentes: todas menos a i-ésima
covariaveis <- data[, -j, drop = FALSE]
# Ajustar modelo de regressão
modelo <- lm(var_resposta ~ ., data = as.data.frame(covariaveis))
# Coeficiente de determinação R2
r2 <- summary(modelo)$r.squared
# Calcular o VIF
vifs[j] <- 1 / (1 - r2)
}
# Nomear os VIFs de acordo com as colunas originais
names(vifs) <- colnames(data)
return(vifs)
}
# Aplicando a função aos preditores (x1, x2, x3)
vifs = calcular_vif(dados[, c("x1", "x2", "x3")])
round(vifs,3)## x1 x2 x3
## 87461.395 87461.196 1.001library(car) # Pacote necessário para calcular o VIF
all.equal(vifs, vif(ajuste))## [1] TRUEwhich(vifs>5)## x1 x2
## 1 2- Decida qual covariável será removida do ajuste.
# Decidindo qual variável remover com base no VIF e significância
# Identificar a variável com maior VIF
vif_maior <- names(which.max(vifs))
# Ajustar o modelo sem a variável de maior VIF
modelo_reduzido <- update(ajuste, paste(".~.-", vif_maior))
# Resumo do modelo reduzido
summary(modelo_reduzido)##
## Call:
## lm(formula = y ~ x2 + x3, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2302 -0.6478 -0.0301 0.6969 3.5280
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.97990 0.03113 63.61 <2e-16 ***
## x2 7.01728 0.03334 210.46 <2e-16 ***
## x3 5.03687 0.03069 164.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9833 on 997 degrees of freedom
## Multiple R-squared: 0.9864, Adjusted R-squared: 0.9864
## F-statistic: 3.615e+04 on 2 and 997 DF, p-value: < 2.2e-16# Recalcular o VIF no modelo reduzido
vif_reduzido <- vif(modelo_reduzido)
vif_reduzido## x2 x3
## 1.000236 1.0002365 Seleção de Modelos através do R2 ajustado
\[R^2 = 1 - \frac{SSE}{SST}\] \[R^2_a = 1 - \frac{(n-1) SSE}{(n-p-1) SST}\] \[\underbrace{\sum_{i=1}^n{(Y_i - \bar{Y})^2}}_{SST} = \underbrace{\sum_{i=1}^n{(\hat{Y}_i-\bar{Y})^2}}_{SSR} + \underbrace{\sum_{i=1}^n{(Y_i - \hat{Y}_i)^2}}_{SSE}\]
- Ajuste os modelos abaixo e para cada ajuste calcule o R2 ajustado.
\[M_1: \quad Y_i = \beta_0 + X_{i1} \beta_1 + e_i\] \[M_2: \quad Y_i = \beta_0 + X_{i3} \beta_3 + e_i\] \[M_3: \quad Y_i = \beta_0 + X_{i1} \beta_1 + X_{i3} \beta_3 + e_i\] \[M_4: \quad Y_i = \beta_0 + X_{i1} \beta_1 + X_{i2} \beta_2 + X_{i3} \beta_3 + e_i\]
#ajustando os modelos
ajuste1 = lm(y~x1, data=dados)
ajuste2 = lm(y~x3, data=dados)
ajuste3 = lm(y~x1+x3, data=dados)
ajuste4 = lm(y~x1+x2+x3, data=dados)
#R2a - coeficiente de determinação ajustado
R2a = round(c(summary(ajuste1)$adj.r.squared, summary(ajuste2)$adj.r.squared, summary(ajuste3)$adj.r.squared, summary(ajuste4)$adj.r.squared),7)
R2a## [1] 0.6186706 0.3814889 0.9863440 0.9863717which.max(R2a)## [1] 4#calculando manualmente R2 e R2 ajustado
SSE1 = sum((y-ajuste1$fitted.values)^2)
SSE2 = sum((y-ajuste2$fitted.values)^2)
SSE3 = sum((y-ajuste3$fitted.values)^2)
SSE4 = sum((y-ajuste4$fitted.values)^2)
SSR1 = sum((ajuste1$fitted.values - mean(y))^2)
SSR2 = sum((ajuste2$fitted.values - mean(y))^2)
SSR3 = sum((ajuste3$fitted.values - mean(y))^2)
SSR4 = sum((ajuste4$fitted.values - mean(y))^2)
SST1 = SSR1 + SSE1
SST2 = SSR2 + SSE2
SST3 = SSR3 + SSE3
SST4 = SSR4 + SSE4
#nbetasj = pj + 1
nbetas1 = length(ajuste1$coefficients)
nbetas2 = length(ajuste2$coefficients)
nbetas3 = length(ajuste3$coefficients)
nbetas4 = length(ajuste4$coefficients)
R21 = 1-SSE1/SST1 #R2 do ajuste1
R2a1 = 1-(n-1)*SSE1/((n-nbetas1)*SST1) #R2 ajustado do ajuste 1
R22 = 1-SSE2/SST2 #R2 do ajuste1
R2a2 = 1-(n-1)*SSE2/((n-nbetas2)*SST2) #R2 ajustado do ajuste 2
R23 = 1-SSE3/SST3 #R2 do ajuste1
R2a3 = 1-(n-1)*SSE3/((n-nbetas3)*SST3) #R2 ajustado do ajuste 3
R24 = 1-SSE4/SST4 #R2 do ajuste1
R2a4 = 1-(n-1)*SSE4/((n-nbetas4)*SST4) #R2 ajustado do ajuste 4
#R2 - coeficiente de determinação
rbind(round(c(summary(ajuste1)$r.squared, summary(ajuste2)$r.squared, summary(ajuste3)$r.squared, summary(ajuste4)$r.squared),3),
round(c(R21, R22, R23, R24),3))## [,1] [,2] [,3] [,4]
## [1,] 0.619 0.382 0.986 0.986
## [2,] 0.619 0.382 0.986 0.986#R2a - coeficiente de determinação ajustado
rbind(round(c(summary(ajuste1)$adj.r.squared, summary(ajuste2)$adj.r.squared, summary(ajuste3)$adj.r.squared, summary(ajuste4)$adj.r.squared),3),
round(c(R2a1, R2a2, R2a3, R2a4),3))## [,1] [,2] [,3] [,4]
## [1,] 0.619 0.381 0.986 0.986
## [2,] 0.619 0.381 0.986 0.986Curiosidade
O modelo utilizado para a simulação dos dados foi:
\[Y_i = 2 + 3 X_{i1} + 4 X_{i2} + 5 X_{i3} + e_i, \quad e_i \stackrel{iid}{\sim} N(0, 1),\] em que as covariáveis \(X_{i1}\) e \(X_{i2}\) foram geradas de forma altamente correlacionada, sendo: \[X_{i2} = X_{i1} + \nu_i, \quad \nu_i \stackrel{iid}{\sim} N(0, 10^{-5}).\]
Ou seja, \(X_{i2}\) é praticamente uma cópia de \(X_{i1}\), com uma pequena perturbação aleatória.
Podemos reescrever o modelo original em função apenas de \(X_{i1}\) e \(X_{i3}\) da seguinte forma:
\[ \begin{aligned} Y_i &= 2 + 3 X_{i1} + 4 (X_{i1} + \nu_i) + 5 X_{i3} + e_i \\ &= 2 + 7 X_{i1} + 5 X_{i3} + (4 \nu_i + e_i) \\ &= 2 + 7 X_{i1} + 5 X_{i3} + \varepsilon_i, \end{aligned} \]
sendo:
\[ \varepsilon_i = 4 \nu_i + e_i. \]
Note que \(\varepsilon_i\) também segue uma distribuição normal, pois é uma combinação linear de normais independentes. Mais especificamente:
\[ \varepsilon_i \sim N\left(0,\ 4^2 \times 10^{-5} + 1^2\right) = N\left(0,\ 1{,}00016\right). \]
Essa transformação ilustra como a presença de colinearidade entre \(X_1\) e \(X_2\) pode, na prática, tornar um dos coeficientes “não identificável” ao analisarmos os efeitos conjuntos, já que os dois preditores estão praticamente fornecendo a mesma informação.
Além disso, esse exemplo evidencia como pequenos erros de medição ou variações quase imperceptíveis podem gerar problemas de multicolinearidade severa na regressão linear.
c(ajuste3$coefficients, summary(ajuste3)$sigma^2)## (Intercept) x1 x3
## 1.9807293 7.0180033 5.0363902 0.96871556 Seleção de Modelos através de AIC e BIC
Seja \(L\) a função de verossimilhança do modelo de regressão linear \(Y_i \sim N(\boldsymbol{X}_i\boldsymbol{\beta}, \sigma^2)\), \(i=1, \ldots, n\), aplicada no ponto de máxima verossimilhança. Então \[l = \ln(L) = -\frac{n}{2}\ln(2\pi \hat{\sigma}^2)-\frac{1}{2\hat{\sigma}^2}\sum_{i=1}^n{(Y_i - \boldsymbol{X}_i^T\boldsymbol{\hat{\beta}})^2},\] sendo \(\hat{\sigma}^2 = \frac{\sum_{i=1}^n{(Y_i - \boldsymbol{X}_i^T\boldsymbol{\hat{\beta}})^2}}{n-p-1} = \frac{SSE}{n-p-1}\) e \(\boldsymbol{\hat{\beta}} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{Y}\). Essa função pode ser reescrita em termos de SSE da seguinte forma \[l = -\frac{n}{2} \ln(2\pi) - \frac{n}{2} \ln\left(\frac{SSE}{n-p-1}\right)-\frac{n-p-1}{2 \; SSE} SSE \Rightarrow -2l = n \ln(2\pi) + n \ln\left(\frac{SSE}{n-p-1}\right) + (n-p-1).\]
Sejam \[AIC = -2 \ln(L) + 2 K,\] \[BIC = -2 \ln(L) + \ln(n) K,\] sendo \(K\) o total de parâmetros do modelo (dimensão de \(\boldsymbol{\beta}\) + 1 por conta da variância). Note que quando aplicarmos a função de log-verossimilhança para diferentes ajustes de MRLM, os termos \(n \ln(2\pi) + n\) serão os mesmos para os diferentes ajustes. Logo, podemos usar uma versão simplificada destas medidas calculando \[AIC^* = n \ln(SSE/n) + 2 K,\] \[BIC^* = n \ln(SSE/n) + \ln(n) K,\]
- Para cada ajuste realizado anteriormente, calcule AIC, BIC, AIC\(^*\) e BIC\(^*\).
#AIC - calculando pela funcao de verossimilhanca e pela funcao do R
AICs = c(AIC(ajuste1), AIC(ajuste2), AIC(ajuste3), AIC(ajuste4))
l1 = -2*sum(dnorm(y, ajuste1$fitted.values, summary(ajuste1)$sigma, log=TRUE))
l2 = -2*sum(dnorm(y, ajuste2$fitted.values, summary(ajuste2)$sigma, log=TRUE))
l3 = -2*sum(dnorm(y, ajuste3$fitted.values, summary(ajuste3)$sigma, log=TRUE))
l4 = -2*sum(dnorm(y, ajuste4$fitted.values, summary(ajuste4)$sigma, log=TRUE))
rbind(c(l1 + 2 * (nbetas1+1),
l2 + 2 * (nbetas2+1),
l3 + 2 * (nbetas3+1),
l4 + 2 * (nbetas4+1)),
AICs)## [,1] [,2] [,3] [,4]
## 6139.577 6623.229 2811.093 2810.065
## AICs 6139.575 6623.227 2811.088 2810.057which.min(AICs)## [1] 4#BIC - calculando pela funcao de verossimilhanca e pela funcao do R
BICs = c(BIC(ajuste1), BIC(ajuste2), BIC(ajuste3), BIC(ajuste4))
rbind(c(l1 + log(n) * (nbetas1+1),
l2 + log(n) * (nbetas2+1),
l3 + log(n) * (nbetas3+1),
l4 + log(n) * (nbetas4+1)),
BICs)## [,1] [,2] [,3] [,4]
## 6154.301 6637.952 2830.724 2834.604
## BICs 6154.299 6637.950 2830.719 2834.596which.min(BICs)## [1] 3#Calculando a versao simplificada do AIC
AICs_simp = AICs - n *log(2*pi) - n
rbind(c(n*log(SSE1) - n*log(n)+2*(nbetas1+1),
n*log(SSE2) - n*log(n)+2*(nbetas2+1),
n*log(SSE3) - n*log(n)+2*(nbetas3+1),
n*log(SSE4) - n*log(n)+2*(nbetas4+1)),
AICs_simp)## [,1] [,2] [,3] [,4]
## 3301.698 3785.35 -26.78881 -27.8198
## AICs_simp 3301.698 3785.35 -26.78881 -27.8198which.min(AICs)## [1] 4which.min(AICs_simp)## [1] 4#Calculando a versao simplificada do BIC
BICs_simp = BICs - n *log(2*pi) - n
rbind(c(n*log(SSE1) - n*log(n)+log(n)*(nbetas1+1),
n*log(SSE2) - n*log(n)+log(n)*(nbetas2+1),
n*log(SSE3) - n*log(n)+log(n)*(nbetas3+1),
n*log(SSE4) - n*log(n)+log(n)*(nbetas4+1)),
BICs_simp)## [,1] [,2] [,3] [,4]
## 3316.422 3800.073 -7.157792 -3.281025
## BICs_simp 3316.422 3800.073 -7.157792 -3.281025which.min(BICs)## [1] 3which.min(BICs_simp)## [1] 3