Bounded mixed regression models for the study of the bronchial and lung cancer mortality rate in Brazil
Giovanni Piccirilli
2024-04-03
The dataGitset used in this application is avaiable at https://github.com/Giovannipcl/Bronchial-and-lung-cancer-mortality-rate-in-Brazil.
The following packages are loaded for the application.
library(gamlss)
library(ggplot2)Dataset
This paper addresses the evolution of the mortality rate over time in Brazilian states with the use of a dataset containing such a rate separated into groups of sex, Age, and Municipal Human Development Index (MHDI) about Longevity.
More specifically, the aim is to propose a bounded longitudinal regression model considering different alternative distributions for modeling the mortality rate.
We elaborated a data set from Brazilian public sites which had information about 27 Brazilian states over the last 30 years, and the variables are mortality rate of bronchial and lung cancer; Sex (male or female), treated as a factor; four age groups: between 50 and 59 years: \(\textrm{Age}_1\), between 60 and 69 years: \(\textrm{Age}_2\), between 70 and 79 years: \(\textrm{Age}_3\), over 80 years \(\textrm{Age}_4\).
Because the Brazilian national public health system, called Unified Health System (SUS in Portuguese), was created in 1988 and each federal government has implemented different policies through the system, we include a time variable considering the following seven federal government intervals of time (91-94, 95-98, 99-2002, 2003-2006, 2007-2010, 2011-2014, 2015-2018 and 2019-2023.). The response variable is the average mortality rate for each period, and the mortality rate refers to the gross rate per 100 thousand inhabitants. Its expression is given by the total number of events in the population at risk (gross rate = 100,000 \(\times\) (number of occurrences in defined period/reference population for the defined period)).
To obtain the mortality rate in a scale in the (0,1) interval we apply the following transformation: \((y_{obs} - y_{min})/(y_{max} - y_{min})\). Where \(y_{max}\) and \(y_{min}\) are the maximum and minimum rates reported originally.
dataGit = read.csv("https://raw.githubusercontent.com/Giovannipcl/Bronchial-and-lung-cancer-mortality-rate-in-Brazil/main/DataFinal.csv",header = TRUE)
dataGit$id = as.factor(dataGit$id)
dataGit$Sex = as.factor(ifelse(dataGit$sexo == 1, "Male", "Female"))
names(dataGit)[c(4,5,6,11,12)] = c("Time", "SexN", "IDH", "Age", "Longevity")
head(dataGit)## X y id Time SexN IDH MHDI V1 V2 V3 Age Longevity
## 1 1 0.03845559 1 1 1 0.5165076 0.7032559 0 0 0 50-59 High
## 2 2 0.09294294 1 1 1 0.5165076 0.7032559 1 0 0 60-69 High
## 3 3 0.20552606 1 1 1 0.5165076 0.7032559 0 1 0 70-79 High
## 4 4 0.13227486 1 1 1 0.5165076 0.7032559 0 0 1 80- High
## 5 5 0.04291458 2 1 1 0.4790902 0.6460460 0 0 0 50-59 Medium
## 6 6 0.06848291 2 1 1 0.4790902 0.6460460 1 0 0 60-69 Medium
## state Sex
## 1 Acre Male
## 2 Acre Male
## 3 Acre Male
## 4 Acre Male
## 5 Alagoas Male
## 6 Alagoas MaleDescriptive analysis
The following histogram illustrate the response variable (mortality rate) resclaed to (0,1) interval.
p<-ggplot(dataGit, aes(x=y)) +
geom_histogram(color="black", fill="grey")+ theme_classic() + labs(title="Mortality rate histogram plot",x="Mortality rate", y = "Count")
print(p)
The following Figure displays the evolution over time from each state. The mortality rate increases in some states with medium MHDI for older groups, although high MHDI states mean stagnation.
print(ggplot(data = dataGit, aes(x = Time, y = y,group = id))+
facet_grid(Age ~ Sex)+
geom_line(lty= 2)+
labs(y = "Mortality rate",x = "Time",cex=1.2)+
theme_bw() +theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())+ theme(axis.text=element_text(size=16),
axis.title=element_text(size=16,face="bold"), strip.text = element_text(
size = 16)))
The two following boxplot’s depict the groups concerning the response variable.
dataGit_male = dataGit[dataGit$Sex == "Male",]
p <- ggplot(dataGit_male, aes(x = Age, y = y)) + geom_boxplot()
p + theme_bw() +theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())+ theme(axis.text=element_text(size=20),
axis.title=element_text(size=20,face="bold"), strip.text = element_text(
size = 20))
dataGit_female = dataGit[dataGit$Sex == "Female",]
p <- ggplot(dataGit_female, aes(x = Age, y = y)) + geom_boxplot()
p + theme_bw() + ylim(0,1) +theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())+ theme(axis.text=element_text(size=20),
axis.title=element_text(size=20,face="bold"), strip.text = element_text(
size = 20))
Running the models
The proposed model is given by \[ Y_{ij} \sim JSB(\mu_{ij}, \sigma_{ij}, \nu_{ij}), \] where \(y_{ij}\) is the observed bronchial and lung cancer mortality rate divided by 500 of state \(i\), \(i=1,\dots,27\) and intervals \(j=1,\ldots,8\).
Different models were fitted. The first model assumes a random intercept and covariates only for the location parameter. Model 4 considers some second-order interaction, time, sex, and age effects on the dispersion parameter. A random slope is considered, improving the flexibility and inferences of the model.
gen.Family(family = "GU", type = "logit")
gen.Family(family = "RG",type = "logit")
gen.Family(family = "exGAUS",type = "logit")
gen.Family(family = "SN1",type = "logit")
gen.Family(family = "PE",type = "logit")Model 1
m1GU<-gamlss(y~Sex+Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitGU, method = mixed(40,40))
m1RG<-gamlss(y~Sex + V1 + V2 + V3 + Time + random(id) + MHDI , data=dataGit,family=logitRG, method = mixed(40,40))
m1exGAUS<-gamlss(y~Sex + V1 + V2 + V3 + Time + random(id) + MHDI, data=dataGit,family=logitexGAUS, method = mixed(40,40))
m1SN1<-gamlss(y~Sex + V1 + V2 + V3 + Time + random(id) + MHDI, data=dataGit,family=logitSN1, method = mixed(40,40))
m1PE<-gamlss(y~Sex+Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitPE, method = mixed(40,40))
m1BE<-gamlss(y~Sex+Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=BE, method = mixed(40,40))Model 2
m2GU<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitGU, method = mixed(40,40))
m2RG<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI , data=dataGit,family=logitRG, method = mixed(40,40))
m2exGAUS<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitexGAUS, method = mixed(40,40))
m2SN1<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitSN1, method = mixed(40,40))
m2PE<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=logitPE, method = mixed(40,40))
m2BE<-gamlss(y~Sex*Time + V1 + V2 + V3 + random(id) + MHDI, data=dataGit,family=BE, method = mixed(40,40))Model 3
m3GU<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI, sigma.formula = ~ Time, data=dataGit,family=logitGU, method = mixed(40,40))
m3RG<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI , sigma.formula = ~ Time, data=dataGit,family=logitRG, method = mixed(40,40))
m3exGAUS<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI, sigma.formula = ~ Time, data=dataGit,family=logitexGAUS, method = mixed(40,40))
m3SN1<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI, sigma.formula = ~ Time, data=dataGit,family=logitSN1, method = mixed(40,40))
m3PE<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI, sigma.formula = ~ Time, data=dataGit,family=logitPE, method = mixed(40,40))
m3BE<-gamlss(y~Sex + Time + V1 + V2 + V3 + random(id) + MHDI, sigma.formula = ~ Time, data=dataGit,family=BE, method = mixed(40,40))Model 4
m4GU<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=logitGU, trace=TRUE,method = mixed(40,40))
m4RG<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=logitRG, trace=TRUE,method =mixed(40,40))
m4exGAUS<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=logitexGAUS, trace=TRUE,method = mixed(40,40))
m4SN1<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=logitSN1, trace=TRUE,method = mixed(40,40))
m4PE<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=logitPE, trace=TRUE,method = mixed(40,40))
m4BE<-gamlss(y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age,data=dataGit,family=BE, trace=TRUE,method = mixed(40,40))Criteriaaic = matrix(0,4,7)
Criteriaaic[1,] = c("Modelo 1",round(AIC(m1GU),2), round(AIC(m1RG),2), round( AIC(m1exGAUS),2), round(AIC(m1SN1),2), round(AIC(m1PE),2), round(AIC(m1BE),2))
Criteriaaic[2,] = c("Modelo 2",round(AIC(m2GU),2), round(AIC(m2RG),2), round( AIC(m2exGAUS),2), round(AIC(m2SN1),2), round(AIC(m2PE),2), round(AIC(m2BE),2))
Criteriaaic[3,] = c("Modelo 3",round(AIC(m3GU),2), round(AIC(m3RG),2), round( AIC(m3exGAUS),2), round(AIC(m3SN1),2), round(AIC(m3PE),2), round(AIC(m3BE),2))
Criteriaaic[4,] = c("Modelo 4",round(AIC(m4GU),2), round(AIC(m4RG),2), round( AIC(m4exGAUS),2), round(AIC(m4SN1),2), round(AIC(m4PE),2), round(AIC(m4BE),2))
Criteriaaic = as.data.frame(Criteriaaic)
names(Criteriaaic) = c("Modelo","Gumbel", "Reverse Gumbel", "exGAUS", "Skew Normal", "Power Exponential","Beta")
Criteriabic = matrix(0,4,7)
Criteriabic[1,] = c("Modelo 1",round(BIC(m1GU),2), round(BIC(m1RG),2), round( BIC(m1exGAUS),2), round(BIC(m1SN1),2), round(BIC(m1PE),2), round(BIC(m1BE),2))
Criteriabic[2,] = c("Modelo 2",round(BIC(m2GU),2), round(BIC(m2RG),2), round( BIC(m2exGAUS),2), round(BIC(m2SN1),2), round(BIC(m2PE),2), round(BIC(m2BE),2))
Criteriabic[3,] = c("Modelo 3",round(BIC(m3GU),2), round(BIC(m3RG),2), round( BIC(m3exGAUS),2), round(BIC(m3SN1),2), round(BIC(m3PE),2), round(BIC(m3BE),2))
Criteriabic[4,] = c("Modelo 4",round(BIC(m4GU),2), round(BIC(m4RG),2), round( BIC(m4exGAUS),2), round(BIC(m4SN1),2), round(BIC(m4PE),2), round(BIC(m4BE),2))
Criteriabic = as.data.frame(Criteriabic)
names(Criteriabic) = c("Modelo","Gumbel", "Reverse Gumbel", "exGAUS", "Skew Normal", "Power Exponential","Beta")Criteriaaic## Modelo Gumbel Reverse Gumbel exGAUS Skew Normal Power Exponential
## 1 Modelo 1 -5001.58 -5800.47 -6000.75 -5811.27 -6117.35
## 2 Modelo 2 -5131.08 -5876.45 -6065.38 -5885.65 -6174.35
## 3 Modelo 3 -5305.66 -5820.81 -6027.33 -5907.98 -6157.19
## 4 Modelo 4 -6500.47 -6808.64 -6883.84 -6869.67 -6935.3
## Beta
## 1 -5215.78
## 2 -5265.61
## 3 -5322.66
## 4 -6708.65Criteriabic## Modelo Gumbel Reverse Gumbel exGAUS Skew Normal Power Exponential
## 1 Modelo 1 -4820.65 -5615.14 -5809.29 -5621.07 -5923.12
## 2 Modelo 2 -4944.22 -5685.06 -5868.45 -5689.83 -5974.98
## 3 Modelo 3 -5117.12 -5631.14 -5830.67 -5712 -5957.95
## 4 Modelo 4 -6143.29 -6446.59 -6500.02 -6489.22 -6544.5
## Beta
## 1 -5028.63
## 2 -5072.98
## 3 -5130.19
## 4 -6347.12The AIC and BIC criteria are smaller for Model 4 with Power Exponential. Besides that, except for the Gumbel model, all distributions presented smaller criteria than the Beta distribution.
Residuals
Quantile residuals
The randomized quantiles were computed for model 4 with Logistic Power distribution. The QQplot of such residuals was evaluated, and the Shapiro-Wilk test for the normality.
rq = residuals(m4PE)
shapiro.test(rq)##
## Shapiro-Wilk normality test
##
## data: rq
## W = 0.99808, p-value = 0.04par(mfrow = c(1,2))
hist(rq, main = "Histogram of quantile residuals", xlab = "Quantile residuals")
qqnorm(residuals(m4PE),xlim=c(-4,4), ylim = c(-4,5),cex.axis = 1.5,main = "",cex.lab = 1.5,cex = 1.5)
qqline(residuals(m4PE), lwd = 3)
Envelope for quantile residuals
###############################
rq = residuals(m4PE)
B <- 19
n = length(dataGit$y)
e <- matrix(0,n,B)
tr = fitted.values(m4PE)
sr = m4PE$sigma.fv
vr = m4PE$nu.fv
for (i in 1:B){
set.seed(123*i)
#Sim_Y <- rBE(n,tr,sr)
Sim_Y = rlogitPE(n, tr,sr,vr)
X = data.frame(cbind(Sim_Y,dataGit))
head(X)
attach(X)
#m1PE2<-gamlss(Sim_Y~sexo*tempo + V1*tempo + V2*tempo + V3*tempo +tempo + Longevidade*tempo +re(random=~tempo|id), sigma.formula = ~ tempo,data=y,family=BE, trace=TRUE, method = mixed(20,20))
#m1PE<-gamlss(Sim_Y~sexo*tempo + Idade*tempo +tempo + Longevidade*tempo +re(random=~tempo|id), sigma.formula = ~ tempo + sexo*Idade,nu.formula = ~ tempo,data=X,family=logitPE, trace=TRUE,method = RS(200))
m1PE<-gamlss(Sim_Y~Sex*Age +Time + MHDI +re(random=~Time|id), sigma.formula = ~ Time + Sex + Age, nu.formula = ~Age, data=X,family=logitPE, trace=FALSE,method = CG(120))
e[,i] <- sort(residuals(m1PE))
}
e1 <- numeric(n)
e2 <- numeric(n)
for (i in 1:n)
{
eo <- sort(e[i,])
e1[i] <- quantile(eo,0.025)
e2[i] <- quantile(eo,0.975)
}
rq = residuals(m4PE)
med <- apply(e,1,median)
td <- rq
faixa <- range(td,e1,e2)
e1g3 <- e1
e2g3 <- e2
tdg3 <- td
par(mar = c(5, 5, 1, 1)) # Set the margin on all sides to 6
qqnorm(td, xlab="Theoretical quantiles",
ylab="Quantile residuals", ylim=faixa, xlim=c(-3,3), pch=16,cex.axis = 2.1,main = "",cex.lab = 2.1,cex = 1)
par(new=TRUE)
#
#qqnorm(residuals(m1PE),xlim=c(-4,4), ylim = c(-4,4),cex.axis = 2.7,main = "",cex.lab = 2.7,cex = 2.5)
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa, xlim=c(-3,3),lty=1, lwd = 2,main="")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa, xlim=c(-3,3),lty=1,lwd = 2, main="")
par(new=TRUE)
qqnorm(med,axes=F,xlab="",ylab="", type="l",ylim=faixa, xlim=c(-3,3),lty=1, lwd = 2,main="")
Random effects
We can access the random effects for intercept and time for the states.
## State Intercept Time
## 1 Acre -0.087 0.056
## 2 Alagoas -0.568 0.079
## 3 Amap� 0.571 -0.124
## 4 Amazonas 0.64 -0.07
## 5 Bahia -0.695 0.052
## 6 Cear� -0.653 0.13
## 7 Distrito Federal 0.67 -0.121
## 8 Esp�rito Santo 0.319 -0.057
## 9 Goi�s 0.379 -0.036
## 10 Maranh�o -1.189 0.137
## 11 Mato Grosso 0.304 -0.057
## 12 Mato Grosso do Sul 0.521 -0.064
## 13 Minas Gerais 0.069 -0.047
## 14 Paran� 0.804 -0.077
## 15 Para�ba -1.233 0.169
## 16 Par� -0.17 -0.022
## 17 Pernambuco -0.237 0.029
## 18 Piau� -0.916 0.122
## 19 Rio Grande do Norte -0.839 0.111
## 20 Rio Grande do Sul 1.21 -0.047
## 21 Rio de Janeiro 0.61 -0.091
## 22 Rond�nia 0.179 -0.02
## 23 Roraima 0.23 -0.025
## 24 Santa Catarina 0.335 -0.022
## 25 Sergipe -0.141 0.005
## 26 S�o Paulo 0.501 -0.075
## 27 Tocantins -0.625 0.066Predicted values for each state
The following plot shows the estimated curve with the observations (dots) separated by state, sex and age.
mu = predict(m4PE,newdata = dataGit, what = "mu")
sigma = predict(m4PE,newdata = dataGit, what = "sigma")
nu = predict(m4PE,newdata = dataGit, what = "nu")
Sim_Y = lapply(1:(length(dataGit$y)), function(i) rlogitPE(10000, mu[i],exp(sigma[i]),exp(nu[i])))pred_Y = unlist(lapply(Sim_Y,mean))
dataGit$pred = pred_Y
# Choose a "i", in this case i = 18
data1 = dataGit[(dataGit$id == 18),]
print(ggplot(data = data1, aes(x = Time, y = pred, colour = Sex))+
facet_wrap(~Age)+
ggtitle(data1$state)+
geom_line()+ labs(x = "Time",cex=1.2)+
labs(y = "y",cex=1.2)+
geom_point(data=data1,aes(x=Time,y=y, colour = Sex),pch=1) +
theme_bw() +theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())+ theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold")))