3

Quero juntar os três gráficos abaixo usando a função plot. Eu faço o básico iniciando com o comando par(mfrow=c(1,3)) mas nao consigo junta-los.

O problema acredito que seja o fato de estar usando um objeto summary direto da regressão quantilica.

Vocês poderiam me ajudar nessa?

library(quantreg)
data(engel)

x1 <- engel$income
x2 <- 0.5*engel$income*2^(0.3)
x3 <- engel$income*3^(2)

fit1 <- rq(foodexp ~ x1, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)    
fit2 <- rq(foodexp ~ x2, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)    
fit3 <- rq(foodexp ~ x3, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)

fit1_betas <- summary(fit1, se="iid")
fit2_betas <- summary(fit2, se="iid")
fit3_betas <- summary(fit3, se="iid")

par(mfrow=c(1,3)) 

plot(fit1_betas, parm=2, main=expression("Regression": beta[0]), ylab="beta", xlab = "tau", cex=1, pch=20)

plot(fit2_betas, parm=2, main=expression("Regression": beta[0]), ylab="beta", xlab = "tau", cex=1, pch=20)

plot(fit3_betas, parm=2, main=expression("Regression": beta[0]), ylab="beta", xlab = "tau", cex=1, pch=20)

1 Resposta 1

2

Consegui fazer no ggplot2, baseando no código que fiz para a sua outra pergunta: Alterar gráficos dos betas de uma regressão

O código não ficou bonito, mas funcionou. A parte que ficou mais feia é onde coloco aquelas linhas vermelhas, que na verdade são as estimativas de uma regressão linear normal.

Segue:

library(broom)
library(dplyr)
library(ggplot2)
library(quantreg)

data(engel)

x1 <- engel$income
x2 <- 0.5*engel$income*2^(0.3)
x3 <- engel$income*3^(2)

fit1 <- rq(foodexp ~ x1, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)    
fit2 <- rq(foodexp ~ x2, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)    
fit3 <- rq(foodexp ~ x3, tau=seq(0.05, 0.95, by = 0.05), method="br", data=engel)

lm_fit1 <- lm(foodexp ~ x1, data = engel) 
lm_fit2 <- lm(foodexp ~ x2, data = engel) 
lm_fit3 <- lm(foodexp ~ x3, data = engel) 

lm_data <- bind_rows(
  lm_fit1 %>% tidy() %>% mutate(id = "fit1"), 
  lm_fit2 %>% tidy() %>% mutate(id = "fit2"),
  lm_fit3 %>% tidy() %>% mutate(id = "fit3")
  ) %>%
  bind_cols(bind_rows(
    lm_fit1 %>% confint() %>% as.data.frame(),
    lm_fit2 %>% confint() %>% as.data.frame(),
    lm_fit3 %>% confint()%>% as.data.frame()
  )) %>%
  rename(low = `2.5 %`, high = `97.5 %`) %>%
  filter(term != "(Intercept)")

bind_rows(
  tidy(fit1) %>% mutate(id = "fit1"),
  tidy(fit2) %>% mutate(id = "fit2"),
  tidy(fit3) %>% mutate(id = "fit3")
) %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(x = tau, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_point() +
  geom_hline(data = lm_data, aes(yintercept = estimate), colour = "red") +
  geom_hline(data = lm_data, aes(yintercept = low), colour = "red", linetype = "dashed") +
  geom_hline(data = lm_data, aes(yintercept = high), colour = "red", linetype = "dashed") +
  facet_wrap(~id, scales = "free")

Tente entender o passo a passo do que foi feito. Segue também o gráfico obtido.

inserir a descrição da imagem aqui

2
  • Daniel, novamente muito obrigado! Uma correção (?) na ultima linha seria term´ ao invés de id´. Certo??
    – cassius
    14/09/2016 às 13:26
  • 1
    Que bom que ajudou! Acho que funciona com os dois, mas realamente com term fica melhor. 14/09/2016 às 13:28

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .