In [1]:
# rm(list=ls())
options(OutDec = ",")
#===================================================================
# Exemplo de gastos mensais com cartão de crédito.
#===================================================================
# Análise.
#===================================================================
library(MASS)
dataxy <- read.table("../dados/exemplo_14.txt",header=T)
print(dataxy)
MDR Acc Age Income Avgexp Ownrent Selfempl 1 0 1 38 4,52 124,98 1 0 2 0 1 33 2,42 9,85 0 0 3 0 1 34 4,50 15,00 1 0 4 0 1 31 2,54 137,87 0 0 5 0 1 32 9,79 546,50 1 0 6 0 1 23 2,50 92,00 0 0 7 0 1 28 3,96 40,83 0 0 8 0 1 29 2,37 150,79 1 0 9 0 1 37 3,80 777,82 1 0 10 0 1 28 3,20 52,58 0 0 11 0 1 31 3,95 256,66 1 0 12 0 0 42 1,98 0,00 1 0 13 0 0 30 1,73 0,00 1 0 14 0 1 29 2,45 78,87 1 0 15 0 1 35 1,91 42,62 1 0 16 0 1 41 3,20 335,43 1 0 17 0 1 40 4,00 248,72 1 0 18 7 0 30 3,00 0,00 1 0 19 0 1 40 10,00 548,03 1 1 20 3 0 46 3,40 0,00 0 0 21 0 1 35 2,35 43,34 1 0 22 1 0 25 1,88 0,00 0 0 23 0 1 34 2,00 218,52 1 0 24 1 1 36 4,00 170,64 0 0 25 0 1 43 5,14 37,58 1 0 26 0 1 30 4,51 502,20 0 0 27 0 0 22 3,84 0,00 0 1 28 0 1 22 1,50 73,18 0 0 29 0 0 34 2,50 0,00 1 0 30 0 1 40 5,50 1532,77 1 0 31 0 1 22 2,03 42,69 0 0 32 1 1 29 3,20 417,83 0 0 33 1 0 25 3,15 0,00 1 0 34 0 1 21 2,47 552,72 1 0 35 0 1 24 3,00 222,54 0 0 36 0 1 43 3,54 541,30 1 0 37 0 0 43 2,28 0,00 0 0 38 0 1 37 5,70 568,77 1 0 39 0 1 27 3,50 344,47 0 0 40 0 1 28 4,60 405,35 1 0 41 0 1 26 3,00 310,94 1 0 42 0 1 23 2,59 53,65 0 0 43 0 1 30 1,51 63,92 0 0 44 0 1 30 1,85 165,85 0 0 45 0 1 38 2,60 9,58 0 0 46 0 0 28 1,80 0,00 0 1 47 0 1 36 2,00 319,49 0 0 48 0 0 38 3,26 0,00 0 0 49 0 1 26 2,35 83,08 0 0 50 0 1 28 7,00 644,83 1 0 51 0 0 50 3,60 0,00 0 0 52 0 1 24 2,00 93,20 0 0 53 0 1 21 1,70 105,04 0 0 54 0 1 24 2,80 34,13 0 0 55 0 1 26 2,40 41,19 0 0 56 1 1 33 3,00 169,89 0 0 57 0 1 34 4,80 1898,03 0 0 58 0 1 33 3,18 810,39 0 0 59 0 0 45 1,80 0,00 0 0 60 0 1 21 1,50 32,78 0 0 61 2 1 25 3,00 95,80 0 0 62 0 1 27 2,28 27,78 0 0 63 0 1 26 2,80 215,07 0 0 64 0 1 22 2,70 79,51 0 0 65 3 0 27 4,90 0,00 1 0 66 0 0 26 2,50 0,00 0 1 67 0 1 41 6,00 306,03 0 1 68 0 1 42 3,90 104,54 0 0 69 0 0 22 5,10 0,00 0 0 70 0 1 25 3,07 642,47 0 0 71 0 1 31 2,46 308,05 1 0 72 0 1 27 2,00 186,35 0 0 73 0 1 33 3,25 56,15 0 0 74 0 1 37 2,72 129,37 0 0 75 0 1 27 2,20 93,11 0 0 76 1 0 24 4,10 0,00 0 0 77 0 1 24 3,75 292,66 0 0 78 0 1 25 2,88 98,46 0 0 79 0 1 36 3,05 258,55 0 0 80 0 1 33 2,55 101,68 0 0 81 0 0 33 4,00 0,00 0 0 82 1 1 55 2,64 65,25 1 0 83 0 1 20 1,65 108,61 0 0 84 0 1 29 2,40 49,56 0 0 85 3 0 40 3,71 0,00 0 0 86 0 1 41 7,24 235,57 1 0 87 0 0 41 4,39 0,00 1 0 88 0 0 35 3,30 0,00 1 0 89 0 0 24 2,30 0,00 0 0 90 1 0 54 4,18 0,00 0 0 91 2 0 34 2,49 0,00 0 0 92 0 0 45 2,81 0,00 1 0 93 0 1 43 2,40 68,38 0 0 94 4 0 35 1,50 0,00 0 0 95 2 0 36 8,40 0,00 0 0 96 0 1 22 1,56 0,00 0 0 97 1 1 33 6,00 474,15 1 0 98 1 1 25 3,60 234,05 0 0 99 0 1 26 5,00 451,20 1 0 100 0 1 46 5,50 251,52 1 0
In [2]:
#===================================================================
# Definindo as variáveis do modelo.
#===================================================================
avgexp <- dataxy[dataxy[,5]>0,5]
idade <- dataxy[dataxy[,5]>0,3]
renda <- dataxy[dataxy[,5]>0,4]
renda2 <- renda^2
casa <- dataxy[dataxy[,5]>0,6]
n <- length(avgexp)
In [3]:
#===================================================================
# Ajuste da regressão linear por mínimos quadrados.
#===================================================================
fity <- lm(avgexp ~ idade + casa + renda + renda2 )
print(summary(fity))
Call:
lm(formula = avgexp ~ idade + casa + renda + renda2)
Residuals:
Min 1Q Median 3Q Max
-429,03 -130,36 -51,10 53,98 1460,62
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -237,147 199,352 -1,190 0,23841
idade -3,082 5,515 -0,559 0,57814
casa 27,941 82,922 0,337 0,73721
renda 234,347 80,366 2,916 0,00482 **
renda2 -14,997 7,469 -2,008 0,04870 *
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 284,8 on 67 degrees of freedom
Multiple R-squared: 0,2436, Adjusted R-squared: 0,1984
F-statistic: 5,394 on 4 and 67 DF, p-value: 0,0007952
In [4]:
#===================================================================
# Dispersão entre renda e resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
plot(renda,fity$res,pch=15,xlab="Renda",ylab="U",xlim=c(0,12),
ylim=c(-500,2000))
In [5]:
#===================================================================
# Estimador de White e outros.
#===================================================================
X <- cbind(rep(1,n),idade,casa,renda,renda2)
XtXinv <- solve(t(X)%*%X)
S0 <- matrix(0,5,5)
DM2 <- matrix(0,5,5)
for(i in 1:n){
S0 <- S0 + (fity$res[i]^2)*(X[i,]%*%t(X[i,]))
mii <- 1 - t(X[i,])%*%XtXinv%*%X[i,]
DM2 <- DM2 + c((fity$res[i]^2)/mii)*(X[i,]%*%t(X[i,]))
}
S0 <- S0/n
DM2 <- DM2/n
varbw <- n*XtXinv%*%S0%*%XtXinv
varbdm <- n*XtXinv%*%DM2%*%XtXinv
sdbw <- sqrt(diag(varbw))
sdbdm <- sqrt(diag(varbdm))
B <- cbind(NA,mean(idade),mean(casa),mean(renda),NA)
B <- rbind(B,fity$coef)
B <- rbind(B,sqrt(diag(vcov(fity))))
B <- rbind(B,B[2,]/B[3,])
B <- rbind(B,sdbw)
B <- rbind(B,sdbw*sqrt(n/(n-5)))
B <- rbind(B,sdbdm)
row.names(B) <- c("Media Amostral","Coeficiente","Erro Padrao","Razao t",
"Erro Padrao White","D. e M. (1)","D. e M. (2)")
print(B)
(Intercept) idade casa renda renda2 Media Amostral NA 31,2777778 0,3750000 3,437083 NA Coeficiente -237,146514 -3,0818140 27,9409084 234,347027 -14,996844 Erro Padrao 199,351665 5,5147165 82,9223236 80,365950 7,469337 Razao t -1,189589 -0,5588345 0,3369528 2,915999 -2,007788 Erro Padrao White 212,990530 3,3016612 92,1877767 88,866352 6,944563 D. e M. (1) 220,794952 3,4226411 95,5657314 92,122602 7,199027 D. e M. (2) 221,088927 3,4477148 95,6721114 92,083684 7,199538
In [6]:
#===================================================================
# Teste de White.
#===================================================================
w <- fity$res^2
renda4 <- renda^4
idade2 <- idade^2
fitw <- lm(w ~ idade + casa + renda + renda2
+ idade2 + renda4
+ idade:casa + idade:renda + idade:renda2
+ casa:renda + casa:renda2 + renda:renda2)
print("===========================================================")
print(summary(fitw))
R2w <- 1-sum(fitw$res^2)/sum((w-mean(w))^2)
WS <- n*R2w
print("===========================================================")
print(c(WS > qchisq(0.95,13-1)))
[1] "==========================================================="
Call:
lm(formula = w ~ idade + casa + renda + renda2 + idade2 + renda4 +
idade:casa + idade:renda + idade:renda2 + casa:renda + casa:renda2 +
renda:renda2)
Residuals:
Min 1Q Median 3Q Max
-403431 -83907 -8912 36958 1675167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1637390,4 1290979,7 1,268 0,2097
idade 5366,2 48893,8 0,110 0,9130
casa 812036,8 991630,2 0,819 0,4161
renda -2021697,6 1053559,1 -1,919 0,0598 .
renda2 669055,3 365666,7 1,830 0,0724 .
idade2 -424,1 627,5 -0,676 0,5018
renda4 3762,7 2277,4 1,652 0,1038
idade:casa 4661,7 14424,6 0,323 0,7477
idade:renda 11499,9 15614,3 0,736 0,4643
idade:renda2 -1093,3 1568,1 -0,697 0,4884
casa:renda -510192,3 469792,6 -1,086 0,2819
casa:renda2 51835,1 61799,8 0,839 0,4050
renda:renda2 -86805,3 51162,6 -1,697 0,0950 .
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 274600 on 59 degrees of freedom
Multiple R-squared: 0,199, Adjusted R-squared: 0,0361
F-statistic: 1,222 on 12 and 59 DF, p-value: 0,2905
[1] "==========================================================="
[1] FALSE
In [7]:
#===================================================================
# Teste de Breusch-Pagan e Koenker-Basset
#===================================================================
Z <- cbind(rep(1,n),renda,renda2) # duas variaveis
fitLM <- lm(avgexp ~ renda + renda2 )
g <- (fitLM$res^2)/mean(fitLM$res^2) - 1
LM <- 0.5*t(g)%*%Z%*%solve(t(Z)%*%Z)%*%t(Z)%*%g
print("===========================================================")
print(c(LM > qchisq(0.95,2)))
print("===========================================================")
print(1-pchisq(LM,2)) # valor-p
[1] "==========================================================="
[1] TRUE
[1] "==========================================================="
[,1]
[1,] 1,165549e-09
In [8]:
#===================================================================
# Teste de Koenker-Basset.
#===================================================================
V <- mean( ((fitLM$res^2) - mean(fitLM$res^2))^2 )
u <- fitLM$res^2
ubarra <- mean(fitLM$res^2)
LM2 <- (1/V)*t(u-ubarra)%*%Z%*%solve(t(Z)%*%Z)%*%t(Z)%*%(u-ubarra)
print("===========================================================")
print(c(LM2 > qchisq(0.95,2)))
print("===========================================================")
print(1-pchisq(LM2,2)) # p-valor
[1] "==========================================================="
[1] TRUE
[1] "==========================================================="
[,1]
[1,] 0,04214613
In [9]:
#===================================================================
# A diferença entre LM e LM2 sugere falta de normalidade dos dados.
#===================================================================
In [10]:
#===================================================================
# Fim
#===================================================================