twopiece R package with applications

The twopiece R package implements the family of Two-Piece distributions. The family of univariate two–piece distributions is a familly of univariate three-parameter location-scale models, where skewness is introduced by differing scale parameters either side of the location. The twopiece R package contains the implementation of Density, distribution function, quantile function and random generation for the 3 and 4 parameter two piece distribution with 3 parameterisations. Below, we present two examples in the context of data fitting and binary regression with a flexible link function. The package is avalailbe on R-forge, with extensions by David Rossell.

References.

  1. Inference in Two-Piece Location-Scale Models with Jeffreys Priors

Example 1: data fitting

rm(list=ls())
# Loading the twopiece package
library(twopiece)
## Loading required package: flexclust
## Loading required package: grid
## Loading required package: lattice
## Loading required package: modeltools
## Loading required package: stats4
## Loading required package: mvtnorm
## Loading required package: label.switching
# Data: BMI of 100 female Australian athletes
data <- c(20.56,20.67,21.86,21.88,18.96,21.04,21.69,20.62,22.64,19.44,25.75,21.20,22.03,25.44,22.63,21.86,22.27,21.27,23.47,23.19,
23.17,24.54,22.96,19.76,23.36,22.67,24.24,24.21,20.46,20.81,20.17,23.06,24.40,23.97,22.62,19.16,21.15,21.40,21.03,21.77,
21.38,21.47,24.45,22.63,22.80,23.58,20.06,23.01,24.64,18.26,24.47,23.99,26.24,20.04,25.72,25.64,19.87,23.35,22.42,20.42,
22.13,25.17,23.72,21.28,20.87,19.00,22.04,20.12,21.35,28.57,26.95,28.13,26.85,25.27,31.93,16.75,19.54,20.42,22.76,20.12,
22.35,19.16,20.77,19.37,22.37,17.54,19.06,20.30,20.15,25.36,22.12,21.25,20.53,17.06,18.29,18.37,18.93,17.79,17.05,20.31)
summary(data)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.75   20.27   21.82   21.99   23.39   31.93
n <- length(data)
# log-likelihood par=(mu,sigma,gamma,delta)
ll <- function(par){
if(par[2]>0 & par[3]>-1 & par[3]<1 & par[4]>0){
tempval <- dtp4(data,par[1],par[2],par[3],par[4],dt,param="eps",log=T)
return( -sum(tempval) )
}
else return(Inf)
}
# Optimization step
OPT <- optim(c(20,3,-0.25,10),ll)
# AIC, BIC, and Maximum Likelihood Estimators
AIC <- 2*OPT$value + 2*4
BIC <- 2*OPT$value + log(n)*4
MLE <- OPT$par
# Plot of the fitted density
tempf <- function(x) dtp4(x,MLE[1],MLE[2],MLE[3],MLE[4],dt,param="eps")
hist(data, breaks = 20 ,probability=T,xlim=c(15,35))
curve(tempf,15,35,add=T,lwd=2)
box()

Example 2: binary regression

rm(list=ls())
# Loading the twopiece package
library(twopiece)
# Flour beetle data
ni = c(59, 60, 62, 56, 63, 59, 62, 60) # number of beetles exposed
nk = c(6, 13, 18, 28, 52, 53, 61, 60) # number of beetles killed
points = nk/ni
dose = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839) # dose
n = length(ni)
# log-likelihood par=(alpha,beta,gamma)
ll = function(par){
if( par[3]>-1 & par[3]<1 ){
var = nk*ptp3(par[1]+par[2]*dose,0,1,par[3],pnorm,param="eps",log=T) + (ni-nk)*log(1-ptp3(par[1]+par[2]*dose,0,1,par[3],pnorm,param="eps",log=F))
return(-sum(var))
}
else return(-Inf)
}
# Optimization step
OPT <- optim(c(-40,20,0),ll)
# AIC, BIC, and Maximum Likelihood Estimators
AIC <- 2*OPT$value + 2*3
BIC <- 2*OPT$value + log(n)*3
MLE <- OPT$par
# Fitted dose-response curve using the two-piece normal link
DR <- function(d) ptp3(MLE[1]+MLE[2]*d,0,1,MLE[3],pnorm,param="eps",log=F)
plot(dose,points,xlim=c(1.65,1.95),ylim=c(0,1),xlab="Dose",ylab="Response",main="Bliss data")
curve(DR,1.65,1.95,add=T,lwd=2)