dsfa (version 2.0.2)

comper_mv: comper

Description

The comper implements the multivariate composed-error distribution in which the \(\mu_1\), \(\sigma_{V1}\), \(\sigma_{U2}\), \(\mu_2\), \(\sigma_{V2}\), \(\sigma_{U2}\) and \(\delta\) can depend on additive predictors. Useable with mgcv::gam, the additive predictors are specified via a list of formulae.

Usage

comper_mv(
  link = list("identity", "logshift", "logshift", "identity", "logshift", "logshift",
    "glogit"),
  s = c(-1, -1),
  distr = c("normhnorm", "normhnorm", "normal"),
  rot = 0,
  b = 0.01
)

Value

An object inheriting from class general.family of the mgcv package, which can be used in the mgcv and dsfa package.

Arguments

link

seven item list, specifying the links for \(\mu_1\), \(\sigma_{V1}\), \(\sigma_{U2}\), \(\mu_2\), \(\sigma_{V2}\), \(\sigma_{U2}\) and \(\delta\). See details.

s

integer vector of length two; each element corresponds to one marginal.

distr

string vector of length three; the first two elements determine the distribution of the marginals. Available are:
`normhnorm`, Normal-halfnormal distribution
`normexp`, Normal-exponential distribution
The last element determines the distribution of the copula:
`independent`, Independence copula
`normal`, Gaussian copula
`clayton`, Clayton copula
`gumbel`, Gumbel copula
`frank`, Frank copula
`joe`, Joe copula
`amh`, Ali-Mikhail-Haq copula

rot

integer determining the rotation for Archimedian copulas. Can be 90, 180 or 270.

b

positive parameter of the logshift link function.

Details

Used with gam() to fit distributional stochastic frontier model. The function is called with a list containing three formulae:

1. The first formula specifies the response of marginal one on the left hand side and the structure of the additive predictor for \(\mu_1\) parameter on the right hand side. Link function is "identity".

2. The second formula is one sided, specifying the additive predictor for the \(\sigma_{V1}\) on the right hand side. Link function is "logshift", e.g. \(\log \{ \sigma_{V1} \} + b \).

3. The third formula is one sided, specifying the additive predictor for the \(\sigma_{U1}\) on the right hand side. Link function is "logshift", e.g. \(\log \{ \sigma_{U1} \} + b \).

4. The fourth formula specifies the response of marginal two on the left hand side and the structure of the additive predictor for \(\mu_2\) parameter on the right hand side. Link function is "identity".

5. The fifth formula is one sided, specifying the additive predictor for the \(\sigma_{V2}\) on the right hand side. Link function is "logshift", e.g. \(\log \{ \sigma_{V2} \} + b \).

6. The sixth formula is one sided, specifying the additive predictor for the \(\sigma_{U2}\) on the right hand side. Link function is "logshift", e.g. \(\log \{ \sigma_{U2} \} + b \).

7. The seventh formula is one sided, specifying the additive predictor for the \(\delta\) on the right hand side. Link function is "glogit".

The fitted values and linear predictors for this family will be seven column matrices. For more details of the distribution see dcomper().

References

• schmidt2023multivariatedsfa

• wood2017generalizeddsfa

• aigner1977formulationdsfa

• kumbhakar2015practitionerdsfa

• azzalini2013skewdsfa

• schmidt2020analyticdsfa

Examples

# \donttest{
#Set seed, sample size and type of function
set.seed(1337)
N=1000 #Sample size
s<-c(-1,-1) #Set to production function for margin 1 and set to cost function for margin 2

distr_cop="normal"
distr_marg1="normhnorm"
distr_marg2="normhnorm"

#Generate covariates
x1<-runif(N,-1,1); x2<-runif(N,-1,1); x3<-runif(N,-1,1)
x4<-runif(N,-1,1); x5<-runif(N,-1,1); x6<-runif(N,-1,1)
x7<-runif(N,-1,1)

mu1=6+2*x1+(-2/3)*x1^2 #production function parameter 1
sigma_v1=exp(-1.5+sin(2*pi*x2)) #noise parameter 1
sigma_u1=exp(-1) #inefficiency parameter 1
mu2=5*x4^2+4*log(x4+2)^(1/4) #cost function parameter 2
sigma_v2=exp(-1.5) #noise parameter 2
sigma_u2=exp(-1+sin(2*pi*x6)) #inefficiency parameter 2
delta=transform(x=matrix(1+2.5*cos(4*x7)),
      type="glogitinv",
      par=delta_bounds(distr_cop), deriv_order = 0)

#Simulate responses and create dataset
Y<-rcomper_mv(n=N, mu=cbind(mu1,mu2),
              sigma_v=cbind(sigma_v1, sigma_v2),
              sigma_u = cbind(sigma_u1, sigma_u2), s=s,
              delta=delta,
              distr = c(distr_marg1,distr_marg2,distr_cop))
dat<-data.frame(y1=Y[,1],y2=Y[,2], x1, x2, x3, x4, x5, x6, x7)

#Write formulae for parameters
mu_1_formula<-y1~s(x1,bs="ps")
sigma_v1_formula<-~s(x2,bs="ps")
sigma_u1_formula<-~1
mu_2_formula<-y2~s(x4,bs="ps")
sigma_v2_formula<-~1
sigma_u2_formula<-~s(x6,bs="ps")
delta_formula<-~s(x7,bs="ps")

#Fit model
model<-dsfa(formula=list(mu_1_formula,sigma_v1_formula,sigma_u1_formula,
                              mu_2_formula,sigma_v2_formula,sigma_u2_formula,
                              delta_formula),  data=dat,
                 family=comper_mv(s=s, distr=c(distr_marg1,distr_marg2,distr_cop)),
                 optimizer="efs")

#Model summary
summary(model)

#Smooth effects
#Effect of x1 on the predictor of the production function of margin 1
plot(model, select=1) #Estimated function
lines(x1[order(x1)], 2*x1[order(x1)]+(-1/3)*x1[order(x1)]^2-
        mean(2*x1+(-1/3)*x1^2), col=2) #True effect

#Effect of x2 on the predictor of the noise of margin 1
plot(model, select=2) #Estimated function
lines(x2[order(x2)], -1.5+sin(2*pi*x2[order(x2)])-
        mean(-1.5+sin(2*pi*x2)),col=2) #True effect

#Effect of x4 on the predictor of the production function of margin 2
plot(model, select=3) #Estimated function
lines(x4[order(x4)], 3+5*x4[order(x4)]^2+4*log(x4[order(x4)]+2)^(1/4)-
        mean(3+5*x4^2+4*log(x4+2)^(1/4)), col=2) #True effect

#Effect of x6 on the predictor of the inefficiency of margin 2
plot(model, select=4) #Estimated function
lines(x6[order(x6)], -1+sin(2*pi*x6[order(x6)])-
        mean(-1+sin(2*pi*x6)),col=2) #True effect

#Effect of x7 on the predictor of the copula
plot(model, select=5) #Estimated function
lines(x7[order(x7)], 2.5*cos(4*x7[order(x7)])-
        mean(2.5*cos(4*x7)),col=2) #True effect

efficiency(model)
elasticity(model)

#' ### Second example with real data

data(BurkinaFarms)
data(BurkinaFarms_polys)

#Write formulae for parameters
mu_1_formula<-qharv_millet~s(land_millet, bs="ps")+s(labour_millet, bs="ps")+
                           s(material, bs="ps")+s(fert_millet, bs="ps")+
                           s(adm1, bs="mrf",xt=BurkinaFarms_polys)
sigma_v1_formula<-~1
sigma_u1_formula<-~farmtype+s(pest_millet, bs="ps")

mu_2_formula<-qharv_sorghum~s(land_sorghum, bs="ps")+s(labour_sorghum, bs="ps")+
                            s(material, bs="ps")+s(fert_sorghum, bs="ps")+
                            s(adm1, bs="mrf",xt=BurkinaFarms_polys)
sigma_v2_formula<-~1
sigma_u2_formula<-~farmtype+s(pest_sorghum, bs="ps")

delta_formula<-~1

model<-dsfa(formula=list(mu_1_formula, sigma_v1_formula, sigma_u1_formula,
                                mu_2_formula, sigma_v2_formula, sigma_u2_formula,
                                delta_formula),
                        data=BurkinaFarms,
                        family=comper_mv(s=c(-1,-1),
                        distr=c("normhnorm","normhnorm","normal")),
                        optimizer="efs")
plot(model) 
# }