#########################
#Entendendo um GAMLSS####
#########################
library(gamlss)
library(gamlss.cens)
library(gamlss.util)
library(gamlss.dist)

### Exemplo 1 #####
set.seed(30)
femea<-rnorm(60,40,10)
set.seed(30)
macho<-rnorm(60,80,10)
peso<-c(femea,macho)
genero<- c(rep('femea',60),  rep('macho',60))

dados<-data.frame(peso,genero)
head(dados)

#LM
model1=lm(peso~genero)
summary(model1)

#GAMLSS
model2<-gamlss(peso~genero)
summary(model2)
exp(2.39510)



#obtendo os parametros estimados
dados$genero[1]
model2$mu.fv[1]
model2$sigma.fv[1]

dados$genero[100]
model2$mu.fv[100]
model2$sigma.fv[100]

cbind(genero,media=model2$mu.fv,desvio.p=model2$sigma.fv)


hist(peso[genero=='femea'],probability = T)
curve(dnorm(x,model1$coefficients[1],11.06),add=T,col=2,lwd=4)#lm
curve(dnorm(x,model2$mu.fv[1],model2$sigma.fv[1]),add=T,col=4,lwd=4)#gamlss

hist(peso[genero=='macho'],probability = T)
curve(dnorm(x,sum(model1$coefficients),11.06),add=T,col=2,lwd=4)#lm
curve(dnorm(x,model2$mu.fv[120],model2$sigma.fv[120]),add=T,col=4,lwd=4)#gamlss






#Exemplo 2 - Variancia diferentes #####
set.seed(30)
femea<-rnorm(60,40,10)
set.seed(30)
macho<-rnorm(60,80,20)
peso<-c(femea,macho)
dados<-data.frame(peso,genero)
head(dados)


#LM
model1=lm(peso~genero)
summary(model1)

#GAMLSS
model2<-gamlss(peso~genero,sigma.formula =~genero,data=dados)
summary(model2)

cbind(genero,media=model2$mu.fv,desvio.p=model2$sigma.fv)


hist(peso[genero=='femea'],probability = T)
curve(dnorm(x,model1$coefficients[1],17.49),add=T,col=2,lwd=4)#lm
curve(dnorm(x,model2$mu.fv[1],model2$sigma.fv[1]),add=T,col=4,lwd=4)#gamlss

hist(peso[genero=='macho'],probability = T,ylim=c(0,0.026))
curve(dnorm(x,sum(model1$coefficients),17.49),add=T,col=2,lwd=4)#lm
curve(dnorm(x,model2$mu.fv[120],model2$sigma.fv[120]),add=T,col=4,lwd=4)#gamlss



summary(model1)
#SD = 2.25;   3.19
summary(model2)
#sd=  1.416 ; 3.16




#predicoes
model2$mu.fv
model2$sigma.fv
model2$sigma.fv^2


dados_novo<-data.frame(genero=c('macho','femea'))

predito<-predictAll(model2,newdata=dados_novo, type='response')
predito$mu
predito$sigma
predito$sigma^2





### Exemplo 3- Survival ####
library(gamlss.cens)


set.seed(10)
homem<-rWEI(300,2,2)
set.seed(10)
mulher<-rWEI(300,4,5)
hist(homem)
hist(mulher)

genero=c(rep('H',300),rep("M",300))
tempos=c(homem,mulher)
set.seed(20)
delta=rbinom(600,1,0.9)
dados<-data.frame(tempos,delta,genero)
head(dados)

km=survfit(Surv(tempos,delta)~genero)
plot(km,col=c(3,6),lwd=4)

#forma~escala~1
model1=gamlss(Surv(tempos,delta)~1,family = cens(WEI))
summary(model1)
curve(1-pWEI(x,model1$mu.fv[1],model1$sigma.fv[1]),add=T,lwd=3)

#forma~genero  escala~1
model2=gamlss(Surv(tempos,delta)~genero,family = cens(WEI))
summary(model2)
curve(1-pWEI(x,model2$mu.fv[1],model2$sigma.fv[1]),add=T,lwd=2,lty=2)
curve(1-pWEI(x,model2$mu.fv[600],model2$sigma.fv[600]),add=T,lwd=2,lty=2)



#forma~escala~genero
model3=gamlss(Surv(tempos,delta)~genero,sigma.fo=~genero,family = cens(WEI),data=dados)
summary(model3)
curve(1-pWEI(x,model3$mu.fv[1],model3$sigma.fv[1]),add=T,lwd=2,col=1)
curve(1-pWEI(x,model3$mu.fv[600],model3$sigma.fv[600]),add=T,lwd=2,col=1)


wp(model1,ylim.all = 1.5)
wp(model2,ylim.all = 1.5)
wp(model3,ylim.all = 1.5)


#### metodo de seleção
m0=gamlss(Surv(tempos,delta)~1,family = cens(WEI))
m1<-stepGAICAll.A(m0,scope=list(lower=~1,upper=~genero))
summary(m1)


#predicao
gen.cens(WEI)
dados_novo=data.frame(genero=c("H","M"))
predito<-predictAll(model3,newdata=dados_novo, type='response')
predito$mu
predito$sigma

#E(Y)
predito$mu*gamma(1/predito$sigma+1)









#### Exemplo 4 - non-lienar####
library(gamlss)
library(gamlss.util)
set.seed(1)
x1=runif(200,0,10)
plot(x1,6*sin(x1*pi/4),ylab='mu')
set.seed(1)
y=6*sin(x1*pi/4)+rnorm(200)

m1<-gamlss(y~x1) #mu linear
par(mfrow=c(1,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
term.plot(m1,parameter = 'mu')
par(mfrow=c(1,1))
wp(m1)


m2<-gamlss(y~pb(x1)) #mu nao linear
edfAll(m2) #gl>2 indica nao linear
par(mfrow=c(1,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
term.plot(m2,parameter = 'mu')
par(mfrow=c(1,1))
wp(m2)


AIC(m1,m2)
BIC(m1,m2)

#valores esperados
plot(x1,y)
lines(x1[order(x1)],m1$mu.fv[order(x1)],col=1,lwd=3)
lines(x1[order(x1)],m2$mu.fv[order(x1)],col=2,lwd=3)


dados=data.frame(y,x1)

#densidades por pontos de x1
plotSimpleGamlss(y,x1, model=m2, data=dados,
                 x.val=seq(1,10,1),val=1, N=1000)

#percentis
centiles(m2, xvar=x1, lwd=2,cent=c(5,50,95))






#### Exemplo 5 - non-lienar#########

set.seed(1)
x1=runif(200,0,10)
plot(x1,6*sin(x1*pi/4),ylab='mu')
plot(x1,0.5*x1,ylab='sigma')

set.seed(1)
y=6*sin(x1*pi/4)+0.5*x1*rnorm(200)


#mu nao linear e sigma constante
m2<-gamlss(y~pb(x1))
edfAll(m2) #gl>2 indica nao linear
term.plot(m2,parameter = 'mu')
wp(m2)

#mu nao linear e sigma linear
m3<-gamlss(y~pb(x1),sigma.fo=~x1)
edfAll(m3) #gl>2 indica nao linear
par(mfrow=c(2,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
term.plot(m3,parameter = 'mu')
plot(x1,0.5*x1,ylab='sigma')
term.plot(m3,parameter = 'sigma')
par(mfrow=c(1,1))
wp(m3)


AIC(m2,m3)
BIC(m2,m3)

#valores esperados
plot(x1,y)
#variancia constante
lines(x1[order(x1)],m2$mu.fv[order(x1)],col=1,lwd=3)
#variancia~x1(tempo)
lines(x1[order(x1)],m3$mu.fv[order(x1)],col=2,lwd=3)


dados=data.frame(y,x1)

par(mfrow=c(1,2))

#variancia constante
plotSimpleGamlss(y,x1, model=m2, data=dados,
 x.val=seq(1,10,1),val=2, N=1000)

#variancia~X1
plotSimpleGamlss(y,x1, model=m3, data=dados,
                 x.val=seq(1,10,1),val=2, N=1000)

#percentis

#variancia constante
centiles(m2, xvar=x1, lwd=2,cent=c(5,50,95))
#variancia~X1
centiles(m3, xvar=x1, lwd=2,cent=c(5,50,95))
par(mfrow=c(1,1))









#### Exemplo 6 - non-lienar#######

set.seed(1)
x1=runif(200,0,10)
par(mfrow=c(1,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
plot(x1,5*sin(x1*pi/6)+6,ylab='sigma')
par(mfrow=c(1,1))

set.seed(1)
y=6*sin(x1*pi/4)+(5*sin(x1*pi/6)+6)*rnorm(200)

#media nao linear e sigma constante
m2<-gamlss(y~pb(x1))
edfAll(m2) #gl>2 indica nao linear
term.plot(m2,parameter = 'mu')
wp(m2)

#media nao linear e sigma linear
m3<-gamlss(y~pb(x1),sigma.fo=~x1)
edfAll(m3) #2 indica ñ linear
par(mfrow=c(2,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
term.plot(m3,parameter = 'mu')
plot(x1,5*sin(x1*pi/6)+6,ylab='sigma')
term.plot(m3,parameter = 'sigma')
par(mfrow=c(1,1))
wp(m3)


#media nao linear e sigma nao linear
m4<-gamlss(y~pb(x1),sigma.fo=~pb(x1))
edfAll(m4) #gl>2 indica nao linear
par(mfrow=c(2,2))
plot(x1,6*sin(x1*pi/4),ylab='mu')
term.plot(m4,parameter = 'mu')
plot(x1,5*sin(x1*pi/6)+6,ylab='sigma')
term.plot(m4,parameter = 'sigma')
par(mfrow=c(1,1))
wp(m4,ylim.all = 2)

AIC(m2,m3,m4)
BIC(m2,m3,m4)

#Valores esperados
plot(x1,y)
lines(x1[order(x1)],m4$mu.fv[order(x1)],col=2,lwd=3)


dados=data.frame(y,x1)

#media nao linear e sigma constante
plotSimpleGamlss(y,x1, model=m2, data=dados,
   x.val=seq(1,10,1),val=4, N=1000)

#media nao linear e sigma linear
plotSimpleGamlss(y,x1, model=m3, data=dados,
     x.val=seq(1,10,1),val=4, N=1000)

#media nao linear e sigma nao linear
plotSimpleGamlss(y,x1, model=m4, data=dados,
 x.val=seq(1,10,1),val=4, N=1000)

#percentis

#media nao linear e sigma constante
centiles(m2, xvar=x1, lwd=2,cent=c(5,50,95))
centiles(m3, xvar=x1, lwd=2,cent=c(5,50,95))
centiles(m4, xvar=x1, lwd=2,cent=c(5,50,95))




#### Exemplo 7####
library(gamlss.dist)

set.seed(10)
a<-rPE2(200,5,1,2)
#mesocurtica
curve(dPE2(x,5,1,2),2,8,ylim=c(0,0.8),lwd=4,xlab='Y',ylab="Densidade")

set.seed(10)
b<-rPE2(200,5,1,1)
#leptocurtica
curve(dPE2(x,5,1,1),0,10,col='chartreuse1',lwd=4,add=T)

set.seed(10)
c<-rPE2(200,5,1,5)
#platicurtica
curve(dPE2(x,5,1,5),0,10,col='deepskyblue',lwd=4,add=T)

y=c(a,b,c)
tipo<-rep(c('a','b','c'),each=200)

dados<-data.frame(y,tipo)
head(dados)


#Normal
m1<-gamlss(y~tipo,sigma.fo=~tipo,data=dados)
#PE2
m2<-gamlss(y~tipo,sigma.fo=~tipo,nu.fo=~tipo,
        family = PE2,n.cyc=1500)


summary(m1)
summary(m2)

#Normal
m1<-gamlss(y~1,sigma.fo=~tipo,data=dados)
#PE2
m2<-gamlss(y~1,sigma.fo=~1,nu.fo=~tipo,family = PE2)



hist(a,probability = T,main = 'mesocurtica',ylim=c(0,0.8))
curve(dnorm(x,m1$mu.fv[1],m1$sigma.fv[1]),add=T,lwd=3,lty=1)
curve(dPE2(x,m2$mu.fv[1],m2$sigma.fv[1],m2$nu.fv[1]),add=T,lwd=4,lty=2,col='2')

hist(b,probability = T,main='leptocúrtica',ylim=c(0,0.8),breaks = 30)
curve(dnorm(x,m1$mu.fv[400],m1$sigma.fv[400]),add=T,lwd=3)
curve(dPE2(x,m2$mu.fv[400],m2$sigma.fv[400],m2$nu.fv[400]),add=T,lwd=3,lty=2,col=2)

hist(c,probability = T,main='platicúrtica',ylim=c(0,0.90))
curve(dnorm(x,m1$mu.fv[600],m1$sigma.fv[600]),add=T,lwd=3)
curve(dPE2(x,m2$mu.fv[600],m2$sigma.fv[600],m2$nu.fv[600]),col=2,add=T,lwd=4,lty=2)


AIC(m1,m2)
BIC(m1,m2)