---
title: "MCMC-glmm modelling of maturation probability - J-P Verta"
output:
    html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library("MCMCglmm")
library("ggplot2")
library("plyr")
library("ggpubr")
library("ggsignif")
```

Functions to work with probit scale.
```{r}
probitlink = function (theta, bvalue = NULL, inverse = FALSE, deriv = 0, short = TRUE, 
    tag = FALSE) 
{
    if (is.character(theta)) {
        string <- if (short) 
            paste("probitlink(", theta, ")", sep = "")
        else paste("qnorm(", theta, ")", sep = "")
        if (tag) 
            string <- paste("Probit:", string)
        return(string)
    }
    if (!inverse && length(bvalue)) {
        theta[theta <= 0] <- bvalue
        theta[theta >= 1] <- 1 - bvalue
    }
    if (inverse) {
        switch(as.character(deriv), `0` = {
            ans <- pnorm(theta)
            if (is.matrix(theta)) dim(ans) <- dim(theta)
            ans
        }, `1` = {
            1/Recall(theta = theta, bvalue = bvalue, inverse = FALSE, 
                deriv = deriv)
        }, `2` = {
            Junk <- qnorm(theta)
            ans <- -Junk * dnorm(Junk)
            if (is.vector(theta)) ans else if (is.matrix(theta)) {
                dim(ans) <- dim(theta)
                ans
            } else {
                warning("can only handle vectors and matrices;", 
                  " converting to vector")
                ans
            }
        }, `3` = {
            Junk <- qnorm(theta)
            junk <- dnorm(Junk)
            junk * (Junk^2 - 1)
        }, stop("argument 'deriv' unmatched"))
    }
    else {
        switch(as.character(deriv), `0` = {
            ans <- qnorm(theta)
            if (is.matrix(theta)) dim(ans) <- dim(theta)
            ans
        }, `1` = {
            if (is.matrix(theta)) {
                ans <- 1/dnorm(qnorm(theta))
                dim(ans) <- dim(theta)
                ans
            } else {
                1/dnorm(qnorm(as.vector(theta)))
            }
        }, `2` = {
            Junk <- qnorm(theta)
            ans <- Junk/(dnorm(Junk))^2
            if (is.vector(theta)) ans else if (is.matrix(theta)) {
                dim(ans) <- dim(theta)
                ans
            } else {
                warning("can only handle vectors and matrices;", 
                  " converting to vector")
                ans
            }
        }, `3` = {
            Junk <- qnorm(theta)
            junk <- dnorm(Junk)
            (1 + 2 * Junk^2)/junk^3
        }, stop("argument 'deriv' unmatched"))
    }
}

### back transformation logit
probit.inv = function(x) {
 		probitlink(x, inverse=TRUE, bvalue = .Machine$double.eps)
  		}
		
probit= function(x) {
 		probitlink(x, inverse=FALSE, bvalue = .Machine$double.eps)
  		}
```

Load data.
```{r}
load("MCMCglmmMATmodelPLOSpublic.RData")
```

Format data. Length at P1 is the only phenotypic measurement used in the model (apart from maturation status).
```{r}
MaleData_P2_OUL$DNA.ID.at.Tagging = as.factor(MaleData_P2_OUL$DNA.ID.at.Tagging)
MaleData_P2_OUL$VGLL3TOP.RECPRO = as.factor(MaleData_P2_OUL$VGLL3TOP.RECPRO)
MaleData_P2_OUL$ln.length_P1.scaled = scale(log(MaleData_P2_OUL$length.cm_P1))
```


## Test for vgll3 effects (without length)

Define prior (default prior from: de Villemereuil, P., O. Gimenez, B. Doligez, and R. Freckleton. 2013. Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits. Methods in Ecology and Evolution 4:260-275)
```{r}
prior.MatModel_OUL = NULL
prior.MatModel_OUL = list(
	R = list(V = diag(1), fix= 1), 
	G = list(G1 = list(V= diag(1) , nu= 1000, alpha.mu= c(0), alpha.V= diag(1)),
			 G3 = list(V= diag(1) , nu= 1000, alpha.mu= c(0), alpha.V= diag(1)))) 
prior.MatModel_OUL
```

## Collapse heterozygotes and run models with and without length

Create a new column where heterozygotes are collapsed.
```{r}
MaleData_P2_OUL$VGLL3TOP.COLLAPSED = MaleData_P2_OUL$VGLL3TOP.RECPRO
MaleData_P2_OUL$VGLL3TOP.COLLAPSED = revalue(MaleData_P2_OUL$VGLL3TOP.COLLAPSED, c("EL"="EL/LE","LE"="EL/LE"))
unique(MaleData_P2_OUL$VGLL3TOP.COLLAPSED)
```

Run model without length on collapsed heterozygotes.
```{r}
model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed = NULL
model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed = MCMCglmm(
	fixed= mature.bin ~ 1 + VGLL3TOP.COLLAPSED ,
	random= ~ Tank + DNA.ID.at.Tagging,
	rcov=~ units,
	data= MaleData_P2_OUL,
 	trunc= TRUE,
	# pr = TRUE,
	ginverse= list(DNA.ID.at.Tagging = OUL.WARM_Male.Pedigree.GINV.MCMCglmm ),
	family= "threshold",
	prior= prior.MatModel_OUL,
		nitt=  10000000 + 100000,
		burnin=  10000,
		thin=     1000)
```


Check if model converges ("passed" in the stationary test).
```{r}
heidel.diag(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed[["VCV"]])
heidel.diag(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed[["Sol"]])
```

Summary of results.
```{r}
summary(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed)
```

Quality check of model.
```{r}
plot(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol)
plot(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$VCV)
```

Check that autocorrelation structure in the data is OK (below 10% is fine).
```{r}
#autocorr.MCMCglmm(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed)
```

Compare mean maturation probability at average length.
```{r}
#### means on proportion scale:
model = NULL
model = model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed
newdata = NULL
newdata = MaleData_P2_OUL[,which(
	names(MaleData_P2_OUL) == "mature.bin" |
	names(MaleData_P2_OUL) == "VGLL3TOP.COLLAPSED" |
	names(MaleData_P2_OUL) == "Tank" |
	names(MaleData_P2_OUL) == "DNA.ID.at.Tagging" )]
	
library("MCMCglmm")
predict.model = predict.MCMCglmm(
		newdata = newdata,
		object= model, 
		marginal= model$Random$formula, 
		type= "response", 
		interval= "confidence")
# yields mean lower upper, but many repated data -> use unique

predict.model = unique(data.frame(round(predict.model,4), newdata))
names(predict.model)[1:3] = c("mean","lower","upper")
head(predict.model)
meanMat = unique(predict.model[,c("mean","lower","upper","VGLL3TOP.COLLAPSED")])


# change probit contrast to direction EE-LL
data.frame(VGLL3.Contrast = c("EE-LL"),
 contrast = mean(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"])),
 	lower = HPDinterval(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"]))[1],
 	upper =	HPDinterval(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"]))[2])

data.frame(VGLL3.Contrast = c("EE-EL"),
 contrast = mean(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"])),
 	lower = HPDinterval(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]))[1],
 	upper =	HPDinterval(model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]))[2])
```

Compare maturation probabilities between genotypes (others than EE constrasts). Output is 95% credible intervals - if the interval includes zero, the models are not different at *P*<0.05.
```{r}
# EL - LL
model = (model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]  - 
		model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"])
# confidence interval
round(cbind(mean = mean(model),HPDinterval(model)),2)
# approximate p-value
2*length(which(model>0))/length(model)
```

Plot maturation probabilities at P1.
```{r}
ggplot(meanMat,aes(x=VGLL3TOP.COLLAPSED,y=mean,colour=VGLL3TOP.COLLAPSED)) + 
  geom_point(size=10) + 
  scale_color_manual(values=c("blue","darkgrey","red")) +
  geom_errorbar(width=.0, aes(ymin=lower, ymax=upper),size=1.5) +
  ylab("Mean maturation probability") + 
  xlab("vgll3 genotype") +
  theme_classic() +
  theme(axis.text=element_text(size=14,color="black"),axis.title=element_text(size=14,color="black"))
```

Plot maturation frequencies.
```{r}
matTable = table(MaleData_P2_OUL$mature.bin,MaleData_P2_OUL$VGLL3TOP.COLLAPSED)
plotData = as.data.frame.matrix(t(matTable))
colnames(plotData) = c("Immature", "Early")
plotData$pct = apply(plotData,1,function(x){x["Early"]/(x["Early"]+x["Immature"])*100})
plotData$vgll3 = rownames(plotData)

ggplot(plotData,aes(x=vgll3,y=pct)) + 
  geom_bar(stat="identity",fill="black") + 
  ylab("Mature parr at 9 months (%)") + 
  guides(fill=guide_legend(title="Maturation status")) +
  scale_x_discrete(labels=c("vgll3*EE","vgll3*EL", "vgll3*LL")) +
  theme_classic() +
  geom_text(aes(label=paste(Early,"/",Early+Immature,sep='')), vjust=-0.25) +
  theme(axis.text.x=element_text(size=14,color=c("blue","darkgrey", "red"),face="italic"),
        axis.title.x=element_blank(),
        axis.title.y=element_text(size=14,color="black"),
        axis.text.y=element_text(size=14,color="black"),
        legend.position = "none")
```

Model with heterozygotes collapsed and length as co-variate.
```{r}
model.MCMCglmm_OUL.FullRandom.HetsCollapsed = NULL
model.MCMCglmm_OUL.FullRandom.HetsCollapsed= MCMCglmm(
	fixed= mature.bin ~ 1 + VGLL3TOP.COLLAPSED + ln.length_P1.scaled,
	random= ~ Tank
			+ DNA.ID.at.Tagging,
	rcov=~ units,
	data= MaleData_P2_OUL,
 	trunc= TRUE,
	# pr = TRUE,
	ginverse= list(DNA.ID.at.Tagging = OUL.WARM_Male.Pedigree.GINV.MCMCglmm ),
	family= "threshold",
	prior= prior.MatModel_OUL,
		nitt=  10000000 + 100000,
		burnin=  10000,
		thin=     1000)
```

Check if model converges ("passed" in the stationary test).
```{r}
heidel.diag(model.MCMCglmm_OUL.FullRandom.HetsCollapsed[["VCV"]])
heidel.diag(model.MCMCglmm_OUL.FullRandom.HetsCollapsed[["Sol"]])
```

Summary of results.
```{r}
summary(model.MCMCglmm_OUL.FullRandom.HetsCollapsed)
```

Quality check of model.
```{r}
plot(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol)
plot(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$VCV)
```

Check that autocorrelation structure in the data is OK (below 10% is fine).
```{r}
#autocorr.MCMCglmm(model.MCMCglmm_OUL.FullRandom.HetsCollapsed)
```

#	Compare mean maturation probability at average length.
```{r}
#### means on proportion scale:
model = NULL
#model = model.MCMCglmm_OUL.NoLength.NoFeed.HetsCollapsed
model = model.MCMCglmm_OUL.FullRandom.HetsCollapsed
newdata = NULL
newdata = MaleData_P2_OUL[,which(
	names(MaleData_P2_OUL) == "mature.bin" |
	names(MaleData_P2_OUL) == "VGLL3TOP.COLLAPSED" |
	names(MaleData_P2_OUL) == "ln.length_P1.scaled" | # just add length
	names(MaleData_P2_OUL) == "Tank" |
	names(MaleData_P2_OUL) == "DNA.ID.at.Tagging" )]
	
# predict equally spaced along length range > PER VGLL3TOP.COLLAPSED <
newdata$ln.length_P1.scaled = NA
newdata$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EE"] = seq(min(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EE"]), max(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EE"]), length= nrow(newdata[newdata$VGLL3TOP.COLLAPSED == "EE",]))
newdata$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EL/LE"] = seq(min(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EL/LE"]), max(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "EL/LE"]), length= nrow(newdata[newdata$VGLL3TOP.COLLAPSED == "EL/LE",]))
newdata$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "LL"] = seq(min(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "LL"]), max(MaleData_P2_OUL$ln.length_P1.scaled[newdata$VGLL3TOP.COLLAPSED == "LL"]), length= nrow(newdata[newdata$VGLL3TOP.COLLAPSED == "LL",]))

# check whether ranges are similar here:
range(MaleData_P2_OUL$ln.length_P1.scaled) 
range(newdata$ln.length_P1.scaled) 

library("MCMCglmm")
predict.model = predict.MCMCglmm(
		newdata = newdata,
		object= model, 
		marginal= model$Random$formula, 
		type= "response", 
		interval= "confidence")
# yields mean lower upper

predict.model = unique(data.frame(round(predict.model,4), newdata))
names(predict.model)[1:3] = c("mean","lower","upper")
head(predict.model)
meanMatLength = unique(predict.model[,c("mean","lower","upper","VGLL3TOP.COLLAPSED","ln.length_P1.scaled")])

# reverse scaling and log scale (getting length in cm)
RevScalePar = list(
	attr(scale(log(MaleData_P2_OUL$length.cm_P1)),"scaled:center"),
	attr(scale(log(MaleData_P2_OUL$length.cm_P1)),"scaled:scale"))
meanMatLength$"length.cm_P1" = exp( (meanMatLength$"ln.length_P1.scaled"* RevScalePar[[2]] + RevScalePar[[1]]) )

# change probit contrast to direction EE-LL
data.frame(VGLL3.Contrast = c("EE-LL"),
 contrast = mean(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"])),
 	lower = HPDinterval(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"]))[1],
 	upper =	HPDinterval(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"]))[2])

data.frame(VGLL3.Contrast = c("EE-EL"),
 contrast = mean(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"])),
 	lower = HPDinterval(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]))[1],
 	upper =	HPDinterval(model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] - (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]))[2])
```

Compare maturation probabilities between genotypes (others than EE constrasts). Output is 95% credible intervals - if the interval includes zero, the models are not different at *P*<0.05.
```{r}
# EL - LL
model = (model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDEL/LE"]  - 
		model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"(Intercept)"] + model.MCMCglmm_OUL.FullRandom.HetsCollapsed$Sol[,"VGLL3TOP.COLLAPSEDLL"])
# confidence interval
round(cbind(mean = mean(model),HPDinterval(model)),2)
# approximate p-value
2*length(which(model>0))/length(model)
```

Plot the predicted maturation probabilities for given length.
```{r}
ggplot(meanMatLength,aes(length.cm_P1,mean,color=VGLL3TOP.COLLAPSED)) +
  #geom_ribbon(aes(ymin=lower.MAT, ymax=upper.MAT), alpha=0.2) + 
  geom_line(size=2) + 
  theme_classic() + 
  scale_color_manual(values=c("blue","darkgrey", "red")) +
  xlab("Length at 8 months post-hatch (cm)") + 
  ylab("Maturation probability") +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black"),
        plot.title = element_text(size=14, face="bold",hjust=0.5),
        legend.position = "none") 
```

## Figures

Figure 1 with average maturation probability (statistically more solid than comparing maturation frequencies).
```{r}
# comparison p-values
my_comparisons1  = data.frame(start=c("EE"), 
                            end=c("LL"), 
                            y=c(0.4,0.4),
                            label=c("P=0.006"))

my_comparisons2  = data.frame(start=c("EE"), 
                            end=c("EL/LE"), 
                            y=c(0.36,0.36),
                            label=c("NS"))

my_comparisons3  = data.frame(start=c("EL/LE"), 
                            end=c("LL"),
                            y=c(0.31,0.31),
                            label=c("P=0.018"))

# add number of inds
meanMat$N = c(length(which(MaleData_P2_OUL$VGLL3TOP.COLLAPSED == "EE")),
              length(which(MaleData_P2_OUL$VGLL3TOP.COLLAPSED == "EL")),
              length(which(MaleData_P2_OUL$VGLL3TOP.COLLAPSED == "LL")))

probBars = ggplot(meanMat) + 
  geom_point(aes(x=VGLL3TOP.COLLAPSED,y=mean,colour=VGLL3TOP.COLLAPSED),size=5) + 
  geom_errorbar(aes(x=VGLL3TOP.COLLAPSED,y=mean,ymin=lower, ymax=upper,colour=VGLL3TOP.COLLAPSED),width=0) +
  scale_colour_manual(values=c("blue","darkgrey","red")) +
  ylab("Mean maturation probability") + 
  xlab("") +
  theme_classic() +
  geom_signif(data=my_comparisons1,
              aes(xmin=start, xmax=end, annotations=label, y_position=y),
              textsize = 3, size=0.5, tip_length = 0.01,
              manual=TRUE) +
  geom_signif(data=my_comparisons2,
              aes(xmin=start, xmax=end, annotations=label, y_position=y),
              textsize = 3, size=0.5, tip_length = 0.01,
              manual=TRUE) +
  geom_signif(data=my_comparisons3,
              aes(xmin=start, xmax=end, annotations=label, y_position=y),
              textsize = 3, size=0.5, tip_length = 0.01,
              manual=TRUE) +
  #geom_text(aes(y=-0.03,label=paste("N=",N,sep='')), vjust=-0.25) +
  scale_x_discrete(labels=c("EE","EL","LL")) +
  theme(axis.text.x=element_text(size=14,color="black",face="italic"),
        axis.title=element_text(size=14,color="black"),
        axis.text.y=element_text(size=14,color="black"),
        legend.position = "none") 

probModel = ggplot(meanMatLength,aes(length.cm_P1,mean,color=VGLL3TOP.COLLAPSED)) +
 # geom_ribbon(aes(ymin=lower.MAT, ymax=upper.MAT), alpha=0.1) + 
  geom_line(size=2) + 
  theme_classic() + 
  scale_color_manual(values=c("blue","darkgrey", "red")) +
  xlab("Length at 8 months of age (cm)") + 
  ylab("Maturation probability") +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black"),
        plot.title = element_text(size=14, face="bold",hjust=0.5),
        legend.position = "none") 

ggarrange(probBars, probModel, 
          labels = c("A","B"),
          ncol = 2, nrow = 1,
          widths = c(2, 3))
```

