3

Meu interesse é realizar um histograma da variável resposta com a densidade estimada conjuntamente com a densidade empírica de três modelos de regressão ajustados pelo pacote gamlss, onde esses modelos possuem covariáveis. A ideia é obter um gráfico “próximo” ao exemplo mostrado abaixo.

inserir a descrição da imagem aqui

Minha tentativa de realizar tal procedimento foi, no máximo, traçar a densidade empírica dessas distribuições, utilizando os seguintes códigos:

library(gamlss)
data(Orange)
mod1 = gamlss(circumference ~ age, 
              family=NO(mu.link = log), data = Orange)
mod2 = gamlss(circumference ~ age, 
              family=GA, data = Orange)
mod3 = gamlss(circumference ~ age, 
              family=WEI, data = Orange)
ymod1 <- predict(mod1)
ymod2 <- predict(mod2)
ymod3 <- predict(mod3)
dens1 <- density(exp(ymod1))
dens2 <- density(exp(ymod2))
dens3 <- density(exp(ymod3))

x11()
par(cex = 1.1)
hist(Orange$circumference, freq = FALSE,
     xlab = "Circunferência", ylab = "Densidade", main=" ")
lines(dens1, col = "red", lwd = 2)
lines(dens2, col = "blue", lwd = 2)
lines(dens3, col = "green", lwd = 2)
legend(x = "topright", legend = c("Normal", "Gama", "Weibull"), 
       lty = c(1,1,1),
       col = c("red", "blue", "green"))

inserir a descrição da imagem aqui

No entanto, observe que, como mencionei, as densidades realizadas acima são empíricas, pois (a densidade de tais distribuições não são bimodais). Acredito que seja necessário utilizar os parâmetros estimados, mas não sei como fazer isso e/ou talvez expressar para uma única observação tal fato levando em consideração a estrutura do preditor linear e as especificidades algébricas em termos dos parâmetros das distribuições (Ex: parâmetro de forma da distribuição Weibull é dado por outra expressão de acordo com a reparametrização utilizada no R), o mesmo para a distribuição Gama.

1 Resposta 1

3

Pode calcular média e desvio dos valores preditos e usá-los para estimar os parâmetros das funções de densidade de probabilidade, plotando-as para a extensão dos dados.

Não uso o pacote gamlss. Usarei glm para ilustrar o procedimento:

# Modelos
mod.n = glm(circumference ~ age, gaussian("log"), Orange)
mod.g = glm(circumference ~ age, Gamma("log"), Orange)

# Valores previstos
pred.n <- predict(mod.n, type = "r")
pred.g <- predict(mod.g, type = "r")
# type = "response" fornece os valores na escala da variável resposta

# Extensão dos dados de X
range(Orange$circumference)
#> [1]  30 214

# Densidade estimada pelos valores previstos
dens.n <- dnorm(0:240,
  mean = mean(pred.n),
  sd = sd(pred.n))
dens.g <- dgamma(0:240,
  shape = (mean(pred.g)/sd(pred.g))^2,
  rate = mean(pred.g)/sd(pred.g)^2)

hist(Orange$circumference, freq = FALSE, breaks = seq(0, 240, 20))
lines(dens.n, col = "red")
lines(dens.g, col = "blue")
legend("topright", legend = c("Normal", "Gamma"), lty = c(1,1), col = c("red", "blue"))

inserir a descrição da imagem aqui

Note que é importante fornecer a mesma extensão de X para o histograma e para as densidades.

Weibull

glm não possui suporte para Weibull. Mas assumindo o objeto pred.w com os valores preditos de um modelo ajustado com Weibull, aqui está como estimar os parâmetros λ (scale) e k (shape) a partir da média e desvio padrão:

k = (sd(pred.w)/mean(pred.w))^-1.086
s = mean(pred.w)/gamma(1 + 1/k)

dens.w <- dweibull(0:240, shape = k, scale = s)

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 .