library(MASS)	

	
	####################################################.
	####################################################.	
	##FUNCTION INDIVIDUAL RANDOM EFFECTS MODEL (farem)##.
	##FOR COMPUTE THE GRADIENT, HESSIAN AND SIGMA^2   ##
	####################################################.
	####################################################.

	farem<-function(betanew,data) { 

	#We suposse the next structure in the dataset:Id (Identification of individual),. 
	#Group (Group number), Y (Individual outcome: 1 death, 0 alive), O (Population's.
	#observed deaths), n (Population at risk), Z1 (covariate Z1), Z2 (covariate Z2),.
	#...,Zp (covariate Zp).

	Dataprove=data

	#K is the number of groups. We suppose that groups are ordered and they have all.
	#the correlatives numbers. For example:1,2,3 and not 1,3 (There is no number 2).
	#N is the number of observations and p is the number of covariates.
	N<-dim(Dataprove)[1]
	K<-Dataprove[N,2]
	p<-dim(Dataprove)[2]-5
	
	ngr<-matrix(,nrow=1,ncol=K)
	for (i in 1:K) ngr[1,i]<-dim(subset(Dataprove,Dataprove[2]==i))[1]

	gamma0<-betanew[1]
	beta=betanew[-1]

	#Individual outcome.	
	Yki<-matrix(,nrow=N,ncol=1)
	Yki[,1]<-Dataprove[,3]

	#Individual mean.	
	muki<-matrix(,nrow=N,ncol=1)
	muki[,1]<-exp(gamma0+as.vector(as.matrix(Dataprove[,6:(5+p)])%*%beta))

	#Matrix D for the IRM.	
	Dki<-matrix(,nrow=N,ncol=p+1)
	Dki[,1]<-muki[,1]
	for (j in 1:p){
	Dki[,j+1]<-as.numeric(Dataprove[,5+j])*muki[,1]
	}

	#Inverse variance-covariance matrix for the IRM.

		##First, we compute sigma square.		
		muki2<-matrix(,nrow=N,ncol=1)
		muki2[,1]<-muki[,1]^2

		Yaver<-matrix(,nrow=1,ncol=K)
		muk<-matrix(,nrow=1,ncol=K)
		phik<-matrix(,nrow=K,ncol=1)
		sigma2rek<-matrix(,nrow=K,ncol=1)
		
		ini<-1
		end<-ngr[1]
		for (i in 1:K) {
			Yaver[1,i]<-sum(Yki[ini:end])/ngr[i]
			muk[1,i]<-sum(muki[ini:end])/ngr[i]
			phik[i,1]<-sum(muki2[ini:end])/ngr[i]
			sigma2rek[i,]<-max((Yaver[1,i]*(Yaver[1,i]*ngr[i]-2*ngr[i]*muk[1,i]-1)+2*((t(muki[ini:end,1])%*%Yki[ini:end,1])/ngr[i]))/(ngr[i]*(muk[1,i]^2)-phik[i,1])+1,-100)
			ini<-end+1
			end<-ngr[i+1]+end}

		sigma2re<-sum(sigma2rek[1:K])/K

		#We compute the expression for one part (transpose(muk)*Inverse(Deltak)*muk).
		sumk<-matrix(,nrow=1,ncol=K)
		ini<-1
		end<-ngr[1]
		for (i in 1:K) {sumk[1,i]<-sum((muki[ini:end]^2)/(muki[ini:end]*(1-(1+sigma2re)*muki[ini:end])))
				ini<-end+1
				end<-ngr[i+1]+end}
		
	   #Finally, we define the elements of the inverse of V.	
	   Vki<-list(matrix(,nrow=K,ncol=1))

	   for (j in 1:K){
	     Vki[[j]]<-matrix(0,nrow=ngr[1,j],ncol=ngr[1,j])

	     for (i in 1:ngr[1,j]) {			
		yy=1-(1+sigma2re)*muki[i]
                Vki[[j]][i,(i:ngr[1,j])]<-(-sigma2re*(1/yy)*(1/(1-(1+sigma2re)*muki[(i:ngr[1,j])])))*(1/(1+sigma2re*sumk[1,j]))
 	     	Vki[[j]][i,i]<-(1/(muki[i]*yy))-(sigma2re*(1/yy)^2)*(1/(1+sigma2re*sumk[1,j]))
	        }
		Vki[[j]]=Vki[[j]]+t(Vki[[j]])-diag(diag(Vki[[j]]))
            }

	###########################.
	#Gradient, Hessian for IRM#.
	###########################.

	#Vectors of individual responses for each group. 	
	Ykilist<-list(matrix(,nrow=K,ncol=1))
	
	#Vectors of mean for each group.
	mukilist<-list(matrix(,nrow=K,ncol=1))
	
	ini<-1
	end<-ngr[1]
	for (i in 1:K) {Ykilist[i]<-list(matrix(Dataprove[ini:end,3],nrow=ngr[i],ncol=1))
			mukilist[i]<-list(matrix(muki[ini:end],nrow=ngr[i],ncol=1))
		        ini<-end+1
			end<-ngr[i+1]+end}
	
	#Vector diference response and mean.
	Ykiminusmuki<-list(matrix(,nrow=K,ncol=1))
	for (i in 1:K) {Ykiminusmuki[[i]]<-(Ykilist[[i]]-mukilist[[i]])}

	#Matrix Dk for each group.
	Dkilist<-list(matrix(,nrow=K,ncol=1))
        for (i in 1:K) Dkilist[[i]]<-matrix(,nrow=ngr[1,i],ncol=p+1)
	ini<-1
	end<-ngr[1]
 	for (j in 1:K){
		for (n in 1:(p+1)){
		    Dkilist[[j]][,n]<-Dki[ini:end,n]}		      
		    ini<-end+1
		    end<-ngr[j+1]+end}

	#Gradient.  
	ElementKGr<-list(matrix(,nrow=K,ncol=1))
	Grlist<-list(matrix(0,nrow=(p+1),ncol=1))
	for (i in 1:K) {
		ElementKGr[[i]]<-t(Dkilist[[i]])%*%Vki[[i]]%*%Ykiminusmuki[[i]]
		Grlist[[1]]<-Grlist[[1]]+ElementKGr[[i]]}

	Gr<-matrix(,nrow=(p+1),ncol=1)
        for (i in 1:(p+1)) {Gr[i,1]<-Grlist[[1]][i]}
	
	#Hessian.
	ElementKHs<-list(matrix(,nrow=K,ncol=1))
	Hslist<-list(matrix(0,nrow=(p+1),ncol=(p+1)))
	for (i in 1:K) {
		ElementKHs[[i]]<--1*t(Dkilist[[i]])%*%Vki[[i]]%*%Dkilist[[i]]
		Hslist[[1]]<-Hslist[[1]]+ElementKHs[[i]]}

	Hs<-matrix(,nrow=(p+1),ncol=(p+1))
        for (i in 1:(p+1)){
	    for (j in 1:(p+1)) {Hs[i,j]<-Hslist[[1]][i,j]}}

	#We return the Gradient (Gr), Hessian(Hs) and sigma^2(sigma2re).
	

	list(Gradient=Gr,Hessian=Hs,sigma2re=sigma2re)}


	#####################################################.
	#####################################################.	
	##FUNCTION AGGREGATED RANDOM EFFECTS MODEL (fagrem)##.
	##FOR COMPUTE THE GRADIENT, HESSIAN AND SIGMA^2    ##
	#####################################################.
	#####################################################.

	
	fagrem<-function(betanew,data) { 

	#We suposse the next structure in the dataset:Id (Identification of individual),. 
	#Group (Group number), Y (Individual outcome: 1 death, 0 alive), O (Population's.
	#observed deaths), n (Population at risk), Z1 (covariate Z1), Z2 (covariate Z2),.
	#...,Zp (covariate Zp).

	Dataprove=data

	#K is the number of groups. We suppose that groups are ordered and they have all.
	#the correlatives numbers. For example:1,2,3 and not 1,3 (There is no number 2).
	#N is the number of observations and p is the number of covariates.
	N<-dim(Dataprove)[1]
	K<-Dataprove[N,2]
	p<-dim(Dataprove)[2]-5

	ngr<-matrix(,nrow=1,ncol=K)
	for (i in 1:K) ngr[1,i]<-dim(subset(Dataprove,Dataprove[2]==i))[1]

	gamma0<-betanew[1]
	beta<-betanew[-1]

	#Individual mean (muki).	
	muki<-matrix(,nrow=N,ncol=1)
	muki[,1]<-exp(gamma0+as.vector(as.matrix(Dataprove[,6:(5+p)])%*%beta))

	#Individual matrix D.	
	Dki<-matrix(,nrow=N,ncol=p+1)
	Dki[,1]<-muki[,1]
	for (j in 1:p){
	Dki[,j+1]<-as.numeric(Dataprove[,5+j])*muki[,1]
	}

	#Variance for the ARM.
	muki2<-matrix(,nrow=N,ncol=1)
	muki2[,1]<-muki[,1]^2

	#Outcome for the ARM as defined in Sheppard and Prentice (Biometrics,1995).
	Y<-matrix(,nrow=1,ncol=K)

	muk<-matrix(,nrow=1,ncol=K)
	phik<-matrix(,nrow=K,ncol=1)
	#Matrix D for the ARM. 
	Dk<-matrix(,nrow=K,ncol=p+1)
	
	#First, we compute sigma square.
	sigma2amk<-matrix(,nrow=K,ncol=1)
		
		ini<-1
		end<-ngr[1]
		for (i in 1:K) {
			Y[1,i]<-((Dataprove[ini,4])/(Dataprove[ini,5]))
			muk[1,i]<-sum(muki[ini:end])/ngr[i]
			phik[i,1]<-sum(muki2[ini:end])/ngr[i]
			for (j in 1:(p+1)) {Dk[i,j]<-sum(Dki[ini:end,j])/ngr[i]}
			sigma2amk[i,]<-max(((Y[,i]-muk[,i])^2-(muk[,i]-phik[i,])/(Dataprove[ini,5]))/(muk[,i]^2-phik[i,]*(1/(Dataprove[ini,5]))),-100)
			ini<-end+1
			end<-ngr[i+1]+end}

	sigma2am<-sum(sigma2amk[1:K])/K

	#Finally, we define the variance.
	Vk<-matrix(,nrow=1,ncol=K)
	ini<-1
	end<-ngr[1]
	for (i in 1:K) {Vk[1,i]<-sigma2am*((muk[i]^2)-(phik[i,]/(Dataprove[ini,5])))+(muk[i]-phik[i,])*(1/(Dataprove[ini,5]))
			ini<-1
			end<-ngr[1]}
	
	###########################.
	#Gradient, Hessian for ARM#.
	###########################.
		
	#Gradient.  
	Dkt<-t(Dk)
	ElementKaGr<-list(matrix(,nrow=K,ncol=1))
	Gralist<-list(matrix(0,nrow=(p+1),ncol=1))
	for (i in 1:K) {ElementKaGr[[i]]<-Dkt[,i]*(1/Vk[i])*(Y[i]-muk[i])
	Gralist[[1]]<-Gralist[[1]]+ElementKaGr[[i]]}

	Gra<-matrix(,nrow=(p+1),ncol=1)
        for (i in 1:(p+1)) {Gra[i,1]<-Gralist[[1]][i]}

	#Hessian.
	ElementKaHs<-list(matrix(,nrow=K,ncol=1))
	Hsalist<-list(matrix(0,nrow=(p+1),ncol=(p+1)))
	for (i in 1:K) {ElementKaHs[[i]]<--1*matrix(Dkt[,i],nrow=(p+1),ncol=1)%*%(1/Vk[i])%*%Dk[i,]
			Hsalist[[1]]<-Hsalist[[1]]+ElementKaHs[[i]]}

	Hsa<-matrix(,nrow=(p+1),ncol=(p+1))
        for (i in 1:(p+1)){
	    for (j in 1:(p+1)) {Hsa[i,j]<-Hsalist[[1]][i,j]}}

	#We return the Gradient (Gra), Hessian(Hsa) and sigma^2(sigma2am).
	
	list(Gradienta=Gra,Hessiana=Hsa,sigma2am=sigma2am)}


	########################################################.
	########################################################.	
	##FUNCTION POPULATION-BASED ESTIMATING EQUATION (fpbm)##.
	##FOR COMPUTE THE GRADIENT, HESSIAN AND SIGMAS ^ 2    ##
	########################################################.
	########################################################.

	fpbm<-function(betanew,data) { 

	#We suposse the next structure in the dataset:Id (Identification of individual),. 
	#Group (Group number), Y (Individual outcome: 1 death, 0 alive), O (Population's.
	#observed deaths), n (Population at risk), Z1 (covariate Z1), Z2 (covariate Z2),.
	#...,Zp (covariate Zp).

	Dataprove=data

	#K is the number of groups. We suppose that groups are ordered and they have all.
	#the correlatives numbers. For example:1,2,3 and not 1,3 (There is no number 2).
	#N is the number of observations and p is the number of covariates.
	N<-dim(Dataprove)[1]
	K<-Dataprove[N,2]
	p<-dim(Dataprove)[2]-5

	ngr<-matrix(,nrow=1,ncol=K)
	for (i in 1:K) ngr[1,i]<-dim(subset(Dataprove,Dataprove[2]==i))[1]
	
	gamma0<-betanew[1]
	beta<-betanew[-1]

	
	###########################################.
	#Gradient, Hessian for the individual part#.
	###########################################.

	#Individual outcome.	
	Yki<-matrix(,nrow=N,ncol=1)
	Yki[,1]<-Dataprove[,3]

	#Individual mean.	
	muki<-matrix(,nrow=N,ncol=1)
	muki[,1]<-exp(gamma0+as.vector(as.matrix(Dataprove[,6:(5+p)])%*%beta))

	#Individual matrix D.	
	Dki<-matrix(,nrow=N,ncol=p+1)
	Dki[,1]<-muki[,1]
	for (j in 1:p){
	Dki[,j+1]<-as.numeric(Dataprove[,5+j])*muki[,1]
	}

	#Inverse variance-covariance matrix individual part.
		#First, we compute sigma square for the individual part.

		muki2<-matrix(,nrow=N,ncol=1)
		muki2[,1]<-muki[,1]^2

		Yaver<-matrix(,nrow=1,ncol=K)
		muk<-matrix(,nrow=1,ncol=K)
		phik<-matrix(,nrow=K,ncol=1)
		sigma2rek<-matrix(,nrow=K,ncol=1)
		
		ini<-1
		end<-ngr[1]
		for (i in 1:K) {
			Yaver[1,i]<-sum(Yki[ini:end])/ngr[i]
			muk[1,i]<-sum(muki[ini:end])/ngr[i]
			phik[i,1]<-sum(muki2[ini:end])/ngr[i]
			sigma2rek[i,]<-max((Yaver[1,i]*(Yaver[1,i]*ngr[i]-2*ngr[i]*muk[1,i]-1)+2*((t(muki[ini:end,1])%*%Yki[ini:end,1])/ngr[i]))/(ngr[i]*(muk[1,i]^2)-phik[i,1])+1,-100)
			ini<-end+1
			end<-ngr[i+1]+end}

		sigma2re<-sum(sigma2rek[1:K])/K

		#We compute the expression for one part (transpose(muk)*Inverse(Deltak)*muk).
		sumk<-matrix(,nrow=1,ncol=K)
		ini<-1
		end<-ngr[1]
		for (i in 1:K) {sumk[1,i]<-sum((muki[ini:end]^2)/(muki[ini:end]*(1-(1+sigma2re)*muki[ini:end])))
				ini<-end+1
				end<-ngr[i+1]+end}
		
	   #Finally we define the elements of the inverse of V.	
	   Vki<-list(matrix(,nrow=K,ncol=1))

	   for (j in 1:K){
	     Vki[[j]]<-matrix(0,nrow=ngr[1,j],ncol=ngr[1,j])

	     for (i in 1:ngr[1,j]) {			
		yy=1-(1+sigma2re)*muki[i]
	        Vki[[j]][i,(i:ngr[1,j])]<-(-sigma2re*(1/yy)*(1/(1-(1+sigma2re)*muki[(i:ngr[1,j])])))*(1/(1+sigma2re*sumk[1,j]))
                Vki[[j]][i,i]<-(1/(muki[i]*yy))-(sigma2re*(1/yy)^2)*(1/(1+sigma2re*sumk[1,j]))
 	     }
		Vki[[j]]=Vki[[j]]+t(Vki[[j]])-diag(diag(Vki[[j]]))
            }

	#Vectors of individual responses for each group. 
	Ykilist<-list(matrix(,nrow=K,ncol=1))
	
	#Vectors of mean for each group.
	mukilist<-list(matrix(,nrow=K,ncol=1))
	
	ini<-1
	end<-ngr[1]
	for (i in 1:K) {Ykilist[i]<-list(matrix(Dataprove[ini:end,3],nrow=ngr[i],ncol=1))
			mukilist[i]<-list(matrix(muki[ini:end],nrow=ngr[i],ncol=1))
		        ini<-end+1
			end<-ngr[i+1]+end}
	
	#Vector diference response and mean.
	Ykiminusmuki<-list(matrix(,nrow=K,ncol=1))
	for (i in 1:K) {Ykiminusmuki[[i]]<-(Ykilist[[i]]-mukilist[[i]])}

	#Matrix Dk for each group.
	Dkilist<-list(matrix(,nrow=K,ncol=1))
        for (i in 1:K) Dkilist[[i]]<-matrix(,nrow=ngr[1,i],ncol=p+1)
	ini<-1
	end<-ngr[1]
 	for (j in 1:K){
		for (n in 1:(p+1)){
		    Dkilist[[j]][,n]<-Dki[ini:end,n]}		      
		    ini<-end+1
		    end<-ngr[j+1]+end}

	#Gradient Individual part.  
	ElementKGr<-list(matrix(,nrow=K,ncol=1))
	Grlist<-list(matrix(0,nrow=(p+1),ncol=1))
	for (i in 1:K) {
		ElementKGr[[i]]<-t(Dkilist[[i]])%*%Vki[[i]]%*%Ykiminusmuki[[i]]
		Grlist[[1]]<-Grlist[[1]]+ElementKGr[[i]]}

	Gr<-matrix(,nrow=(p+1),ncol=1)
        for (i in 1:(p+1)) {Gr[i,1]<-Grlist[[1]][i]}
	
	##Hessian individual part.
	ElementKHs<-list(matrix(,nrow=K,ncol=1))
	Hslist<-list(matrix(0,nrow=(p+1),ncol=(p+1)))
	for (i in 1:K) {
		ElementKHs[[i]]<--1*t(Dkilist[[i]])%*%Vki[[i]]%*%Dkilist[[i]]
		Hslist[[1]]<-Hslist[[1]]+ElementKHs[[i]]}

	Hs<-matrix(,nrow=(p+1),ncol=(p+1))
        for (i in 1:(p+1)){
	    for (j in 1:(p+1)) {Hs[i,j]<-Hslist[[1]][i,j]}}
	


	###########################################.
	#Gradient, Hessian for the aggregated part#.
	###########################################.


	#Outcome aggregated part.
	Ybar<-matrix(,nrow=1,ncol=K)

	#Matrix D for the aggregated part. 
	Dk<-matrix(,nrow=K,ncol=p+1)

        ini<-1
	end<-ngr[1]
	for (i in 1:K) {Ybar[1,i]<-((Dataprove[ini,4]-sum(Yki[ini:end]))/(Dataprove[ini,5]-ngr[i]))
			for (j in 1:(p+1)) {Dk[i,j]<-sum(Dki[ini:end,j])/ngr[i]}
			ini<-end+1
			end<-ngr[i+1]+end}

	#Sigma square aggregated part.
	sigma2pbk<-matrix(,nrow=K,ncol=1)
	ini<-1
	end<-ngr[1]
	for  (i in 1:K) {sigma2pbk[i,]<-max(((Ybar[,i]-muk[,i])^2-(muk[,i]-phik[i,])/(Dataprove[ini,5]-ngr[i]))/(muk[,i]^2-phik[i,]*(1/(Dataprove[ini,5]-ngr[i]))),0)
			ini<-end+1
			end<-ngr[i+1]+end}

	sigma2pb<-sum(sigma2pbk[1:K])/K

	#Variance aggregated part.
	Dkt<-t(Dk)	
	Vkbar<-matrix(,nrow=1,ncol=K)
	ElementKarGr<-list(matrix(,nrow=K,ncol=1))
	Grarlist<-list(matrix(0,nrow=(p+1),ncol=1))
	ini<-1
	end<-ngr[1]
	for (i in 1:K) {Vkbar[1,i]<-sigma2pb*((muk[i]^2)-(phik[i,]/(Dataprove[ini,5]-ngr[i])))+(muk[i]-phik[i,])*(1/(Dataprove[ini,5]-ngr[i]))
			ElementKarGr[[i]]<-Dkt[,i]*(1/Vkbar[i])*(Ybar[i]-muk[i])
			Grarlist[[1]]<-Grarlist[[1]]+ElementKarGr[[i]]
			ini<-end+1
			end<-ngr[i+1]+end}

	#Gradient aggregated part.
	Grar<-matrix(,nrow=(p+1),ncol=1)
	for (i in 1:(p+1)) {Grar[i,1]<-Grarlist[[1]][i]}

	#Hessian aggregated part.

	ElementKarHs<-list(matrix(,nrow=K,ncol=1))
	Hsarlist<-list(matrix(0,nrow=(p+1),ncol=(p+1)))
	for (i in 1:K) {ElementKarHs[[i]]<--1*matrix(Dkt[,i],nrow=(p+1),ncol=1)%*%(1/Vkbar[i])%*%Dk[i,]
			Hsarlist[[1]]<-Hsarlist[[1]]+ElementKarHs[[i]]}

	Hsar<-matrix(,nrow=(p+1),ncol=(p+1))
        for (i in 1:(p+1)){
	    for (j in 1:(p+1)) {Hsar[i,j]<-Hsarlist[[1]][i,j]}}

	######################################################################.
	#Gradient, Hessian for the individual and aggregated part combination#.
	######################################################################.

	#Gradient PBEE.
	
	Grpb<-(Gr+Grar)

	#Hessian PBEE.

	Hspb<-(Hs+Hsar)	

	#We return the Gradient (Grpb), Hessian(Hspb), sigma^2 individual. 
	#part (sigma2re) and sigma^2 aggregated part (sigma2pb).

	list(Gradientpb=Grpb,Hessianpb=Hspb,sigma2re=sigma2re,sigma2pb=sigma2pb)}


	
	#######################################################################.
	#######################################################################.
	##PARAMETER ESTIMATES USING THE NEWTON-RAPHSON ALGORITHM FOR THE IRM ##.
	##(function farem2), ARM (function fagrem2) AND PBEE (function fpbm2)##.
	#######################################################################.
	#######################################################################.

	#We suposse the next structure in the dataset:Id (Identification of individual),. 
	#Group (Group number), Y (Individual outcome: 1 death, 0 alive), O (Population's.
	#observed deaths), n (Population at risk), Z1 (covariate Z1), Z2 (covariate Z2),.
	#...,Zp (covariate Zp).

	farem2<-function(data,tol,maxiter=100,betainitial){

	betanew<-betainitial
	betaold<-betanew+1
	itercount=0
	Hdim=length(betanew)^2

	#Solutions for the IRM.

	while(max(abs((betaold-betanew)/betaold))>tol){
	     q<-farem(betanew,data)	     
	     betaold<-betanew
	     if(sum(is.finite(q$Hessian))<Hdim) itercount=999
	     if(itercount>maxiter) break
	     betanew<-betanew-ginv(q$Hessian)%*%q$Gradient	     
	     itercount=itercount+1
	     if(is.na(q$sigma2re)) itercount=100
	     if(itercount>maxiter) break
	     if(q$sigma2re>50) itercount=100
	     if(itercount>maxiter) break
	     }
	list(betanew=betanew,sigma2re=q$sigma2re, iterN=itercount)
        }

	fagrem2<-function(data,tol,maxiter=100,betainitial){

	betanewa<-betainitial
	betaold<-betanewa+1
	itercount=0
	Hdim=length(betanewa)^2

	#Solutions for the ARM.

	while(max(abs((betaold-betanewa)/betaold))>tol){
	     d<-fagrem(betanewa,data)	     
	     betaold<-betanewa
	     if(sum(is.finite(d$Hessiana))<Hdim) itercount=999
	     if(itercount>maxiter) break
	     betanewa<-betanewa-ginv(d$Hessiana)%*%d$Gradienta	     
	     itercount=itercount+1
	     if(is.na(d$sigma2am)) itercount=100
	     if(itercount>maxiter) break
	     if(d$sigma2am>50) itercount=100
	     if(itercount>maxiter) break
	     }
	list(betanewa=betanewa,sigma2am=d$sigma2am, iterN=itercount)
        }

	fpbm2<-function(data,tol,maxiter=100,betainitial){
	
	betanewpb<-betainitial
	betaold<-betanewpb+1
	itercount=0
	Hdim=length(betanewpb)^2

	#Solutions for the PBEE.

	while(max(abs((betaold-betanewpb)/betaold))>tol){
	     e<-fpbm(betanewpb,data)	     
	     betaold<-betanewpb
	     if(sum(is.finite(e$Hessianpb))<Hdim) itercount=999
	     if(itercount>maxiter) break
	     betanewpb<-betanewpb-ginv(e$Hessianpb)%*%e$Gradientpb	     
	     itercount=itercount+1
	     if(is.na(max(e$sigma2re,e$sigma2pb))) itercount=100
	     if(itercount>maxiter) break
	     if(max(e$sigma2re,e$sigma2pb)>50) itercount=100
	     if(itercount>maxiter) break
	     }
	list(betanewpb=betanewpb,sigma2re=e$sigma2re,sigma2pb=e$sigma2pb, iterN=itercount)
        }
	
	###############################################.
	###############################################.	
	##FUNCTION FOR SIMULATE THE DATA (fgenerate2)##.
	##IN THE SIMPLE CONFOUNDING CASE (SCC)       ##
	###############################################.
	###############################################.
	
	fgenerate2<-function(group,populationsize,samplesize,variance){

	#group=number of groups, populationsize= population size.
	#samplesize= sample size, variance=within group variance.
		
	K<-group
	nk<-populationsize
	mk<-samplesize	
	varwithin<-variance	

	#Covariate X1ki & X2ki.  They are correlated 0.3 at the community.
	#and individual levels (see Prentice & Sheppard).

	zK=mvrnorm(K,c(0,0),matrix(c(1,.3,.3,1),2,2))

	cov=0.3*sqrt(varwithin)*1
	covm=matrix(c(varwithin,cov,cov,1),2,2)

	X1ki=matrix(0,nk,K)
	X2ki=matrix(0,nk,K)

	for(i in 1:K){
	 znk=mvrnorm(nk,zK[i,],covm)
	 X1ki[,i]=znk[,1]
	 X2ki[,i]=znk[,2]
	}
	
	#Country specific frailties were generated as independent.
	#realized  values from a gamma  distribution with  mean 1.
	#and variance sigma^2. The mean of a gamma is shape*scale.
	#and the variance is shape*(scale^2).

	meanhk<-1
	varhk<-0.05
	shape<-(meanhk^2)/(varhk)
	scale<-(varhk)/(meanhk)
	hk<-rgamma(K,shape=shape,scale=scale)[]
	hk=t(matrix(rep(hk,nk),nrow=K,ncol=nk))

	#The disease events, yki, were generated by determining.
	#wether  a  uniform  random  variable  wass  less  than.
	#hk*exp(gamma0+beta1*X1ki+beta2*X2ki).
	
	gamma0<--3
	beta1<-0.2
	beta2<-0.2

	yki<-matrix(,nrow=nk,ncol=K)
	unif<-matrix(runif(nk*K,0,1),nrow=nk,ncol=K)
	
        yki=ifelse(unif<hk*exp(gamma0+beta1*X1ki+beta2*X2ki),1,0)

	#######################################.
	#selection of random sample of size mk#.
	#and organize data to apply functions #.
	#farem, fagrem and fpbm               #
	#######################################.

	datalist<-list(matrix(,nrow=K,ncol=1))
	sampledatalist<-list(matrix(,nrow=K,ncol=1))
	ini<-1
	end<-mk
	data<-matrix(,nrow=mk*K,ncol=5)
	for (i in 1:K){

	datalist[[i]]<-cbind(matrix(c(1:nk),nrow=nk,ncol=1),matrix(yki[,i],nrow=nk,ncol=1),matrix(X1ki[,i],nrow=nk,ncol=1),matrix(X2ki[,i],nrow=nk,ncol=1),matrix(c(i),nrow=nk,ncol=1))
	sampledatalist[[i]]<-datalist[[i]][as.matrix(sample(datalist[[i]][,1],mk)),]
	data[ini:end,]<-sampledatalist[[i]][]
			ini<-end+1
			end<-mk*(i+1)
	}

	O<-matrix(apply(yki,2,sum),nrow=K,ncol=1)

	ini<-1
	end<-mk
	datapop<-matrix(,nrow=mk*K,ncol=1)
	for (i in 1:K) {datapop[ini:end,1]<-O[i,]
		        ini<-end+1
		        end<-mk*(i+1)}

	datadatapop<-cbind(data,datapop,c(nk))
	datafin<-data.frame(id=matrix(datadatapop[,1]),group=matrix(datadatapop[,5]),YIND=matrix(datadatapop[,2]),O=matrix(datadatapop[,6]),n=matrix(datadatapop[,7]),X1ki=matrix(datadatapop[,3]),X2ki=matrix(datadatapop[,4]))}


	######################.
	##Simulation results##.
	######################.


	fsimulationB<-function(seed,Niter,sigma2,K,N){
	set.seed(seed)
	count1=0
	result1=matrix(0,nrow=Niter,ncol=16)

	while(count1<Niter){
	tempdata1=fgenerate2(K,2000,N,sigma2)  
	tol<-0.001                               
	gg<-glm(YIND~X1ki+X2ki,data=tempdata1,family=binomial)
	betaini<-as.vector(c(matrix(gg$coefficients[1]),matrix(gg$coefficients[2]),matrix(gg$coefficients[3])))	
	tempdata1a1<-farem2(tempdata1,tol,maxiter=50,betaini)
	tempdata1a2<-fagrem2(tempdata1,tol,maxiter=50,betaini)
	tempdata1a3<-fpbm2(tempdata1,tol,maxiter=50,betaini)
	count1=count1+1
	result1[count1,1:5]=matrix(c(t(tempdata1a1$betanew),tempdata1a1$sigma2re,tempdata1a1$iterN),nrow=1,ncol=5)
	result1[count1,6:10]=matrix(c(t(tempdata1a2$betanewa),tempdata1a2$sigma2am,tempdata1a2$iterN),nrow=1,ncol=5)
	result1[count1,11:16]=matrix(c(t(tempdata1a3$betanewpb),tempdata1a3$sigma2re,tempdata1a3$sigma2pb,tempdata1a3$iterN),nrow=1,ncol=6)
	print(count1)	
	}

	list(result=result1,sigma2=sigma2,K=K,N=N)
	}


	#100 groups-100 sample size in each group.

	finalresultB100100B.25=fsimulationB(123,1000,.25,100,100)
	save(list=c("finalresultB100100B.25",".Random.seed"),file="100100B025.RData")
	savedseed=.Random.seed

	finalresultB100100B.5=fsimulationB(savedseed,1000,.5,100,100)
	save(list=c("finalresultB100100B.5",".Random.seed"),file="100100B05.RData")
	savedseed=.Random.seed

	finalresultB100100B1=fsimulationB(savedseed,1000,1,100,100)
	save(list=c("finalresultB100100B1",".Random.seed"),file="100100B1.RData")
	savedseed=.Random.seed

	finalresultB100100B2=fsimulationB(savedseed,1000,2,100,100)
	save(list=c("finalresultB100100B2",".Random.seed"),file="100100B2.RData")
	savedseed=.Random.seed

	finalresultB100100B4=fsimulationB(savedseed,1000,4,100,100)
	save(list=c("finalresultB100100B4",".Random.seed"),file="100100B4.RData")
	savedseed=.Random.seed

	finalresultB100100B8=fsimulationB(savedseed,1000,8,100,100)
	save(list=c("finalresultB100100B8",".Random.seed"),file="100100B8.RData")
	savedseed=.Random.seed

	finalresultB100100B16=fsimulationB(savedseed,1000,16,100,100)
	save(list=c("finalresultB100100B16",".Random.seed"),file="100100B16.RData")


	#100 groups-50 sample size in each group.

	finalresultB10050B.25=fsimulationB(123,1000,.25,100,50)
	save(list=c("finalresultB10050B.25",".Random.seed"),file="10050B025.RData")
	savedseed=.Random.seed

	finalresultB10050B.5=fsimulationB(savedseed,1000,.5,100,50)
	save(list=c("finalresultB10050B.5",".Random.seed"),file="10050B05.RData")
	savedseed=.Random.seed

	finalresultB10050B1=fsimulationB(savedseed,1000,1,100,50)
	save(list=c("finalresultB10050B1",".Random.seed"),file="10050B1.RData")
	savedseed=.Random.seed

	finalresultB10050B2=fsimulationB(savedseed,1000,2,100,50)
	save(list=c("finalresultB10050B2",".Random.seed"),file="10050B2.RData")
	savedseed=.Random.seed

	finalresultB10050B4=fsimulationB(savedseed,1000,4,100,50)
	save(list=c("finalresultB10050B4",".Random.seed"),file="10050B4.RData")
	savedseed=.Random.seed

	finalresultB10050B8=fsimulationB(savedseed,1000,8,100,50)
	save(list=c("finalresultB10050B8",".Random.seed"),file="10050B8.RData")
	savedseed=.Random.seed

	finalresultB10050B16=fsimulationB(savedseed,1000,16,100,50)
	save(list=c("finalresultB10050B16",".Random.seed"),file="10050B16.RData")


	#50 groups-100 sample size in each group.

	finalresultB50100B.25=fsimulationB(123,1000,.25,50,100)
	save(list=c("finalresultB50100B.25",".Random.seed"),file="50100B025.RData")
	savedseed=.Random.seed

	finalresultB50100B.5=fsimulationB(savedseed,1000,.5,50,100)
	save(list=c("finalresultB50100B.5",".Random.seed"),file="50100B05.RData")
	savedseed=.Random.seed

	finalresultB50100B1=fsimulationB(savedseed,1000,1,50,100)
	save(list=c("finalresultB50100B1",".Random.seed"),file="50100B1.RData")
	savedseed=.Random.seed

	finalresultB50100B2=fsimulationB(savedseed,1000,2,50,100)
	save(list=c("finalresultB50100B2",".Random.seed"),file="50100B2.RData")
	savedseed=.Random.seed

	finalresultB50100B4=fsimulationB(savedseed,1000,4,50,100)
	save(list=c("finalresultB50100B4",".Random.seed"),file="50100B4.RData")
	savedseed=.Random.seed

	finalresultB50100B8=fsimulationB(savedseed,1000,8,50,100)
	save(list=c("finalresultB50100B8",".Random.seed"),file="50100B8.RData")
	savedseed=.Random.seed

	finalresultB50100B16=fsimulationB(savedseed,1000,16,50,100)
	save(list=c("finalresultB50100B16",".Random.seed"),file="50100B16.RData")

	
	#50 groups-50 sample size in each group.

	finalresultB5050B.25=fsimulationB(123,1000,.25,50,50)
	save(list=c("finalresultB5050B.25",".Random.seed"),file="5050B025.RData")
	savedseed=.Random.seed

	finalresultB5050B.5=fsimulationB(savedseed,1000,.5,50,50)
	save(list=c("finalresultB5050B.5",".Random.seed"),file="5050B05.RData")
	savedseed=.Random.seed

	finalresultB5050B1=fsimulationB(savedseed,1000,1,50,50)
	save(list=c("finalresultB5050B1",".Random.seed"),file="5050B1.RData")
	savedseed=.Random.seed

	finalresultB5050B2=fsimulationB(savedseed,1000,2,50,50)
	save(list=c("finalresultB5050B2",".Random.seed"),file="5050B2.RData")
	savedseed=.Random.seed

	finalresultB5050B4=fsimulationB(savedseed,1000,4,50,50)
	save(list=c("finalresultB5050B4",".Random.seed"),file="5050B4.RData")
	savedseed=.Random.seed

	finalresultB5050B8=fsimulationB(savedseed,1000,8,50,50)
	save(list=c("finalresultB5050B8",".Random.seed"),file="5050B8.RData")
	savedseed=.Random.seed

	finalresultB5050B16=fsimulationB(savedseed,1000,16,50,50)
	save(list=c("finalresultB5050B16",".Random.seed"),file="5050B16.RData")