### This program shows a data set that does not appear to have a solution to the generalized estimating equations
### The GEE are Poisson with exchangeable correlation structure
###   (The data are simulated data, simulated according to 
###   y_{i1},...y_{i n_i} ~ Poisson(10), i=1,..,20
###   and the number in each cluster n_i ~ ceiling(X_i), where X_i ~ Gamma(mean=10,var=20) 
###   and treat=-1 if i<=10 and treat=1 if  i>10)


### Uses GEE function from YAGS (by VJ CAREY) but modified to input a parameter estimate and not do any iterations
### see myags.onestep.fit and myags.onestep below
library(yags)

"dat106"<-
structure(.Data = list(Y = c(14, 5, 11, 17, 8, 9, 4, 12, 8, 7, 7, 9, 9, 8, 5, 
	17, 10, 8, 8, 14, 8, 12, 9, 11, 5, 6, 5, 6, 10, 11, 7, 10, 10, 6, 11, 
	11, 14, 7, 9, 12, 8, 9, 11, 9, 8, 10, 5, 11, 10, 10, 8, 7, 11, 9, 6, 13,
	7, 11, 3, 11, 22, 9, 6, 11, 10, 5, 9, 11, 9, 13, 14, 7, 13, 11, 10, 12, 
	10, 5, 9, 14, 10, 4, 12, 8, 15, 6, 5, 11, 4, 14, 17, 11, 12, 6, 14, 14, 
	9, 8, 13, 7, 9, 8, 4, 9, 10, 8, 10, 6, 8, 10, 11, 10, 9, 9, 14, 12, 9, 
	6, 13, 7, 8, 12, 14, 11, 15, 5, 15, 7, 7, 8, 12, 10, 8, 11, 5, 13, 7, 
	11, 16, 7, 10, 14, 8, 9, 12, 12, 5, 8, 11, 11, 17, 14, 10, 9, 9, 9, 3, 
	10, 11, 11, 14, 13, 9, 13, 11, 8, 11, 9, 9, 7, 13, 17, 6, 14, 6, 8, 10, 
	10, 15, 8, 11, 13, 12, 10, 11, 9, 7, 10, 8, 18, 3, 10, 12, 11, 11, 11, 
	4, 11, 10, 5, 5, 13, 11, 10, 7, 5, 8, 11, 10, 15, 10, 12, 9, 9, 14, 10, 
	11, 11, 8, 12, 4, 6, 11, 12, 7), strata = c(1, 1, 1, 1, 1, 1, 1, 2, 2, 
	2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 
	4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 
	7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 
	8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 
	10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 
	11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 
	13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 
	14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 
	16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 
	18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 
	20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20), treat = c(-1, -1, 
	-1
	, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1
	, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1
	, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
	-1
	, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
	1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
	1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
	1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), 
	row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", 
	"12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", 
	"24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", 
	"36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", 
	"48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", 
	"60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", 
	"72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", 
	"84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", 
	"96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106",
	"107", "108", "109", "110", "111", "112", "113", "114", "115", "116", 
	"117", "118", "119", "120", "121", "122", "123", "124", "125", "126", 
	"127", "128", "129", "130", "131", "132", "133", "134", "135", "136", 
	"137", "138", "139", "140", "141", "142", "143", "144", "145", "146", 
	"147", "148", "149", "150", "151", "152", "153", "154", "155", "156", 
	"157", "158", "159", "160", "161", "162", "163", "164", "165", "166", 
	"167", "168", "169", "170", "171", "172", "173", "174", "175", "176", 
	"177", "178", "179", "180", "181", "182", "183", "184", "185", "186", 
	"187", "188", "189", "190", "191", "192", "193", "194", "195", "196", 
	"197", "198", "199", "200", "201", "202", "203", "204", "205", "206", 
	"207", "208", "209", "210", "211", "212", "213", "214", "215", "216", 
	"217", "218", "219", "220", "221", "222", "223", "224", "225"), class
	 = "data.frame")



"myags.onestep.fit" <- function(X, Y, ID, weights, cor.met, family,
	alpfun, scalefun, wcorigen, tol, prog.error,b0)
{
## modified from yags.pmat.fit	
	time1<-proc.time()
	MX <- pmat2mat(X)
	xn <- dimnames(X[[1]])[[2]]
	dimnames(MX) <- list(NULL, xn)
	initrun <- glm(unlist(Y) ~ MX - 1, 
		family = family)
	cur.rank <- initrun$rank
	del <- 1
	assign("glm.links", yags.links, f
		 = 1)
		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))
		print(paste("alpha=",alpha))	
		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))
		#	print(paste("Info=",S2))
		U.old<- sum(t(D.i) %*% 	Rinv.i %*% PR)
		print(paste("Score=",U.old))
		print(paste("sumscore=",sum(abs(U.old))))
	#list(coefficients = b0, score=U.old,sumscore=sum(abs(U.old)))
	sum(abs(U.old))
}



"myags.onestep"<-
function(formula, id, weights = NULL, cor.met = NULL, family, alpfun = NULL, 
	scalefun = LZ.scalefun, wcorigen = identni, tol = 0.0001, contrasts = 
	NULL, corstr = c("independence", "exchangeable", "ar1", "unstructured"),startbeta=0)
{
##	modified from yags
	prog.error <- NULL
	call <- match.call()
	m <- match.call(expand = F)
	
	m$contrasts <- m$tol <- m$wcorigen <- m$scalefun <- m$alpfun <- m$
		family <- m$cor.met <- m$id <- m$corstr <- m$maxiter <- m$startbeta<-NULL
	m[[1]] <- as.name("model.frame")
	m <- eval(m, sys.parent())
	Terms <- attr(m, "terms")
	y <- as.matrix(model.extract(m, response))	
	# 1.5.1.1 B : add weight-processing specifically for binomial matrices
# as in glm()
	if(is.null(weights))
		weights <- rep(1, nrow(y))
	else if(length(weights) != nrow(y))
		stop("lengths of y and weights incompatible")
	if(ncol(y) == 2) {
		n <- apply(y, 1, sum)
		y <- y[, 1, drop = F]/n
		weights <- n/weights
	}
	else weights <- 1/weights
	x <- model.matrix(Terms, m, contrasts)
	Y <- load.clustered.outcome(y, id)
	X <- load.clustered.design(x, id)
	if(is.null(alpfun) & length(corstr) > 1) {
		warning("no alpfun provided and corstr has several elements; corstr[1] used"
			)
		corstr <- corstr[1]
	}
	if(!is.null(alpfun) & length(corstr) == 1)
		warning("both alpfun and corstr provided; corstr will be ignored"
			)
	else if(corstr[1] == "exchangeable" & is.null(alpfun)) {
		alpfun <- LZ.exchalp
		wcorigen <- excoriput
	}
	else if(corstr[1] == "ar1" & is.null(alpfun)) {
		if(is.null(cor.met))
			stop("must supply times (cor.met) with ar1")
		alpfun <- prop.ar1alp
		wcorigen <- ar1.coriput
	}
	else if(corstr[1] == "unstructured" & is.null(alpfun)) {
		if(is.null(cor.met))
			stop("must supply times (cor.met) with unstructured")
		alpfun <- prop.unstruc.alp
		wcorigen <- unstr.coriput
	}
	else if(corstr[1] == "independence" & is.null(alpfun)) {
		alpfun <- function(a, b, c, d, p)
		NULL
		wcorigen <- identni
	}
	ans <- myags.onestep.fit(X, Y, id, weights, cor.met, family, alpfun, 
		scalefun, wcorigen, tol, prog.error,startbeta)
#	ans$call <- call
#	class(ans) <- c("myags", "glm", "lm")
	ans
}




CORR<-"exchangeable"
Y<-dat[,"Y"]
strata<-dat[,"strata"]
treat<-dat[,"treat"]		

beta0<-c(10:35)/10
beta1<-(-10:10)/5
out<-matrix(NA,length(beta0),length(beta1))
icnt<-0
for (i in beta0){
	icnt<-icnt+1
	jcnt<-0
	for (j in beta1){
		jcnt<-jcnt+1
	out[icnt,jcnt]<-myags.onestep(Y~treat,id=strata,family=poisson,
scalefun=SCALEFUN,corstr=CORR,startbeta=c(i,j))	

	}
}

dimnames(out)<-list(as.character(beta0),as.character(beta1))
out
out1<-out
out1[out>300]<-300
persp(beta0,beta1,-out1,zlab="-sum | U_i |",xlab="Intercept",ylab="Treat")
title("Minus the Sum of Abs. Value of terms of Estimating Equation\n(bounded below by -300)")