Lista de Regressão Quantilica

Lista de Regressão Quantilica

Questão 3

• No R faça gráficos da função de distribuição acumulada e função quartil para dados simulados Poisson e Exponencial

• Com dados de simulação da distribuição Poisson

Código

set.seed(2023)
y <- rpois(2000,4)
pmid <- midecdf(y)
xmid <- midquantile(y, probs = pmid$y)
par(mfrow = c(1,2))
plot(pmid, xlab = "y", ylab = "CDF", jumps = TRUE)
points(pmid$x, pmid$y, pch = 15)
plot(xmid, xlab = "p", ylab = "Quantile", jumps = TRUE)
points(xmid$x, xmid$y, pch = 15)

• Com dados de simulação da distribuição Exponencial

Código

set.seed(2023)
y <- rexp(2000,4)
pmid <- midecdf(y)
xmid <- midquantile(y, probs = pmid$y)
par(mfrow = c(1,2))
plot(pmid, xlab = "y", ylab = "CDF", jumps = TRUE)
points(pmid$x, pmid$y, pch = 15)
plot(xmid, xlab = "p", ylab = "Quantile", jumps = TRUE)
points(xmid$x, xmid$y, pch = 15)

Questão 5

• Usando os dados “faithful” do R que contém 272 observações do tempo de espera em minutos (waiting) entre erupções e a duração em minutos das erupções (eruptions) no Parque Nacional Yellowstone, responda as questões:

• 1. Faça um gráfico da variável waiting (y) contra eruptions (x). Você consegue identificar algum padrão? Justifique?

Código

data("faithful")
plot_ly(faithful, x = faithful$eruptions, y = faithful$waiting, text = paste("Eruptions:", faithful$eruptions, "\nwaiting:", faithful$waiting)) %>% 
  layout(title = "Diagrama de dispersão entre waiting e eruptions", xaxis = list(title = "Eruptions"), yaxis = list(title = "Waiting"))

• Sim, podemos observar um padrão, é visto uma correlação alta, positiva entre as duas variáveis. E também dois grupos bem divididos, acima de 3 e abaixo de 3 com relação a duração em minutos das erupções.

• 2. Crie uma variável dummy para o tempo de duração das erupções igual ou maior do que 3 minutos. Ajuste um modelo de regressão quantílica para explicar o tempo de espera entre as erupções a partir da dummy criada. Considere os quantis 0.10, 0.25, 0.50, 0.75, 0.90. O que você pode concluir? Existe comportamento diferente entre os quantis? Justifique.

Código

banco <- faithful %>% 
  mutate(dummy = ifelse(faithful$eruptions>=3,1,0))

Podemos observar que ficou assim, pois criação de dummy

Código

table(banco$dummy) %>% 
  knitr::kable(col.names = c("Valores", "Quantidade"))

Valores	Quantidade
0	97
1	175

Ajustando os modelos

Código

fit_10 <- rq(waiting~dummy, tau = .10, data = banco)
fit_25 <- rq(waiting~dummy, tau = .25, data = banco)
fit_50 <- rq(waiting~dummy, tau = .50, data = banco)
fit_75 <- rq(waiting~dummy, tau = .75, data = banco)
fit_90 <- rq(waiting~dummy, tau = .90, data = banco)

Código

summary(fit_10)

Call: rq(formula = waiting ~ dummy, tau = 0.1, data = banco)

tau: [1] 0.1

Coefficients:
            coefficients   lower bd       upper bd      
(Intercept)   4.700000e+01 -1.797693e+308   4.857671e+01
dummy         2.600000e+01   2.338834e+01  1.797693e+308

Código

summary(fit_25)

Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
nonunique

Call: rq(formula = waiting ~ dummy, tau = 0.25, data = banco)

tau: [1] 0.25

Coefficients:
            coefficients lower bd upper bd
(Intercept) 50.00000     46.21107 55.28893
dummy       26.00000     22.10400 30.39600

Código

summary(fit_50)

Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
nonunique

Call: rq(formula = waiting ~ dummy, tau = 0.5, data = banco)

tau: [1] 0.5

Coefficients:
            coefficients lower bd upper bd
(Intercept) 54.00000     50.37215 57.62785
dummy       26.00000     22.98056 30.01944

Código

summary(fit_75)

Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
nonunique

Call: rq(formula = waiting ~ dummy, tau = 0.75, data = banco)

tau: [1] 0.75

Coefficients:
            coefficients lower bd upper bd
(Intercept) 59.00000     54.71107 62.78893
dummy       25.00000     19.20800 27.89600

Código

summary(fit_90)

Call: rq(formula = waiting ~ dummy, tau = 0.9, data = banco)

tau: [1] 0.9

Coefficients:
            coefficients lower bd upper bd
(Intercept) 63.00000     57.42329 65.35342
dummy       25.00000     23.58834 30.41166

Resultado

Depois de realizado o ajuste, podemos observar que não há influência da variável dummy com a variação dos quantis, ou seja, podemos observar que essa variável dummy não inluencia, concluindo assim que não existe diferença entre os quantis.

Concluímos então, que não há diferença entre os quantis. Se utilizado a variável dummy para ser utilizada. Existe apenas uma diferença no intercepto. O que de acordo com a visualização gráfica dos dados era esperado.

O que também pode ser comprovado, pelo gráfico abaixo:

Código

fm <- rq(waiting~dummy, data = banco, tau = 1:9/10)
plot(summary(fm))

Questão 6

• No R faça o que se pede

• 1. Usando o script a seguir gere dados com variância não constante e faça o gráfico usando o pacote “ggplot2”

Código

x <- seq(0,100,length.out = 100)
sig <- 0.1+0.05*x
b_0 <- 6
b_1 <- 0.1
set.seed(2023)
e <- rnorm(100,0,sd = sig)
y <- b_0+b_1*x+e
dat <- data.frame(x,y)
ggplot(dat, aes(x,y))+geom_point()

• 2. Usando a função “geom smoth” faça a regressão de y em x, plote a linha ajustada e adicione o intervalo de confiança.

Código

p <- ggplot(dat, aes(x,y))+
  geom_point()+
  geom_smooth(method = "lm")

p

`geom_smooth()` using formula = 'y ~ x'

• 3. Ajuste um modelo de regressão quantílica aos dados usando a opção “boot” para estimar o erro padrão. Avalie o modelo e faça o gráfico com as retas de regressão para os quantis avaliados.

Código

taus <- 1:9/10
qreg_fit <- rq(x~y, data = dat, tau = taus)
summary(qreg_fit, se = "boot")

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.1

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -23.72121   4.09946   -5.78642   0.00000
y             4.46852   0.51391    8.69507   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.2

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -31.54572   4.72023   -6.68308   0.00000
y             5.76228   0.58922    9.77952   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.3

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -30.26846   5.71263   -5.29851   0.00000
y             5.82402   0.52110   11.17644   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.4

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -26.74213   7.69220   -3.47652   0.00076
y             5.92117   0.64312    9.20696   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.5

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -19.34827   9.62120   -2.01100   0.04707
y             5.55993   0.85103    6.53317   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.6

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -11.37484  14.99429   -0.75861   0.44991
y             5.36082   1.20457    4.45041   0.00002

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.7

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept)  3.23945 22.96509    0.14106  0.88811
y            5.15432  1.62782    3.16640  0.00206

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.8

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept) 40.66269 19.44322    2.09136  0.03908
y            2.83560  1.37691    2.05939  0.04211

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.9

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept) 76.94906 18.73963    4.10622  0.00008
y            0.87712  1.27620    0.68729  0.49352

Código

p <- ggplot(dat, aes(x,y)) +
  geom_point()+
  geom_quantile(quantiles = taus)

p

Smoothing formula not specified. Using: y ~ x

• 4. Compare o resultado obtido com a regressão quantílica com o resultado obtido da regressão linear. O que você pode concluir? Você recomendaria a regressão quantílica nessa situação?

Resposta

Nessa situação através da análise gráfica, podemos observar que a regressão quantílica é mais recomendada, por capturar bem a variância não constante, e por não ser extremamente influenciada por pontos discrepantes, como é a regressão linear.

• 5. O pacote “quantreg” tem um método gráfico para visualizar a mudança nos coeficientes ao longo dos quantis. A linha vermelha apresenta a estimativa de MQO. O argumento parm = “x” indica que iremos avaliar o intercepto. Faça esses gráficos para avaliar o intercepto e o coeficiente da variável independente. O que você pode concluir?

Código

plot(summary(qreg_fit))

Resposta

A partir da análise gráfica acima, podemos observar que de fato a regressão quantilíca consegue captar melhor a variação do intercepto e do coeficiente ao longo da variação do quantil, ou seja, existe valores diferentes para quantis diferentes. Sendo assim, recomendável a utilização de tal.

• 6. Vamos considerar agora erros normais com variância constante. Faça os mesmos gráficos do item anterior e compare os resultados

Código

x <- seq(0,100,length.out = 100)
b_0 <- 6
b_1 <- 0.1
set.seed(2023)
e <- rnorm(100,0,sd = 1)
y <- b_0+b_1*x+e
dat <- data.frame(x,y)

taus <- 1:9/10
qreg_fit <- rq(x~y, data = dat, tau = taus)
summary(qreg_fit, se = "iid")

Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.1

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -55.51176   5.73985   -9.67128   0.00000
y             8.45602   0.50051   16.89490   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.2

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -59.00619   4.47942  -13.17274   0.00000
y             9.15340   0.39060   23.43428   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.3

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -53.95340   4.18666  -12.88697   0.00000
y             8.90387   0.36507   24.38942   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.4

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -50.59019   4.72714  -10.70207   0.00000
y             8.76598   0.41220   21.26632   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.5

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -49.55559   4.00920  -12.36047   0.00000
y             9.05694   0.34960   25.90685   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.6

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -51.42752   3.38468  -15.19422   0.00000
y             9.41435   0.29514   31.89802   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.7

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -51.29791   3.53484  -14.51209   0.00000
y             9.58977   0.30823   31.11210   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.8

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -48.71957   7.26245   -6.70842   0.00000
y             9.62350   0.63328   15.19639   0.00000

Call: rq(formula = x ~ y, tau = taus, data = dat)

tau: [1] 0.9

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -47.93468  10.36248   -4.62579   0.00001
y             9.97068   0.90359   11.03448   0.00000

Código

plot(summary(qreg_fit))

Resposta

A partir da análise gráfica acima, podemos observar que poderiamos utilizar a regressão linear sem muita perca de significância e de resultado, pois, podemos observar a proximidade dos valores, tanto do intercepto como também do coeficiente.

Questão 7

• Os dados “engel”(quntreg) consistem de observações sobre as despesas com alimentos e renda familiar de 235 trabalhadores da Bélgica (Ernst Engel, 1857): “http://www.econ.uiuc.edu/ roger/courses/LSE/data/

• 1. Faça o gráfico para verificar a relação entre as variáveis. O que você pode concluir? Parece haver heterocedasticidade?

Descrição das variáveis

Income: Annual household income in Belgian francs

Foodexp: Annual household food expenditure in Belgian francs

Código

library(quantreg)

Código

data(engel)
DT::datatable(engel, caption = "Banco de Dados")

Código

plot_ly(x = engel$income, y = engel$foodexp, text = paste("Income:", engel$income, "\nFoodexp", engel$foodexp)) %>% 
  layout(title = "Gráfico de Dispersão", xaxis = list(title = "Income"), yaxis = list(title = "Foodexp"))

Resposta

Podemos observar através da análise gráfica que existe sim problema com a heterocedasticidade, pois, a variância vai aumentando.

• 2. Faça a regressão linear e a regressão quantílica na mediana e acrescente as retas dos ajustes ao gráfico. O que você pode concluir? Você usou erros “i.ni.d” ou “i.i.d” (use o comando anova para testar a hipótese de igualdade)? Justifique sua resposta.

Código

fit_lm <- lm(foodexp~income, data = engel)
fit_5 <- rq(foodexp~income, data = engel, tau = 0.5)

Código

plot(x = engel$income, y = engel$foodexp, cex = .25, type = "n", xlab = "Income", ylab = "Foodexp")+
  points(engel$income, engel$foodexp, cex = .5, col = "blue")+
  abline(rq(foodexp~income, tau = .5, data = engel), col = "blue")+
  abline(lm(foodexp~income, data = engel), lty = 2, col = "red")+
  abline(rq(foodexp~income, tau = .25, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .10, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .75, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .90, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .95, data = engel), col = "gray")

integer(0)

Código

summary(fit_5, se = "iid")

Call: rq(formula = foodexp ~ income, tau = 0.5, data = engel)

tau: [1] 0.5

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept) 81.48225 13.23908    6.15468  0.00000
income       0.56018  0.01192   46.99766  0.00000

Código

summary(fit_5, se = "nid")

Call: rq(formula = foodexp ~ income, tau = 0.5, data = engel)

tau: [1] 0.5

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept) 81.48225 19.25066    4.23270  0.00003
income       0.56018  0.02828   19.81032  0.00000

Resposta

De acordo com as análises acima, podemos observar que a utilização de regressão quantilíca é mais adequada nesse caso, e que os erros a serem utilizado é o erro iid, como podemos analisar.

• 3. Faça a regressão quantílica para os quantis .05, .1, .25, .75, .90, .95 e acrescente as retas ajustadas ao gráfico. O que você pode concluir? Você recomendaria a regressão quantílica nesse caso? Avalie o ajuste do modelo.

Código

fit_05 <- rq(foodexp~income, data = engel, tau = .05)
fit_1 <- rq(foodexp~income, data = engel, tau = .1)
fit_25 <- rq(foodexp~income, data = engel, tau = .25)
fit_75 <- rq(foodexp~income, data = engel, tau = .75)
fit_90 <- rq(foodexp~income, data = engel, tau = .90)
fit_95 <- rq(foodexp~income, data = engel, tau = .95)

Código

plot(x = engel$income, y = engel$foodexp, cex = .25, type = "n", xlab = "Income", ylab = "Foodexp")+
  points(engel$income, engel$foodexp, cex = .5, col = "blue")+
  abline(rq(foodexp~income, tau = .05, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .25, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .10, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .75, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .90, data = engel), col = "gray")+
  abline(rq(foodexp~income, tau = .95, data = engel), col = "gray")

integer(0)

Resposta

De acordo com a visualização gráfica acima, conclui-se que a utilização de regressão quantilíca é adequado e que os erros possivelmente podem ser considerados como iid

• 4. Faça agora a regressão quantílica dos logariitmos na base 10 das variáveis y e x. Refaça o gráfico considerando os quantis do item anterior. O que você pode concluir? Você usou erros “i.ni.d” ou “i.i.d”? Justifique sua resposta.

Código

fit_05 <- rq(log10(foodexp)~log10(income), data = engel, tau = .05)
fit_1 <- rq(log10(foodexp)~log10(income), data = engel, tau = .1)
fit_25 <- rq(log10(foodexp)~log10(income), data = engel, tau = .25)
fit_75 <- rq(log10(foodexp)~log10(income), data = engel, tau = .75)
fit_90 <- rq(log10(foodexp)~log10(income), data = engel, tau = .90)
fit_95 <- rq(log10(foodexp)~log10(income), data = engel, tau = .95)

Código

plot(x = engel$income, y = engel$foodexp, cex = .25, type = "n", xlab = "Income", ylab = "Foodexp")+
  points(engel$income, engel$foodexp, cex = .5, col = "blue")+
  abline(rq(log10(foodexp)~log10(income), tau = .05, data = engel), col = "gray")+
  abline(rq(log10(foodexp)~log10(income), tau = .25, data = engel), col = "gray")+
  abline(rq(log10(foodexp)~log10(income), tau = .10, data = engel), col = "gray")+
  abline(rq(log10(foodexp)~log10(income), tau = .75, data = engel), col = "gray")+
  abline(rq(log10(foodexp)~log10(income), tau = .90, data = engel), col = "gray")+
  abline(rq(log10(foodexp)~log10(income), tau = .95, data = engel), col = "gray")

integer(0)

Resposta

No caso acima, podemos observar que o ajuste não ficou adequado aos dados, podemos observar que as retas não conseguiram captar bem a variação dos dados, então nesse caso, o ideal seria utilizar os erros iid, pois é observado que quase não há variação nos coeficientes, ou seja, se realizado o teste podemos comprovar que não há variação.

Questão 8

Considere o banco de dados “barro” do R (Barro e Lee, 1994). São dados sobre as taxas de crescimento do PIB para vários países. A variável dependente é y.net (variação anual do PIB per capita).

Código

data(barro)