envelope.binomialcr<-function(form=form,Fam=Fam,k=k,alfa=alfa,total=total){ 
d<-length(total)
alfa1<-ceiling(k*alfa)
alfa2<-ceiling(k*(1-alfa)) 
glm1<-glm(formula=form,family=Fam,maxit=50) 
X<-model.matrix(glm1)
w<-glm1$weights
W<-diag(w)
H<-sqrt(W)%*%X%*%solve(t(X)%*%W%*%X)%*%t(X)%*%sqrt(W)
h<-diag(H)
rd<-resid(glm1,type="deviance")   
phi<-1
td<-rd*sqrt(phi/(1-h)) 
res<-matrix(0,d,k) 
y<-numeric(d) 
for(i in 1:k){ 
for(j in 1:d){ 
dif<-runif(total[j])-fitted(glm1)[j] 
dif[dif>=0]<-0 
dif[dif<0]<-1 
y[j]<-sum(dif)
             } 
xmat<-cbind(y,total-y) 
fit<-glm(xmat~-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)
res[,i]<-sort(resid(fit,type="deviance")*sqrt(phi/(1-h)))
             } 
e1<-numeric(d) 
e2<-numeric(d) 
for(i in 1:d){ 
eo<-sort(res[i,]) 
e1[i]<-eo[alfa1] 
e2[i]<-eo[alfa2]
             } 
xb<-apply(res,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)
}