envelope.inverse.gaussian<-function(form=form,Fam=Fam,k=k,alfa=alfa){
alfa1<-ceiling(k*alfa)
alfa2<-ceiling(k*(1-alfa))
glm1<-glm(formula=form,family=Fam,maxit=50)
X<-model.matrix(glm1)
w<-fit$weights
W<-diag(w)
H<-sqrt(W)%*%X%*%solve(t(X)%*%W%*%X)%*%t(X)%*%sqrt(W)
h<-diag(H)
n<-nrow(X)
rd<-resid(glm1,type="deviance")  
phi<-n/(deviance(glm1))
td<-rd*sqrt(phi/(1-h))
re<-matrix(0,n,k)
for(i in 1:k){
nresp<-rig(n,fitted(glm1),l=1)
fit<-glm(nresp~-1+X,family=Fam,maxit=50)
w<-fit$weights
W<-diag(w)
H<-sqrt(W)%*%X%*%solve(t(X)%*%W%*%X)%*%t(X)%*%sqrt(W)
h<-diag(H)
re[,i]<-sort(resid(fit,type="deviance")*sqrt(phi/(1-h)))
             }
e1<-numeric(n)
e2<-numeric(n)
for(i in 1:n){
eo<-sort(re[i,])
e1[i]<-eo[alfa1]
e2[i]<-eo[alfa2]
             }
xb<-apply(re,1,mean)
faixa<-range(td,e1,e2)
par(pty="s")
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2)
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2)
par(new=T)
qqnorm(xb,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(td,xlab="Quantis da N(0,1)",ylab="Residuo Componente do desvio",ylim=faixa)
}