Problema para plotar intervalo de confiança no ggplot2 para modelo GLM de Poisson com junção de níveis
Prezados Membros, Estou tentando plotar sem sucesso o intervalo em um modelo GLM de Poisson onde fiz a junção de níveis para a variável categórica Gender. Alguém baseado no CRM abaixo poderia dar uma luz? #Pacotes library(ggplot2) library(dplyr) library(tidyverse) #Dados inventados com uma variável resposta (Consumption) e duas explicativas, sendo uma categórica e outra quali Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80) Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE) Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512) df3 <- data.frame(Consumption, Gender, Income) df3 # GLM de Poisson fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson) summary(fm1) # ANOVA do modelo ajustado anova(fm1,test="Chi") #Comparo a variável Gender sort(tapply(df3$Consumption,df3$Gender,mean)) Gender2<-df3$Gender levels(Gender2) levels(Gender2)[2]<-"Fem_Her" levels(Gender2)[3]<-"Fem_Her" levels(Gender2) fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson) anova(fm1,fm2,test="Chi") # 0.7824 Female/Hermaphrodite são iguais então eu junto #Faço a predição sobre o modelo final e dos intervalos de confiança pred <- predict(fm2, type="response", se.fit = TRUE) df3 = cbind(df3, pred = pred$fit) df3 = cbind(df3, se = pred$se.fit) df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se) df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se) df3 = cbind(df3, Gender2) df<-df3 %>% dplyr::group_by(Income, Gender2) %>% dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE)) df<-as.data.frame(df) #Faço o plot usando o ggplot2 df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=type, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) + geom_line(mapping=aes(x=Income, y=pred)) # Obrigado, Alexandre -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ======================================================================
O problema está na linha ggplot(mapping=aes(x=type, y=value, color = Gender2)) + Não faz muito sentido criar um gráfico de dispersão com uma variável categórica de apenas um nível como variável preditora. Rode o código abaixo que ele deve ficar mais próximo do desejado. df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=Income, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity", colour = "black") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2), size = 2) + geom_line(mapping=aes(x=Income, y=pred), colour = "black") Eu achei esquisito que o nível Fem_Her está muito distante do nível Male no gráfico, mas como a income dos Male é em torno de 20 vezes a dos Fem_Her, então o plot faz sentido. Ah, troquei a cor da linha dos valores preditos porque ela estava sendo confundida com alguns pontos. Talvez preto não seja a melhor opção, então escolha uma cor que se harmonize melhor à tua paleta de cores. Atenciosamente, -- Marcus Nunes Professor Adjunto Universidade Federal do Rio Grande do Norte Centro de Ciências Exatas e da Terra Departamento de Estatística Laboratório de Estatística Aplicada marcus@marcusnunes.me https://marcusnunes.me/ http://shiny.estatistica.ccet.ufrn.br On Wed, Jun 19, 2019 at 11:45 AM ASANTOS por (R-br) < r-br@listas.c3sl.ufpr.br> wrote:
Prezados Membros,
Estou tentando plotar sem sucesso o intervalo em um modelo GLM de Poisson onde fiz a junção de níveis para a variável categórica Gender. Alguém baseado no CRM abaixo poderia dar uma luz?
#Pacotes library(ggplot2) library(dplyr) library(tidyverse)
#Dados inventados com uma variável resposta (Consumption) e duas explicativas, sendo uma categórica e outra quali Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80) Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE) Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512) df3 <- data.frame(Consumption, Gender, Income) df3
# GLM de Poisson fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson) summary(fm1)
# ANOVA do modelo ajustado anova(fm1,test="Chi")
#Comparo a variável Gender sort(tapply(df3$Consumption,df3$Gender,mean)) Gender2<-df3$Gender levels(Gender2) levels(Gender2)[2]<-"Fem_Her" levels(Gender2)[3]<-"Fem_Her" levels(Gender2) fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson) anova(fm1,fm2,test="Chi") # 0.7824 Female/Hermaphrodite são iguais então eu junto
#Faço a predição sobre o modelo final e dos intervalos de confiança
pred <- predict(fm2, type="response", se.fit = TRUE) df3 = cbind(df3, pred = pred$fit) df3 = cbind(df3, se = pred$se.fit) df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se) df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se) df3 = cbind(df3, Gender2)
df<-df3 %>% dplyr::group_by(Income, Gender2) %>% dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE)) df<-as.data.frame(df)
#Faço o plot usando o ggplot2 df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=type, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) + geom_line(mapping=aes(x=Income, y=pred))
#
Obrigado,
Alexandre
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
Bom dia Marcus e obrigado pela ajuda, Mas ainda tem um problema, pois agora com a inserção de ggplot(mapping=aes(x=Income, y=value, color = Gender2)), as linhas para o modelo predito ziguezagueiam ao invés de se comportarem como curvas. A ideia de utilizar a função gather() é justamente funcionar com a função order() quando eu fazia gráficos com a função plot() padrão e aparentemente isto esta bem difícil de ser resolvido. -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ====================================================================== Em 19/06/2019 11:08, Marcus Nunes escreveu:
O problema está na linha
ggplot(mapping=aes(x=type, y=value, color = Gender2)) +
Não faz muito sentido criar um gráfico de dispersão com uma variável categórica de apenas um nível como variável preditora. Rode o código abaixo que ele deve ficar mais próximo do desejado.
df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=Income, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity", colour = "black") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2), size = 2) + geom_line(mapping=aes(x=Income, y=pred), colour = "black")
Eu achei esquisito que o nível Fem_Her está muito distante do nível Male no gráfico, mas como a income dos Male é em torno de 20 vezes a dos Fem_Her, então o plot faz sentido.
Ah, troquei a cor da linha dos valores preditos porque ela estava sendo confundida com alguns pontos. Talvez preto não seja a melhor opção, então escolha uma cor que se harmonize melhor à tua paleta de cores.
Atenciosamente, -- Marcus Nunes Professor Adjunto Universidade Federal do Rio Grande do Norte Centro de Ciências Exatas e da Terra Departamento de Estatística Laboratório de Estatística Aplicada marcus@marcusnunes.me <mailto:marcus@marcusnunes.me> https://marcusnunes.me/ http://shiny.estatistica.ccet.ufrn.br
On Wed, Jun 19, 2019 at 11:45 AM ASANTOS por (R-br) <r-br@listas.c3sl.ufpr.br <mailto:r-br@listas.c3sl.ufpr.br>> wrote:
Prezados Membros,
Estou tentando plotar sem sucesso o intervalo em um modelo GLM de Poisson onde fiz a junção de níveis para a variável categórica Gender. Alguém baseado no CRM abaixo poderia dar uma luz?
#Pacotes library(ggplot2) library(dplyr) library(tidyverse)
#Dados inventados com uma variável resposta (Consumption) e duas explicativas, sendo uma categórica e outra quali Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80) Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE) Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512) df3 <- data.frame(Consumption, Gender, Income) df3
# GLM de Poisson fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson) summary(fm1)
# ANOVA do modelo ajustado anova(fm1,test="Chi")
#Comparo a variável Gender sort(tapply(df3$Consumption,df3$Gender,mean)) Gender2<-df3$Gender levels(Gender2) levels(Gender2)[2]<-"Fem_Her" levels(Gender2)[3]<-"Fem_Her" levels(Gender2) fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson) anova(fm1,fm2,test="Chi") # 0.7824 Female/Hermaphrodite são iguais então eu junto
#Faço a predição sobre o modelo final e dos intervalos de confiança
pred <- predict(fm2, type="response", se.fit = TRUE) df3 = cbind(df3, pred = pred$fit) df3 = cbind(df3, se = pred$se.fit) df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se) df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se) df3 = cbind(df3, Gender2)
df<-df3 %>% dplyr::group_by(Income, Gender2) %>% dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE)) df<-as.data.frame(df)
#Faço o plot usando o ggplot2 df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=type, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) + geom_line(mapping=aes(x=Income, y=pred))
#
Obrigado,
Alexandre
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br <mailto:e-mails%3Aalexandresantosbr@yahoo.com.br> alexandre.santos@cas.ifmt.edu.br <mailto:alexandre.santos@cas.ifmt.edu.br> Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 <http://orcid.org/0000-0001-8232-6722> - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 <http://www.researchgate.net/profile/Alexandre_Santos10> LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 <http://br.linkedin.com/in/alexandre-dos-santos-87961635> Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ <http://www.mendeley.com/profiles/alexandre-dos-santos6/> ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br <mailto:R-br@listas.c3sl.ufpr.br> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
A variância das variáveis é muito grande entre os homens. Não ficaria melhor transformar em log? daniel
On 19 Jun 2019, at 11:45, ASANTOS por (R-br) <r-br@listas.c3sl.ufpr.br> wrote:
Prezados Membros,
Estou tentando plotar sem sucesso o intervalo em um modelo GLM de Poisson onde fiz a junção de níveis para a variável categórica Gender. Alguém baseado no CRM abaixo poderia dar uma luz?
#Pacotes library(ggplot2) library(dplyr) library(tidyverse)
#Dados inventados com uma variável resposta (Consumption) e duas explicativas, sendo uma categórica e outra quali Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80) Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE) Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512) df3 <- data.frame(Consumption, Gender, Income) df3
# GLM de Poisson fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson) summary(fm1)
# ANOVA do modelo ajustado anova(fm1,test="Chi")
#Comparo a variável Gender sort(tapply(df3$Consumption,df3$Gender,mean)) Gender2<-df3$Gender levels(Gender2) levels(Gender2)[2]<-"Fem_Her" levels(Gender2)[3]<-"Fem_Her" levels(Gender2) fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson) anova(fm1,fm2,test="Chi") # 0.7824 Female/Hermaphrodite são iguais então eu junto
#Faço a predição sobre o modelo final e dos intervalos de confiança
pred <- predict(fm2, type="response", se.fit = TRUE) df3 = cbind(df3, pred = pred$fit) df3 = cbind(df3, se = pred$se.fit) df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se) df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se) df3 = cbind(df3, Gender2)
df<-df3 %>% dplyr::group_by(Income, Gender2) %>% dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE)) df<-as.data.frame(df)
#Faço o plot usando o ggplot2 df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=type, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) + geom_line(mapping=aes(x=Income, y=pred))
#
Obrigado,
Alexandre
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
Já houve boas interações, mas seu exemplo me faz questionar: você escreveu ". . . uma variável resposta (Consumption) e duas explicativas, *sendo uma categórica e outra quali* . . . " se a variável Income é qualitativa, porque você deixa o programa tratá-la como numérica? HTH -- Cesar Rabak On Wed, Jun 19, 2019 at 11:45 AM ASANTOS por (R-br) < r-br@listas.c3sl.ufpr.br> wrote:
Prezados Membros,
Estou tentando plotar sem sucesso o intervalo em um modelo GLM de Poisson onde fiz a junção de níveis para a variável categórica Gender. Alguém baseado no CRM abaixo poderia dar uma luz?
#Pacotes library(ggplot2) library(dplyr) library(tidyverse)
#Dados inventados com uma variável resposta (Consumption) e duas explicativas, sendo uma categórica e outra quali Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80) Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE) Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512) df3 <- data.frame(Consumption, Gender, Income) df3
# GLM de Poisson fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson) summary(fm1)
# ANOVA do modelo ajustado anova(fm1,test="Chi")
#Comparo a variável Gender sort(tapply(df3$Consumption,df3$Gender,mean)) Gender2<-df3$Gender levels(Gender2) levels(Gender2)[2]<-"Fem_Her" levels(Gender2)[3]<-"Fem_Her" levels(Gender2) fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson) anova(fm1,fm2,test="Chi") # 0.7824 Female/Hermaphrodite são iguais então eu junto
#Faço a predição sobre o modelo final e dos intervalos de confiança
pred <- predict(fm2, type="response", se.fit = TRUE) df3 = cbind(df3, pred = pred$fit) df3 = cbind(df3, se = pred$se.fit) df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se) df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se) df3 = cbind(df3, Gender2)
df<-df3 %>% dplyr::group_by(Income, Gender2) %>% dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE)) df<-as.data.frame(df)
#Faço o plot usando o ggplot2 df3 %>% tidyr::gather(type, value, Consumption) %>% ggplot(mapping=aes(x=type, y=value, color = Gender2)) + geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") + geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) + geom_line(mapping=aes(x=Income, y=pred))
#
Obrigado,
Alexandre
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
participantes (4)
-
ASANTOS -
Cesar Rabak -
Daniel Guimarães Tiezzi -
Marcus Nunes