glm.rest<-
function(x=X,y=Y,w = rep(1, length(x[, 1])), start = NULL, offset = 0,
        family= gaussian(), maxit=10, epsilon = 0.0001, trace = F, 
        null.dev = NULL, qr = F, ...)
{
pass<-1
X<-x
Y<-y
n<-nrow(X)
p<-ncol(X)
Fam<-as.family(family)
if(any(offset) && dimnames(X)[[2]][1]=="(Intercept)"&&p==1) 
     {
	deviance<-list(deviance=glm.fit(X[,"(Intercept)", drop = F], Y, w, 
                      offset = offset, family = family, maxit = maxit,trace=trace
			, epsilon = epsilon, null.dev
			= NULL)$deviance )
     } else 
{
fr<- glm.fit(x=X, y=Y, w = w, start = NULL,
             offset = offset, family= Fam, maxit = pass, epsilon = epsilon,
             trace = F, null.dev = T, qr = qr, ...)
we<-fr$weights
W<-diag(as.vector(we))
betai<-matrix(fr$coef,p,1)
const<-solve(t(X)%*%W%*%X)%*%t(Cres)%*%
solve(Cres%*%solve(t(X)%*%W%*%X)%*%t(Cres))%*%(sol-Cres%*%betai)
beta<-betai+const
k<-0
if(trace) cat("GLMREST    linear loop ", k, " \n", sep = "",
": coef ", format(round(as.vector(beta), 4))," \n", sep = "")
while((any(abs(beta-betai)>epsilon))&&k<maxit){
neta<-X%*%beta+offset
fr<- glm.fit(x=X, y=Y, w = w, start = neta, offset = offset,
             family= Fam, maxit = pass, epsilon = epsilon, trace = F, 
             null.dev = NULL,  qr = qr, ...)
betat<-matrix(fr$coef,p,1)
we<-fr$weights
W<-diag(as.vector(we))
betai<-beta
const<-solve(t(X)%*%W%*%X)%*%t(Cres)%*%
solve(Cres%*%solve(t(X)%*%W%*%X)%*%t(Cres))%*%(sol-Cres%*%betat)
beta<-betat+const
k<-k+1
if(trace) cat("GLMREST    linear loop ", k, " \n", sep = "",
": coef ", format(round(as.vector(beta), 4))," \n", sep = "")
}
if(maxit==k)
 warning(paste("linear convergence not obtained in ", k," iterations."))
neta<-X%*%beta+offset
mu<-Fam$inverse(neta)
df.residual<-fr$df.residual
if (Fam$family[[1]]=="Binomial") {
     if(is.matrix(Y)) {
		if(dim(Y)[2] > 2)
			stop("only binomial response matrices (2 columns)")
		n <- drop(Y%*% c(1, 1))
                y<-Y[,1]
              	      }
	else {
		if(is.category(Y))
			y <- Y != levels(Y)[1]
		else y <- as.vector(Y)
		n <- rep(1, length(Y))
             }
        y<-y/n
        w<-w*n
                                   }
we<-eval(Fam$weight,local=T)
if(!any(is.na(mu))) {
	             devr<- Fam$deviance(mu, y, w)
		    }
			else {
			      devr <- NA
			     }
	if(nrow(X) > p)
		phi <- devr/nrow(X)
	else phi <- NA
	famname <- Fam$family["name"]
	if(is.null(famname))
		famname <- "Gaussian"
		dispersion <- switch(famname,
			Poisson = 1,
			Binomial = 1,
                        Gamma = (sum(((y-mu)/mu)**2))/df.residual,
			phi)
names(dispersion) <- famname
cov<-solve(t(X)%*%W%*%X)
cov<-cov%*%(diag(1,ncol(X))-(t(Cres)%*%
solve(Cres%*%cov%*%t(Cres))%*%Cres%*%cov))
se<-sqrt(diag(cov))*dispersion
coefs<-as.vector(beta)
lp<-as.vector(neta)-offset
fv<-as.vector(mu)
wei<-as.vector(we)
work<-as.vector((y-mu)*Fam$deriv(mu))
this.call <- match.call()
y<-as.vector(y)
null.dev<-fr$null.deviance
dn <- labels(x)
xn <- dn[[2]]
yn <- dn[[1]]
names(coefs) <- xn
names(work) <- yn
names(fv)<-yn
names(lp)<-yn
names(wei)<-yn
names(se)<-xn
dimnames(cov)<-list(xn, xn)
	fit <- list(coefficients = coefs,residuals=work,fitted.values=fv,
                    standard.error=se,cov.cond.unsc=cov,rank = fr$rank,
                    assign = attr(X,"assign"),df.residual = df.residual,
                   weights=wei)
if(length(attributes(w)) | any(w != w[1])) fit$prior.weights<-w
if(fr$rank < p) {
		if(df.residual > 0)
			fit$assign.residual <- fr$assign.residual
		fit$R.assign <- fr$R.assign
		fit$x.assign <- attr(X, "assign")
	}
	if(qr)
		fit$qr <-qr(X) 
		c(fit, list(family = Fam$family, linear.predictors = lp, 
		deviance = devr,null.deviance = null.dev, call =this.call,
                iter = k, y = y,contrasts = attr(X, "contrasts"), 
                dispersion=dispersion))
}
}