##----------------------------------------------------------------------------- ## Definições da sessão. require(lattice) require(latticeExtra) require(car) RSiteSearch("VIF regression") require(faraway) ## tem a função vif() require(fmsb) ## tem a função VIF() ## Função vif da página da Professora Dra Sueli Giolo. source("http://people.ufpr.br/~giolo/CE071/Exemplos/vif.R") ## Several packages in R provide functions to calculate VIF: vif in ## package HH, vif in package car, VIF in package fmsb, vif in package ## faraway, and vif in package VIF. The number of packages that provide ## VIF functions is surprising given that they all seem to accomplish ## the same thing. ## ## http://beckmw.wordpress.com/2013/02/05/collinearity-and-stepwise-vif-selection ## http://courses.ttu.edu/isqs5349-westfall/images/5349/multicollinearity_99.htm ## Considere o conjunto de dados turtles.txt em ## http://westfall.ba.ttu.edu/isqs5349/Rdata/turtles.txtx ##----------------------------------------------------------------------------- ## Dados sobre o preço de leitão (price) de relógios antigos (do avô) em ## função da idade do relógio (age) e do número de potenciais ## compradores (bidders). prelink <- "http://www.leg.ufpr.br/~walmes/data/business_economics_dataset" da <- read.table(paste(prelink, "/EXAMPLES/FTC.DAT", sep=""), header=FALSE) str(da) names(da) <- c("tar", "nicotine", "weight", "co") str(da) ## tar: conteúdo de alcatrão; ## nicotine: conteúdo de nicotina; ## weight: peso; ## co: monoxido de carbono; ## Os valores no data.frame são dos valores de alcatrão, nicotina e ## monoxido de carbono (mg) e peso (g) para uma amostra de 25 marcas de ## filtros testados. Deseja-se modelar o monoxido de carbono como função ## das demais variáveis. m0 <- lm(co~tar+nicotine+weight, data=da) summary(m0) pairs(da) par(mfrow=c(2,2)) plot(m0) layout(1) residualPlots(m0) m0 <- lm(co~tar+nicotine+weight, data=da[-3,]) summary(m0) pairs(da) cor(da[-3,]) vcov(m0) cov2cor(vcov(m0)) ## Retorna um único valor. ## VIF(m0) ## Um valor para cada variável. vif(m0) ## Calculando pela expressão do livro. r2 <- summary(lm(tar~nicotine+weight, data=da[-3,]))$r.squared 1/(1-r2) ## Calculando pelas correlações entre variáveis. X <- with(da[-3,], cbind(tar, nicotine, weight)) cX <- cor(X) icX <- solve(cX) vifs <- diag(icX); vifs m0 <- lm(co~tar+weight, data=da[-3,]) summary(m0) vif(m0) ## vif_i ~= 1: variável i não envolvida em multicolinearidade; ## vif_i > 5: variável i envolvida em multicolinearidade, \beta_i com ## variância alta e fracamente estimado. ##----------------------------------------------------------------------------- str(swiss) pairs(swiss) help(swiss, help_type="html") m0 <- lm(Fertility~., data=swiss) vif(m0) str(rock) pairs(rock) m0 <- lm(perm~., data=rock) vif(m0) str(longley) help(longley, help_type="html") pairs(longley) m0 <- lm(Employed~., data=longley) vif(m0) ##----------------------------------------------------------------------------- RSiteSearch("ridge regression") require(ridge) lR <- linearRidge(co~tar+nicotine+weight, data=da) lR <- linearRidge(Employed~., data=longley) summary(lR) plot(lR) require(genridge) longley.y <- longley[, "Employed"] longley.X <- data.matrix(longley[, c(2:6,1)]) lambda <- c(0, 0.005, 0.01, 0.02, 0.04, 0.08) lridge <- ridge(longley.y, longley.X, lambda=lambda) traceplot(lridge) # same, using formula interface lridge <- ridge(Employed ~ GNP + Unemployed + Armed.Forces + Population + Year + GNP.deflator, data=longley, lambda=lambda) coef(lridge) traceplot(lridge) traceplot(lridge, X="df") pairs(lridge, radius=0.5) ##-----------------------------------------------------------------------------