"myags.pmat.fit" <- function(X, Y, ID, weights, cor.met, family,
	alpfun, scalefun, wcorigen, tol, 
	maxiter, prog.error)
{
	time1<-proc.time()
	MX <- pmat2mat(X)
	xn <- dimnames(X[[1]])[[2]]
	dimnames(MX) <- list(NULL, xn)
	initrun <- glm(unlist(Y) ~ MX - 1, 
		family = family)
	b0 <- initrun$coef
	cur.rank <- initrun$rank
	del <- 1
	assign("glm.links", yags.links, f
		 = 1)
	iter <- 0
	while(del > tol & iter < maxiter) 
		{
		iter <- iter + 1
		print(paste("Starting iteration",
			iter, "with coefs")
			)
		print(as.numeric(b0))
		eta <- MX %*% b0	# 
		mu <- family()$inverse(eta)
		# 
		d.mu <- family()$deriv(mu)	
	# 
# 1.5.1.1 E use weights
		A <- as.double(family()$
			variance(mu)) * 
			weights	# 
		if(length(A) == 1) A <- rep(
				A, length(
				unlist(Y)))
		# 
		D.i <- load.clustered.design(
			(1/sqrt(A)) * 
			as.double(family()$
			dinv.deta(eta)) * 
			MX, ID)	
	# D scaled by V(mu)
		PR <- (1/sqrt(A)) * (unlist(
			Y) - mu)	
	# pearson resid
		if(!is.null(cor.met) & !
			is.matrix(cor.met)
			)
			TT <- split.preserveord(
				cor.met, ID
				)
		else if(!is.null(cor.met) & 
			is.matrix(cor.met)
			)
			TT <- msplit(
				cor.met, ID
				)
		else TT <- NULL
		alpha <- alpfun(PR, ID, 
			cor.met, scalefun, 
			p = ncol(MX))
		PR <- load.clustered.outcome(
			PR, ID)
		cor.els <- list()
		NC <- length(PR)
		for(i in 1:NC)
			cor.els[[i]] <- 
				list(r = PR[[
				i]], 
				cor.met = 
				TT[[i]], 
				alpha = 
				alpha)
		Rinv.i <- fill(wcorigen, 
			cor.els, c("pmat", 
			"block", "diag"))
		S2 <- solve(sum(t(D.i) %*% 
			Rinv.i %*% D.i))
		delb <- S2 %*% sum(t(D.i) %*% Rinv.i %*% PR)
		del <- max(abs(delb))
		b1 <- b0 + delb
		b0 <- b1
	}
# wrap up using final beta values
	if(iter == maxiter) {
		warning("maximum number of iterations consumed"
			)
		prog.error <- c(prog.error, 
			"niter==maxiter")
	}
	
	eta <- MX %*% b0
	mu <- family()$inverse(eta)
	Rvec <- unlist(Y) - mu	
	# raw residual
	PR <- (1/sqrt(A)) * Rvec
	alpha <- alpfun(PR, ID, cor.met, 
		scalefun, p = ncol(MX))
	phi <- scalefun(PR, ID, p = ncol(MX
		))
	PR <- load.clustered.outcome(PR, ID
		)
	for(i in 1:NC)
		cor.els[[i]] <- list(r = PR[[
			i]], cor.met = TT[[
			i]], alpha = alpha)
	vvt <- function(x)
	matrix(x$r %*% t(x$r), nr = length(
		x$r))
	eet <- fill(vvt, cor.els, c("pmat", 
		"block", "diag"))
	Rinv.i <- fill(wcorigen, cor.els, c(
		"pmat", "block", "diag"))
	S2 <- solve(sum(t(D.i) %*% Rinv.i %*% 
		D.i))
	S5 <- sum(t(D.i) %*% Rinv.i %*% eet %*%
		Rinv.i %*% D.i)
	b1 <- as.numeric(b1)
	names(b1) <- xn
	time2<-proc.time()
##############################################
### modifications by m.p. fay:
### add code to produce modified estimates
    NP<- dim(S5)[1]	
	omega<-array(0,c(NC,NP,NP))
	scoremat<-matrix(0,NC,NP)
	for (i in 1:NC){
		omega[i,,]<- (1/phi)* t(D.i)[[i]] %*% Rinv.i[[i]] %*% D.i[[i]]
		scoremat[i,]<- (1/phi)* t(D.i)[[i]] %*% Rinv.i[[i]] %*% PR[[i]]
	}

### add u and omega onto output list
	list(coefficients = b1, 
		naive.parmvar = phi * S2, 
		robust.parmvar = S2 %*% S5 %*% 
		S2, 		alpha = alpha, phi = 
		phi, linear.predictors = 
		eta, fitted.values = mu, 
		residuals = Rvec, iter = 
		iter, family = family, rank
		 = cur.rank, errorcode = 
		prog.error,u=scoremat,omega=omega)		
#### end modifications by m.p.fay
##############################################
}