Example

We consider fitting various probability models to the following toxicity data:

13.265  8.273 21.223 16.236 17.021  3.285

The plots in Figure 1 show that only one of the 15 fitted distributions (3-parameter gamma) is rejected by goodness-of-fit tests.

Figure 1. Q-Q Plots for a variety of theoretical distributions

Let’s look at some of the estimated quantiles.

Clearly, some of these are poor estimates even though the probability model was deemed to be adequate.

Given we have no theoretical basis for selecting one distribution over another, it makes sense to ‘combine’ (= weight) the results where the weights reflect our degree of confidence in the model fit.

Outline

\[\begin{align} & M=\text{ number of candidate models} \\ & y=\text{ response variable} \\ & \theta \text{ = parameter of interest} \\ \end{align}\]

What we’re interested in determining is \(p\left( \theta \left| y \right. \right)=\text{ posterior probability for }\theta \text{ given data}\)

It is readily shown that \(p\left( \theta \left| y \right. \right)=\sum\limits_{m=1}^{M}{p\left( m\left| y \right. \right)\ }p\left( \theta \left| y \right.,m \right)\) where: \[\begin{align} & p\left( m\left| y \right. \right)=\text{ posterior probability for model }m\text{ given data} \\ & p\left( \theta \left| y \right.,m \right)=\text{ posterior probability for }\theta \text{ given model }m\text{ and data} \\ \end{align}\]

Furthermore, \(p\left( m\left| y \right. \right)\propto p\left( m \right)p\left( y\left| m \right. \right)\) where

\[\begin{align} & p\left( m \right)=\text{ prior probability for model }m\text{ } \\ & p\left( y\left| m \right. \right)=\text{marginal likelihodd for model }m\text{ } \\ \end{align}\]

We can simulate from the required posterior probability density as follows:

1. Select data-generating model at random from our set of models using prior model probabilities;
2. Generate parameter values using specified prior distributions;
3. Generate data from selected model with chosen parameter values;
4. Repeat many times.

A weighted average is simply:

\[\bar{\theta }\text{=}\sum\limits_{m=1}^{M}{{{w}_{m}}}{{{\hat{\theta }}}_{m}}\text{ with }0<{{w}_{m}}<1;\text{ and }\sum\limits_{m=1}^{M}{{{w}_{m}}}=1\]

Let’s use R to do this ‘manually’.

require(fitdistrplus)

## Loading required package: fitdistrplus

## Loading required package: MASS

## Loading required package: survival

dat<-c(13.265,  8.273, 21.223, 16.236, 17.021 , 3.285)
dist<-list("norm","lnorm","logis","gamma")
fit<-lapply(dist,function(x) fitdist(data=dat,x ))
aic<-data.frame(dist=unlist(dist),AIC=NA)
aic$AIC<-unlist(lapply(fit,function(x) x$aic))
aic$delta<-aic$AIC-min(aic$AIC)
aic$wt<-exp(-0.5*aic$delta)/sum(exp(-0.5*aic$delta))
aic$HC1<-unlist(lapply(dist,function(x) do.call(paste("q",x,sep=""),list(0.01,fit[[1]][[1]][1],fit[[1]][[1]][2]))))
aic$HC5<-unlist(lapply(dist,function(x) do.call(paste("q",x,sep=""),list(0.05,fit[[1]][[1]][1],fit[[1]][[1]][2]))))
aic$HC10<-unlist(lapply(dist,function(x) do.call(paste("q",x,sep=""),list(0.10,fit[[1]][[1]][1],fit[[1]][[1]][2]))))
aic[,-1]<-round(aic[,-1],3)
aic

##    dist    AIC delta    wt     HC1    HC5    HC10
## 1  norm 42.381 0.000 0.363  -0.571  3.468   5.622
## 2 lnorm 44.520 2.139 0.125   0.565 32.080 276.304
## 3 logis 42.809 0.428 0.293 -14.018 -4.234   0.194
## 4 gamma 43.395 1.014 0.219   1.054  1.326   1.489

So, our model-averaged SSD (MASSD) is a weighted combination of the normal, log-normal, logistic, and gamma distributions in the following proportions: 0.363, 0.125, 0.293, and 0.219 respectively.

What about the estimated HC_x values?

You might be tempted to think that the estimated HC_x can be obtained by simply applying the weights to the separate HC_x values obtained from each of the distributions used in forming the MASSD. THIS WOULD BE WRONG!

It’s a trivial task to demonstrate this. For example, the HC₅ from our model-averaged fit is 5.47 whereas the weighted HC₅ is 3.59 (see later). We can substitute these values back into our model-averaged SSD (calculations not shown) and determine the (theoretical) fraction affected:

## Model-averaged HC5 fraction affected = 5.00%

## Weighted HC5 fraction affected =  2.31%

Obtaining HC_x values from a model-averaged SSD

Recall that an HC_x is a quantile from the fitted probability model. Now, the SSD is in fact what statisticians refer to as a probability density function or pfd. We denote this as \({f_X}\left( x \right)\).

So, for example, the familiar normal pdf with \(mean = \mu\) and \(variance = \sigma^2\) is:\[{f_X}\left( x \right) = \frac{1}{{\sigma \sqrt {2\pi } }}{e^{ – \frac{1}{2}{{\left( {\frac{{x – \mu }}{\sigma }} \right)}^2}}}\]

The cumulative distribution function or cdf is a mathematical expression (if it exists) corresponding to the cumulative probability \(P\left[ {X \le r} \right]\) where r is any number and is written:\[P\left[ {X \le r} \right] = {F_X}\left( r \right) = \int_{ – \infty }^r {{f_X}\left( x \right)} \;dx\]

Now, the \({p^{th}}\) quantile (let’s call it \({\xi _p}\)) is such that:\[P\left[ {X \le {\xi _p}} \right] = p\]

So, for example, the median (or \({\xi _{0.5}}\)) is found by finding the value \({{\xi _{0.5}}}\) such that:\[P\left[ {X \le {\xi _{0.5}}} \right] = 0.5\].

From our definition of \({F_X}\left( \cdot \right)\) this means solving the integral equation:\[\int_{ – \infty }^{{\xi _{0.5}}} {{f_X}\left( x \right)} \;dx = 0.5\]

Fortunately, there’s an easy way in R – we use the inverse cdf or quantile functions.

So, for example, if we wanted to find the \({5^{th}}\) percentile of a normal distribution having mean=10 and standard deviation=12 :

qnorm(0.05,10,12)

## [1] -9.738244

or the \({80^{th}}\) percentile from the same distribution:

qnorm(0.8,10,12)

## [1] 20.09945

Now, returning to our model averaged SSD.

As explained above, the model averaged SSD is just a weighted combination of k individual pdf’s. That is:

\[\begin{array}{*{20}{c}} {{f_Y}\left( y \right) = \sum\limits_{i = 1}^k {{w_i}{f_i}\left( y \right)} }&{\begin{array}{*{20}{c}} {{\rm{with}}}&{\begin{array}{*{20}{c}} {0 \le {w_i} \le 1;}&{{\rm{and}}} \end{array}}&{\sum\limits_{i = 1}^k {{w_i} = 1} } \end{array}} \end{array}\]

The model-averaged cumulative distribution function is similarly defined:

\[\begin{array}{*{20}{c}} {{F_Y}\left( y \right) = \sum\limits_{i = 1}^k {{w_i}{F_i}\left( y \right)} }&{\begin{array}{*{20}{c}} {{\rm{with}}}&{\begin{array}{*{20}{c}} {0 \le {w_i} \le 1;}&{{\rm{and}}} \end{array}}&{\sum\limits_{i = 1}^k {{w_i} = 1} } \end{array}} \end{array}\]

So, the correct way of determining a model-averaged HC_x is to solve:\[{F_Y}\left( {H{C_x}} \right) = x\] for \({H{C_x}}\).

This is no longer a simple matter of invoking one of R's intrinsic functions – although it can still be done with a single line of R code.

This is done by finding the root of the equation \({F_Y}\left( {H{C_x}} \right) – x\). In other words, solving:

\[{F_Y}\left( {H{C_x}} \right) – x = 0\]

The R function which does this is called uniroot which has the following form:

uniroot(f, interval, ...,
        lower = min(interval), upper = max(interval),
        f.lower = f(lower, ...), f.upper = f(upper, ...),
        extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
        tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)

where

To see how this works, consider the previous example.

param<-lapply(fit, function (x) x$estimate)  # extract the estimated parameters for each component distribution

## Next form the model-averaged SSD

massd<-function(x,p){
  
F<- 0.363*pnorm(x,param[[1]][1],param[[1]][2])+0.125*plnorm(x,param[[2]][1],param[[2]][2])
       + 0.293*plogis(x,param[[3]][1],param[[3]][2])+0.219*pgamma(x,param[[4]][1],param[[4]][2])
return(F-p) #  p is probability (eg: 0.05 for an HC5) - to be supplied       
}

uniroot(massd,c(0,10),p=0.05)$root

## [1] 5.472128

p<-0.05
 0.363*qnorm(p,param[[1]][1],param[[1]][2])+0.125*qlnorm(p,param[[2]][1],param[[2]][2])+
        0.293*qlogis(p,param[[3]][1],param[[3]][2])+0.219*qgamma(p,param[[4]][1],param[[4]][2])

## [1] 3.590141