ssym (version 1.3)

ssym.l: Fitting Semiparametric Symmetric (or Log-Symmetric) Regression Models

Description

ssym.l is used to fit a semiparametric regression model suitable for data set analysis in which the conditional distribution of the response (or log-response) is symmetric and continuous. In this setup, both location and dispersion parameters of the response (or log-response) variable distribution are explicitly modeled through semiparametric functions.

Usage

ssym.l(response, formula.mu, ncs.mu, start.lambda.mu, lambda.mu, formula.phi, ncs.phi,
start.lambda.phi, lambda.phi, family, xi, epsilon, maxiter, subset, local.influence)

Arguments

response

the response (or log-response) variable.

formula.mu

(optional) a symbolic description of the parametric function to be fitted to the location parameter.

ncs.mu

(optional) an explanatory (continuous) variable to be used in the nonparametric function to be fitted to the location parameter.

start.lambda.mu

(optional) a numeric value of starting estimate for the smoothing parameter associated to the location submodel.

lambda.mu

(optional) a numerical value for the smoothing parameter indicating that it is provided by the user rather than estimated from the data.

formula.phi

(optional) a symbolic description of the parametric function to be fitted to the dispersion parameter.

ncs.phi

(optional) an explanatory (continuous) variable to be used in the nonparametric function to be fitted to the dispersion parameter.

start.lambda.phi

(optional) a numeric value of starting estimate for the smoothing parameter associated to the location submodel.

lambda.phi

(optional) a numerical value for the smoothing parameter indicating that it is provided by the user rather than estimated from the data.

family

a description of the error distribution to be used in the model. Supported families include Normal, Student, Contnormal, Powerexp, Hyperbolic, Slash, Sinh-normal and Sinh-t

xi

a numeric value or numeric vector that represents the extra parameter of the specified error distribution.

epsilon

(optional) positive convergence tolerance. Default value is 0.0000001.

maxiter

(optional) an positive integer giving the maximal number of iterations for the estimating process. Default value is 1000.

subset

(optional) expression indicating individuals or observations to keep (or drop).

local.influence

(optional) logical. If TRUE, local influence measures for the location parameters vector under two perturbation schemes are calculated.

Value

• coefs.mu
  a vector of parameter estimates associated with the nonlinear function fitted to the location of the response (or log-response) variable.
• coefs.phi
  a vector of parameter estimates associated with the semiparametric function fitted to the dispersion of the response (or log-response) variable.
• se.mu
  approximated standard errors asociated with coefs.mu.
• se.phi
  approximated standard errors asociated with coefs.phi.
• v
  final weights of the iterative process.
• lambda.mu
  estimate for the smoothing parameter associated to the location parameter.
• gle.mu
  degrees of freedom associated with the nonparametric part of the location submodel.
• lambda.phi
  estimate for the smoothing parameter associated to the dispersion parameter.
• gle.phi
  degrees of freedom associated with the nonparametric part of the dispersion submodel.
• deviance.mu
  a vector of deviances associated with the location of the response (or log-response) variable.
• deviance.phi
  a vector of deviances associated with the dispersion of the response (or log-response) variable.
• mu.fitted
  a vector of fitted values for the location of the response (or log-response) variable.
• phi.fitted
  a vector of fitted values for the dispersion of the response (or log-response) variable.
• lpdf
  a vector of individual contributions to the log-likelihood function.
• pdfz
  a vector of the cumulative distribution function of each individual.
• cw
  if local.influence=TRUE, a matrix of local influence and total local influence measures (under the case-weight perturbation scheme) asociated with coefs.mu and based on the conformal normal curvature.
• pr
  if local.influence=TRUE, a matrix of local influence and total local influence measures (under the response perturbation scheme) asociated with coefs.mu and based on the conformal normal curvature.

Details

The iterative estimation process for the parameters of interest is based on the Fisher scoring and backfitting algorithms. Because some distributions such as Student-t, contaminated normal, power exponential, slash and symmetric hyperbolic may be obtained as a scale mixture of normal distributions, the EM algorithm is applied in those cases to obtain a more efficient iterative process for the parameter estimation. Further, because the Sinh-t distribution can be obtained as a scale mixture of Sinh-normal distributions, the EM algorithm is also applied in that case to obtain a more efficient iterative process for the parameter estimation. The smoothing parameter is chosen using the cross-validation score.

References

Vanegas, L.H. and Paula, G.A. (2014) A Semiparametric Approach for Joint Modeling of Median and Skewness. TEST (submitted) Vanegas, L.H. and Paula, G.A. (2014) Log-symmetric distributions: statistical properties and parameter estimation. Brazilian Journal of Probability and Statistics (submitted)

See Also

ssym.nl

Examples

#########################################################################################
############## Texture of Snacks Data - a log-Student-t semiparametric model ############
#########################################################################################

data(Snacks)
Snacks2 <- Snacks[order(Snacks$snack,Snacks$week,Snacks$texture),]
attach(Snacks2)

snack <- factor(snack)

fit <- ssym.l(log(texture), formula.mu=~snack, ncs.mu=week, formula.phi=~snack,
              family='Student', xi=14, local.influence=TRUE)
summary(fit)

####################### Plot of the fitted model ##############################

par(mfrow=c(1,2))
xl <- "Week"
rx <- range(week)
ry <- range(texture)
h <- fit$coefs.mu[5:length(fit$coefs.mu)]
sa <- ncs.graph(week,h,1000)

plot(week,texture,xlim=rx,ylim=ry,type="p",cex=0.3,lwd=3,ylab="",xlab="")
par(new=TRUE)
plot(sa[,1],exp(sa[,2]),xlim=rx,ylim=ry,type="l",ylab="",xlab=xl,lty=1,main="")
for(i in 2:5){
	 par(new=TRUE)
	 plot(sa[,1],exp(sa[,2] + fit$coefs.mu[i-1]),xlim=rx,ylim=ry,type="l",
	      ylab="",xlab="",lty=i,main="Median")
}
legend(rx[1],ry[2],lty=1:5,bty="n",legend=paste("snack",1:5))

r <- (log(texture) - fit$mu.fitted)^2/fit$xix
ry <- range(r)
plot(week,r,xlim=rx,ylim=ry,type="p",cex=0.3,lwd=3,ylab="",xlab=xl,main="Skewness")
for(i in 1:5){
	abline(h=mean(fit$phi.fitted[snack==i]), lty=i)
}
legend(rx[1],ry[2],lty=1:5,bty="n",legend=paste("snack",1:5))

########################### Residual analysis ##################################

par(mfrow=c(1,2))
m1 <- "Residuals for the median submodel"
res.dev.mu <- sqrt(fit$deviance.mu)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.mu,-3.5),max(res.dev.mu,3.5))
plot(week,res.dev.mu,cex=0.3,lwd=3,ylim=ry,main=m1,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)
m2 <- "Residuals for the skewness submodel"
res.dev.phi <- sqrt(fit$deviance.phi)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.phi,-3.5),max(res.dev.phi,3.5))
plot(week,res.dev.phi,cex=0.3,lwd=3,ylim=ry,main=m2,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)

########################### Sensitivity analysis ##################################

par(mfrow=c(1,2))
m1 <- "Local Influence under case-weight perturbation scheme"
m2 <- "Total Local Influence under case-weight perturbation scheme"
plot(fit$cw[,1], type="h", main=m1, xlab="Index", ylab="")
plot(fit$cw[,2], type="h", main=m2, xlab="Index", ylab="")

###################################################################################
########### Fraction of Cell Volume Data - a Power-exponential model  ###########
###################################################################################

data(Ovocytes)
attach(Ovocytes)
type <- factor(type)

fit <- ssym.l(fraction, formula.mu=~type, ncs.mu=time, formula.phi=~type,
              family="Powerexp", xi=-0.22, local.influence=TRUE)
summary(fit)

####################### Plot of the fitted model ##############################

par(mfrow=c(1,2))
xl <- "Time"
rx <- range(time)
ry <- range(fraction)
h <- fit$coefs.mu[2:length(fit$coefs.mu)]
sa <- ncs.graph(time,h,1000)

plot(time[type=='Mature'],fraction[type=='Mature'],xlim=rx,ylim=ry,type="p",
     cex=0.5,lwd=1,ylab="",xlab="")
par(new=TRUE)
plot(time[type=='Immature'],fraction[type=='Immature'],xlim=rx,ylim=ry,
     type="p",cex=0.5,lwd=2,ylab="",xlab="")
par(new=TRUE)
plot(sa[,1],(sa[,2]),xlim=rx,ylim=ry,type="l",ylab="",xlab=xl,lty=1,main="")
par(new=TRUE)
plot(sa[,1],(sa[,2] + fit$coefs.mu[1]),xlim=rx,ylim=ry,type="l",ylab="",xlab=xl,
     lty=1,main="Mean")
legend(rx[1],ry[2],pt.lwd=c(1,2),bty="n",legend=c("Mature","Immature"),
       pt.cex=0.5,pch=1)

r <- (fraction - fit$mu.fitted)^2/fit$xix
ry <- range(r)
plot(time[type=='Mature'],r[type=='Mature'],xlim=rx,ylim=ry,type="p",cex=0.5,
     lwd=1,ylab="",xlab="")
par(new=TRUE)
plot(time[type=='Immature'],r[type=='Immature'],xlim=rx,ylim=ry,type="p",
     cex=0.5,lwd=2,ylab="",xlab="",main="Dispersion")
abline(h=exp(fit$coefs.phi[1]))
abline(h=exp(fit$coefs.phi[1] + fit$coefs.phi[2]))
legend(rx[1],ry[2],pt.lwd=c(1,2),bty="n",legend=c("Mature","Immature"),
       pt.cex=0.5,pch=1)

########################### Residual analysis ##################################

par(mfrow=c(1,2))
m1 <- "Residuals for the mean submodel"
res.dev.mu <- sqrt(fit$deviance.mu)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.mu,-3.5),max(res.dev.mu,3.5))
plot(time,res.dev.mu,cex=0.3,lwd=3,ylim=ry,main=m1,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)
m2 <- "Residuals for the dispersion submodel"
res.dev.phi <- sqrt(fit$deviance.phi)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.phi,-3.5),max(res.dev.phi,3.5))
plot(time,res.dev.phi,cex=0.3,lwd=3,ylim=ry,main=m2,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)

########################### Sensitivity analysis ##################################

par(mfrow=c(1,2))
m1 <- "Local Influence under case-weight perturbation scheme"
m2 <- "Total Local Influence under case-weight perturbation scheme"
plot(fit$cw[,1], type="h", main=m1, xlab="Index", ylab="")
plot(fit$cw[,2], type="h", main=m2, xlab="Index", ylab="")

#########################################################################################
################## Biaxial Fatigue Data - a Birnbaum-Saunders model   ##################
#########################################################################################

data(Biaxial)
Biaxial2 <- Biaxial[order(Biaxial$Work,Biaxial$Life),]
attach(Biaxial2)

fit <- ssym.l(log(Life), formula.mu=~log(Work), family='Sinh-normal', xi=1.54,
              local.influence=TRUE)
summary(fit)

####################### Plot of the fitted model ##############################

xl <- "Work per cycle"
rx <-range(Work)
ry <- range(Life)
plot(Work,Life,xlim=rx,ylim=ry,type="p",cex=0.3,lwd=3,ylab="",xlab="")
par(new=TRUE)
plot(Work,exp(fit$mu.fitted),xlim=rx,ylim=ry,type="l",ylab="Life",xlab=xl,main="Median")

########################### Residual analysis ##################################

res.dev.mu <- sqrt(fit$deviance.mu)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.mu,-3.5),max(res.dev.mu,3.5))
plot(Work,res.dev.mu,cex=0.3,lwd=3,ylim=ry,main="Residuals",xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)

########################### Sensitivity analysis ##################################

par(mfrow=c(1,2))
m1 <- "Local Influence under case-weight perturbation scheme"
m2 <- "Total Local Influence under case-weight perturbation scheme"
plot(fit$cw[,1], type="h", main=m1, xlab="Index", ylab="")
plot(fit$cw[,2], type="h", main=m2, xlab="Index", ylab="")

###################################################################################
############## Personal Claims Data - a Birnbaum-Saunders-t model  ################
###################################################################################

data(Claims)
Claims2 <- Claims[order(Claims$op_time,Claims$total),]
attach(Claims2)

fit <- ssym.l(log(total), formula.mu=~op_time, formula.phi=~op_time, family='Sinh-t',
              xi=c(0.1,4), local.influence=TRUE)
summary(fit)

####################### Plot of the fitted model ##############################

yl <- "Amount of paid money"
xl <- "Operational time"
rx <- range(op_time)
ry <- range(total)
par(mfrow=c(1,2))
plot(op_time,total,xlim=rx,ylim=ry,type="p",cex=0.3,lwd=3,ylab="",xlab="")
par(new=TRUE)
plot(op_time,exp(fit$mu.fitted),xlim=rx,ylim=ry,type="l",ylab=yl,xlab=xl,
     main="Median")

r <- (log(total) - fit$mu.fitted)^2/fit$xix
ry <- range(r)
plot(op_time,r,xlim=rx,ylim=ry,type="p",cex=0.3,lwd=3,ylab="",xlab="")
par(new=TRUE)
plot(op_time,fit$phi.fitted,xlim=rx,ylim=ry,type="l",ylab="",xlab=xl,
     main="Skewness")

########################### Residual analysis ##################################

par(mfrow=c(1,2))
m1 <- "Residuals for the median submodel"
res.dev.mu <- sqrt(fit$deviance.mu)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.mu,-3.5),max(res.dev.mu,3.5))
plot(op_time,res.dev.mu,cex=0.3,lwd=3,ylim=ry,main=m1,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)
m2 <- "Residuals for the skewness submodel"
res.dev.phi <- sqrt(fit$deviance.phi)*ifelse(fit$z.hat>=0,1,-1)
ry <- c(min(res.dev.phi,-3.5),max(res.dev.phi,3.5))
plot(op_time,res.dev.phi,cex=0.3,lwd=3,ylim=ry,main=m2,xlab=xl,ylab="")
abline(h=-3,lty=3)
abline(h=+3,lty=3)

########################### Sensitivity analysis ##################################

par(mfrow=c(1,2))
m1 <- "Local Influence under case-weight perturbation scheme"
m2 <- "Total Local Influence under case-weight perturbation scheme"
plot(fit$cw[,1], type="h", main=m1, xlab="Index", ylab="")
plot(fit$cw[,2], type="h", main=m2, xlab="Index", ylab="")