# R/envelopes.R In santosneto/rbsmodels: Reparameterized Birnbaum-Saunders regression model

#### Documented in envelope

```#' Simulated Evelope
#' @description A normal plot with simulated envelope of the residual is produced.
#'
#'@usage envelope(model,k=19,alpha=0.05,res="deviance", precision = c("fixed","varying"),
#'
#' @param model object of class \code{gamlss} holding the fitted model.
#' @param k number of replications for envelope construction. Default is 19.
#' @param alpha value giving the size of the envelope. Default is 0.05 which is equivalent to a 95\% band.
#' @param res type of residuals to be extracted. Default is deviance.
#' @param precision If \code{precision = "fixed"} a model with precision fixed is used;
#' else a model with precision variable is used.
#' @param dist The function RBS() defines the RBS distribution.
#'
#'@return A simulated envelope of the class RBS.
#'
#' @author
#'Manoel Santos-Neto \email{[email protected]}, F.J.A. Cysneiros \email{[email protected]}, Victor Leiva \email{[email protected]} and Michelli Barros \email{[email protected]}
#'
#'@references
#'Atkinson, A. C. (1985) Plots, transformations and regression : an introduction to graphical methods of diagnostic regression analysis. Oxford Science Publications, Oxford.
#'
#'
#' @examples
#'
#'## Fixed Precision
#'library(faraway)
#'data(cpd)
#'attach(cpd)
#'model0 = gamlss(actual ~ projected, family=RBS(mu.link="identity"),method=CG())
#'summary(model0)
#'set.seed(2015)
#'envelope(model0)
#'model = gamlss(actual ~ 0+projected, family=RBS(mu.link="identity"),method=CG())
#'summary(model)
#'set.seed(2015)
#'envelope(model)
#'
#'##
#'library(alr3)
#'data(landrent)
#'attach(landrent)
#'resp <- I(Y/X1)
#'y1 <-  split(resp, X4)\$"1"
#'x21 <-  split(X2, X4)\$"1"
#'
#'##Fixed Precision
#'fit0 <- gamlss(y1 ~ x21, family=RBS(mu.link="identity"),method=CG()  )
#'summary(fit0)
#'set.seed(2015)
#'##Varying Precision
#'summary(fit11)
#'set.seed(2015)
#' @export

envelope <- function(model,k=19,alpha=0.05,res="deviance", precision = c("fixed","varying"), dist = RBS(mu.link = "identity",sigma.link = "identity") )
{

if(precision!= "fixed")
{
alfa1 <- ceiling(k*alpha)
alfa2 <- ceiling(k*(1-alpha))
n <- model\$N
td  <- residuals(model,residual=res)
sigma <- model\$sigma.fv
mu <- model\$mu.fv
re <- matrix(0,n,k)
X <- model\$mu.x
Z <- model\$sigma.x
p<- ncol(X)
q<- ncol(Z)

for(i in 1:k)
{

y1 <- mapply(rRBS,n=1,mu=mu,sigma=sigma)
nresp <- y1

if(p == 1) form <- nresp ~ 0 + X
else form <- nresp ~ X[,-1]

if(q == 1) form1 <- nresp ~ 0 + Z
else  form1 <- nresp ~ Z[,-1]

conh0 = gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
model1 <- gamlss(formula=form,sigma.formula = form1 ,family=dist,method=CG(),control = conh0)
rdf <- residuals(model1,residual=res)

re[,i]=sort(rdf)
}
e1 = numeric(n)
e2 = numeric(n)
for(l in 1:n)
{
eo = sort(re[l,])
e1[l] = eo[alfa1]
e2[l] = eo[alfa2]
}
a<-  qqnorm(e1,plot.it=FALSE)\$x
a1<-  qqnorm(e1,plot.it=FALSE)\$y
b<-  qqnorm(e2,plot.it=FALSE)\$x
b1<-  qqnorm(e2,plot.it=FALSE)\$y
r<-  qqnorm(td,plot.it=FALSE)\$x
r1<-  qqnorm(td,plot.it=FALSE)\$y
xx <- c(a,rev(b))
yy <- c(a1,rev(b1))

xb = apply(re,1,mean)
faixa = range(td,e1,e2)
par(mar=c(4., 4., 0.1, 0.1))
plot(r,r1,type = "n", ylim=faixa,axes=FALSE,xlab="",ylab="")
par(new=TRUE)
polygon(xx,yy,col="gray80",border=NA)
par(new=TRUE)
qqnorm(e1,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=1,lwd=1,col="gray")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=1,lwd=1,col="gray")
par(new=TRUE)
qqnorm(xb,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=2,lwd=1,col="black")
par(new=TRUE)
qqnorm(td,xlab="Theorical Quantile",main="",ylab="Empirical Quantile",ylim=faixa,pch=20,cex=1,lwd=1)
} else{
alfa1 <- ceiling(k*alpha)
alfa2 <- ceiling(k*(1-alpha))
n<-model\$N
td  <- residuals(model, residual= res)
sigma <- model\$sigma.fv[1]
mu <- model\$mu.fv
alpha<-sqrt(2/sigma)
bet.a<- (sigma*mu)/(sigma+1)
re <- matrix(0,n,k)
X <- model\$mu.x
p <- ncol(X)

for(i in 1:k)
{

y1 <- mapply(rRBS,n=1,mu=mu,sigma=sigma)
nresp <- y1
if(p == 1) form <- nresp ~ 0 + X
else  form <- nresp ~ X[,-1];

conh0 <- gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
model1 <- gamlss(formula=form ,family=dist,method=CG(),control = conh0)
rdf <- residuals(model1,residual=res)
re[,i] <- sort(rdf)
}
e1 = numeric(n)
e2 = numeric(n)
for(l in 1:n)
{
eo = sort(re[l,])
e1[l] = eo[alfa1]
e2[l] = eo[alfa2]
}

a<-  qqnorm(e1,plot.it=FALSE)\$x
a1<-  qqnorm(e1,plot.it=FALSE)\$y
b<-  qqnorm(e2,plot.it=FALSE)\$x
b1<-  qqnorm(e2,plot.it=FALSE)\$y
r<-  qqnorm(td,plot.it=FALSE)\$x
r1<-  qqnorm(td,plot.it=FALSE)\$y
xx <- c(a,rev(b))
yy <- c(a1,rev(b1))

xb = apply(re,1,mean)
faixa = range(td,e1,e2)
par(mar=c(4., 4., 0.1, 0.1))
plot(r,r1,type = "n", ylim=faixa,axes=FALSE,xlab="",ylab="")
par(new=TRUE)
polygon(xx,yy,col="gray80",border=NA)
par(new=T)
qqnorm(e1,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=1,lwd=1,col="gray")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=1,lwd=1,col="gray")
par(new=TRUE)
qqnorm(xb,axes=F,xlab="",ylab="",main="",type="l",ylim=faixa,lty=2,lwd=1,col="black")
par(new=TRUE)
qqnorm(td,xlab="Theorical Quantile",main="",ylab="Empirical Quantile",ylim=faixa,pch=20,cex=1,lwd=1)
}
}
```
santosneto/rbsmodels documentation built on May 26, 2017, 12:32 a.m.