mxxfpkg:

msurv<-function(study, time, n.risk, surv.rate, confidence)
{

.data <- data.frame(study, time, n.risk, surv.rate)
.data <- .data[order(study, time),]
study <- .data$study
time <- .data$time
n.risk <- .data$n.risk
surv.rate <- .data$surv.rate

# Number of studies
IndiceStudies<-sort(unique(study))
NbStudies<-length(IndiceStudies)

# Number of points in time (number of pooled conditional survival to be estimated)
IndiceTimes<-sort(unique(time))
NbTimes<-length(IndiceTimes)
TimeMax<-max(IndiceTimes)

verif <- function(x)
{
ref.study <- sort(time[study==x])
ref.data  <- IndiceTimes[IndiceTimes <= max(ref.study)]
return(1*(sum(ref.study!=ref.data)==0))
}

verif.data <- data.frame(Sstudy=IndiceStudies, check=sapply(IndiceStudies, FUN="verif"))

if (sum(verif.data$check==0)>0) { 

return(list(
verif.data=verif.data,
summary.fixed=NA,
median.fixed=NA,
mean.fixed=NA,
heterogeneity=NA,
summary.random=NA,
median.random=NA,
mean.random=NA))

}

else {

time.init <- min(time[time>0])
CondSurv <- c(-99, surv.rate[2:length(surv.rate)]/surv.rate[1:(length(surv.rate)-1)])
CondSurv[time==time.init]<-surv.rate[time==time.init] # question : cette fonction ne fonctionne que si c'est par pas de 1 an ?

#### Mise en forme des data

# Arc sine tranformation of the conditional survival (Equations page 7)
ArcSineCondSurv<-asin(sqrt((n.risk*CondSurv+0.25)/(n.risk+0.5)))
VarArcSineCondSurv<-1/(4*(n.risk+0.5))

# ys est la matrice des proba conditionnelles arcsine transformées (1 ligne par étude, 1 colonne par temps, valeurs manquantes remplacées par la moyenne des autres)
ys<-matrix(nrow=NbStudies,ncol=NbTimes)
for (i in 1:NbStudies) ys[i,1:sum(study==IndiceStudies[i])]<-ArcSineCondSurv[study==IndiceStudies[i]]
for (i in 1:NbTimes)  ys[is.na(ys[,i]),i]<-mean(ys[!is.na(ys[,i]),i])

# vars est la matrice des variances des proba conditionnelles arcsine transformées (1 ligne par étude, 1 colonne par temps, valeurs manquantes remplacées par 10^12)
vars<-matrix(nrow=NbStudies,ncol=NbTimes)
for (i in 1:NbStudies) vars[i,1:sum(study==IndiceStudies[i])]<-VarArcSineCondSurv[study==IndiceStudies[i]]
for (i in 1:NbTimes)  vars[is.na(vars[,i]),i]<-10^12



#### Application of the multivariate DerSimonian and Laird methodology (R code from Dan Jackson)
#### Equations page 7
#### Estimation of the pooled arc sine transformed conditional survival

	mat<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	maty<-matrix(0, nrow=NbTimes, ncol=1)
	matfixedeffect<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	matyfixedeffect<-matrix(0, nrow=NbTimes, ncol=1)
	esttrun<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	estmat<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	count<-1
	
	for(i in 1:(NbTimes-1))
	{
		for(j in (i+1):NbTimes) 
		{
			crossses<-(vars[,i]*vars[,j])^0.5
			mean1a<-sum(ys[,i]/crossses)/sum(1/crossses)
			mean2a<-sum(ys[,j]/crossses)/sum(1/crossses)
			Q<-sum((ys[,i]-mean1a)*(ys[,j]-mean2a)/crossses)
			bb<-sum(1/crossses)-sum(1/(crossses^2))/sum(1/crossses)
			est<-Q/bb
			estmat[i,j]<-est
		}
	}
	estmat<-estmat+t(estmat)
	for(i in 1:NbTimes)
	{
		cross_ses<-vars[,i]
		mean1a<-sum(ys[,i]/vars[,i])/sum(1/vars[,i])
		Q<-sum((ys[,i]-mean1a)*(ys[,i]-mean1a)/vars[,i])
		rhos<-rep(1,NbStudies)
		check<-(1-1*(vars[,i]==10^12))
		rhos<-rhos*check
		aa<-sum(rhos)-sum(rhos/vars[,i])/sum(1/vars[,i])
		bb<-sum(1/vars[,i])-sum(1/(vars[,i]^2))/sum(1/vars[,i])
		est<-(Q-aa)/bb
		estmat[i,i]<-est
	}
	eig<-eigen(estmat)
	for(i in 1:NbTimes) esttrun<-esttrun+max(0, eig$values[i])*eig$vectors[,i]%*%t(eig$vectors[,i])
	for(k in 1:NbStudies)
	{
		within<-diag(vars[k,])
		invwithin<-diag(1/vars[k,])
		matfixedeffect<-matfixedeffect+invwithin
		mat1<-within+esttrun
		mat1<-qr.solve(mat1,diag(rep(1,length(mat1[,1]))))
		mat<-mat+mat1
		y<-matrix(ys[k,],nrow=NbTimes,ncol=1)
		mat2<-mat1%*%y
		maty<-maty+mat2
		matyfixedeffect<-matyfixedeffect+invwithin%*%y   #solve(within)%*%y
	}
	covrandomeffect<-solve(mat)					## Variance matrix of the pooled arcsine transformed cond. surv. with random effects
	betahatrandomeffect<-covrandomeffect%*%maty     	## Pooled estimates of arcsine transformed cond. surv. with random effects
	covfixedeffect<-solve(matfixedeffect)			## Variance matrix of the pooled arcsine transformed cond. surv. with fixed effects
	betahatfixedeffect<-covfixedeffect%*%matyfixedeffect  ## Pooled estimates of arcsine transformed cond. surv. with fixed effects


#### Heterogeneity (section 2.3. of the paper)

	HeterogeneityQ<-0
	for (i in 1:NbStudies) 
	{
		IndicHeterogeneity<-vars[i,]!=10^12
		HeterogeneityQ<-HeterogeneityQ+t(ys[i,IndicHeterogeneity]-betahatfixedeffect[IndicHeterogeneity])%*%(diag(1/vars[i,]))[IndicHeterogeneity,IndicHeterogeneity]%*%(ys[i,IndicHeterogeneity]-betahatfixedeffect[IndicHeterogeneity])
	}
	HSquared<-HeterogeneityQ/(sum(as.vector(vars)!=10^12)-NbStudies)
	ISquared<-max(0,100*(HSquared-1)/HSquared)  ## in percentage
	HeterogeneityResults<-c(HeterogeneityQ,HSquared,ISquared)


#### Summarized survival probabilities

	PooledSurvivalFE<-cumprod((sin(betahatfixedeffect))^2)   	## Fixed effects
	PooledSurvivalRE<-cumprod((sin(betahatrandomeffect))^2)   	## Random effects


#### 95% confidence intervals around the summarized survival probabilities


	if (confidence=="Greenwood")
	{
		# Fixed effects
		VarLogPooledSurvivalFE<-cumsum(4*(cos(betahatfixedeffect)/sin(betahatfixedeffect))^2*diag(covfixedeffect))
		PooledSurvivalICinfFE<-exp(log(PooledSurvivalFE)-1.96*sqrt(VarLogPooledSurvivalFE))
		PooledSurvivalICsupFE<-exp(log(PooledSurvivalFE)+1.96*sqrt(VarLogPooledSurvivalFE))
		CIPooledSurvivalFE<-cbind(PooledSurvivalICinfFE,PooledSurvivalICsupFE)

		# Random effects
		DerivativeOfTransformation<-cos(betahatrandomeffect)/sin(betahatrandomeffect)
		VarLogPooledConditionalSurvivalProbaRE<-(DerivativeOfTransformation%*%t(DerivativeOfTransformation))*covrandomeffect
		MatriceTriangular<-diag(rep(1,length(betahatrandomeffect)))
		MatriceTriangular[upper.tri(MatriceTriangular)]<-1
		VarLogPooledSurvivalRE<-4*(t(MatriceTriangular)%*%VarLogPooledConditionalSurvivalProbaRE%*%MatriceTriangular)
		PooledSurvivalICinfRE<-exp(log(PooledSurvivalRE)-1.96*sqrt(diag(VarLogPooledSurvivalRE)))
		PooledSurvivalICsupRE<-exp(log(PooledSurvivalRE)+1.96*sqrt(diag(VarLogPooledSurvivalRE)))
		CIPooledSurvivalRE<-cbind(PooledSurvivalICinfRE,PooledSurvivalICsupRE)
}


	SimulatedArcSinCondSurvFE<-rmvnorm(n=10000,mean=betahatfixedeffect,sigma = covfixedeffect)
	SimulatedArcSinCondSurvRE<-rmvnorm(n=10000,mean=betahatrandomeffect,sigma = covrandomeffect)
	SimulatedPooledSurvFE<-apply((sin(SimulatedArcSinCondSurvFE))^2,1,cumprod)
	SimulatedPooledSurvRE<-apply((sin(SimulatedArcSinCondSurvRE))^2,1,cumprod)

	if (confidence=="MonteCarlo")
	{
		# Fixed effects

		CIPooledSurvivalFE<-t(apply(SimulatedPooledSurvFE,1,quantile,probs=c(0.025,0.975)))

		# Random effects

		CIPooledSurvivalRE<-t(apply(SimulatedPooledSurvRE,1,quantile,probs=c(0.025,0.975)))
	}




#### Median survival time from the summarized survival curve


	# Fixed effects
	if (PooledSurvivalFE[length(PooledSurvivalFE)]>0.50) MedianTimeFE<-NA
	if (PooledSurvivalFE[length(PooledSurvivalFE)]==0.50) MedianTimeFE<-IndiceTimes[length(IndiceTimes)]
	if (PooledSurvivalFE[length(PooledSurvivalFE)]<0.50)
	{
		IndexMin<-max((1:length(PooledSurvivalFE))[PooledSurvivalFE>0.5])
		IndexMax<-min((1:length(PooledSurvivalFE))[PooledSurvivalFE<0.5])
		MedianTimeFE<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(PooledSurvivalFE[IndexMax]-0.5)/(PooledSurvivalFE[IndexMax]-PooledSurvivalFE[IndexMin])
	}
	MedianTimeFEbtp<-rep(NA,10000)
	for (i in 1:10000)
	{
		if (SimulatedPooledSurvFE[length(PooledSurvivalFE),i]<0.50)
		{
			IndexMin<-max((1:length(SimulatedPooledSurvFE[,i]))[SimulatedPooledSurvFE[,i]>0.5])
			IndexMax<-min((1:length(SimulatedPooledSurvFE[,i]))[SimulatedPooledSurvFE[,i]<0.5])
			MedianTimeFEbtp[i]<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(SimulatedPooledSurvFE[IndexMax,i]-0.5)/(SimulatedPooledSurvFE[IndexMax,i]-SimulatedPooledSurvFE[IndexMin,i])
		}
	}
	if (sum(is.na(MedianTimeFEbtp))==0) CIMedianTimeFE<-quantile(MedianTimeFEbtp,probs=c(0.025,0.975))
	if (sum(is.na(MedianTimeFEbtp))>0) CIMedianTimeFE<-c(NA,NA)



	# Random effects
	if (PooledSurvivalRE[length(PooledSurvivalRE)]>0.50) MedianTimeRE<-NA
	if (PooledSurvivalRE[length(PooledSurvivalRE)]==0.50) MedianTimeRE<-IndiceTimes[length(IndiceTimes)]
	if (PooledSurvivalRE[length(PooledSurvivalRE)]<0.50)
	{
		IndexMin<-max((1:length(PooledSurvivalRE))[PooledSurvivalRE>0.5])
		IndexMax<-min((1:length(PooledSurvivalRE))[PooledSurvivalRE<0.5])
		MedianTimeRE<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(PooledSurvivalRE[IndexMax]-0.5)/(PooledSurvivalRE[IndexMax]-PooledSurvivalRE[IndexMin])
	}
	MedianTimeREbtp<-rep(NA,10000)
	for (i in 1:10000)
	{
		if (SimulatedPooledSurvRE[length(PooledSurvivalRE),i]<0.50)
		{
			IndexMin<-max((1:length(SimulatedPooledSurvRE[,i]))[SimulatedPooledSurvRE[,i]>0.5])
			IndexMax<-min((1:length(SimulatedPooledSurvRE[,i]))[SimulatedPooledSurvRE[,i]<0.5])
			MedianTimeREbtp[i]<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(SimulatedPooledSurvRE[IndexMax,i]-0.5)/(SimulatedPooledSurvRE[IndexMax,i]-SimulatedPooledSurvRE[IndexMin,i])
		}
	}
	if (sum(is.na(MedianTimeREbtp))==0) CIMedianTimeRE<-quantile(MedianTimeREbtp,probs=c(0.025,0.975))
	if (sum(is.na(MedianTimeREbtp))>0) CIMedianTimeRE<-c(NA,NA)




#### Mean survival time from the summarized survival curve

	DureeIntervalle<-IndiceTimes[-1]-IndiceTimes[-length(IndiceTimes)]

	# MeanTimeFE<-sum(DureeIntervalle*(PooledSurvivalFE[-length(PooledSurvivalFE)]+PooledSurvivalFE[-1])/2)
	MeanTimeFE <- IndiceTimes[1]*(1+PooledSurvivalFE[1])/2+sum(DureeIntervalle * (PooledSurvivalFE[-length(PooledSurvivalFE)] +PooledSurvivalFE[-1])/2)
	
	Temp<-(SimulatedPooledSurvFE[-length(PooledSurvivalFE),]+SimulatedPooledSurvFE[-1,])/2
	for (i in 1:(length(PooledSurvivalFE)-1)) Temp[i,]<-Temp[i,]*DureeIntervalle[i]
	# CIMeanTimeFE<-quantile(apply(Temp,2,sum),probs=c(0.025,0.975))
    CIMeanTimeFE <- quantile(apply(Temp, 2, sum)+IndiceTimes[1]*(SimulatedPooledSurvFE[1,]+rep(1,10000))/2, probs = c(0.025,0.975))
	#MeanTimeRE<-sum(DureeIntervalle*(PooledSurvivalRE[-length(PooledSurvivalRE)]+PooledSurvivalRE[-1])/2)
	MeanTimeRE <- IndiceTimes[1]*(1+PooledSurvivalRE[1])/2+sum(DureeIntervalle * (PooledSurvivalRE[-length(PooledSurvivalRE)] +PooledSurvivalRE[-1])/2)
	
	Temp<-(SimulatedPooledSurvRE[-length(PooledSurvivalRE),]+SimulatedPooledSurvRE[-1,])/2
	for (i in 1:(length(PooledSurvivalRE)-1)) Temp[i,]<-Temp[i,]*DureeIntervalle[i]
	#CIMeanTimeRE<-quantile(apply(Temp,2,sum),probs=c(0.025,0.975))
    CIMeanTimeRE <- quantile(apply(Temp, 2, sum)+IndiceTimes[1]*(SimulatedPooledSurvRE[1,]+rep(1,10000))/2, probs = c(0.025,0.975))

#### Results

SummarySurvivalFE<-cbind(IndiceTimes,PooledSurvivalFE,CIPooledSurvivalFE)
MedianSurvivalTimeFE<-c(MedianTimeFE,CIMedianTimeFE)
MeanSurvivalTimeFE<-c(MeanTimeFE,CIMeanTimeFE)

SummarySurvivalRE<-cbind(IndiceTimes,PooledSurvivalRE,CIPooledSurvivalRE)
MedianSurvivalTimeRE<-c(MedianTimeRE,CIMedianTimeRE)
MeanSurvivalTimeRE<-c(MeanTimeRE,CIMeanTimeRE)

return(list(
verif.data=verif.data,
summary.fixed=SummarySurvivalFE,
median.fixed=MedianSurvivalTimeFE,
mean.fixed=MeanSurvivalTimeFE,
heterogeneity=HeterogeneityResults,
summary.random=SummarySurvivalRE,
median.random=MedianSurvivalTimeRE,
mean.random=MeanSurvivalTimeRE))
}
}







msurv.strata<-function(study, time, n.risk, surv.rate, confidence, strata)

{

.data <- data.frame(study, time, n.risk, surv.rate)
.data <- .data[order(study, time),]
study <- .data$study
time <- .data$time
n.risk <- .data$n.risk
surv.rate <- .data$surv.rate

#### Determination of the time horizon for which each stratum contains at least two studies

# Number of studies
IndiceS<-sort(unique(study))
NbS<-length(IndiceS)
IndiceT<-sort(unique(time))
NbT<-length(IndiceT)
TimeM<-max(IndiceT)


verif <- function(x)
{
ref.study <- sort(time[study==x])
ref.data  <- IndiceT[IndiceT<=max(ref.study)]
return(1*(sum(ref.study!=ref.data)==0))
}

verif.data <- data.frame(study=IndiceS, check=sapply(IndiceS, FUN="verif"))

if (sum(verif.data$check==0)>0) { 
return(list(
verif.data=verif.data,
summary.fixed=NA,
summary.random=NA,
median.fixed=NA,
mean.fixed=NA,
heterogeneity=NA,
p.value=NA))
}

else {

	strata<-factor(strata)

	Temp<-table(time, strata)
	TimeBound<-max(sort(unique(time))[apply(Temp,1,min)>=2])
	Indic<-((time <= TimeBound) & !is.na(strata))

	study<-study[Indic]
	n.risk<-n.risk[Indic]
	surv.rate<-surv.rate[Indic]
	strata<-strata[Indic]
	time<-time[Indic]

#### Calculus of the conditional survivals

	CondSurv<-rep(-99,length(surv.rate))
	for (i in 1:length(surv.rate))
	{
		if (time[i]==min(time)) {CondSurv[i]<-surv.rate[i]}
		if (time[i]!=min(time)) {CondSurv[i]<-surv.rate[i]/surv.rate[i-1]}
	}

#### Number of strata

NbStrata<-length(unique(strata))


SurvivalFE<-list(NULL)
SurvivalRE<-list(NULL)
MedianTime<-data.frame(matrix(,ncol=7,nrow=(NbStrata+1)))
names(MedianTime)<-c("Strata","MedianFE","Lower95%FE","Upper95%FE","MedianRE","Lower95%RE","Upper95%RE")
MeanTime<-data.frame(matrix(,ncol=7,nrow=(NbStrata+1)))
names(MeanTime)<-c("Strata","MeanFE","Lower95%FE","Upper95%FE","MeanRE","Lower95%RE","Upper95%RE")
Heterogeneity<-data.frame(matrix(,ncol=4,nrow=(NbStrata+1)))
names(Heterogeneity)<-c("Strata","Qstatistic","H index","I-squared%")


####  Meta-analysis in each stratum and in all studies

for (g in 1:(NbStrata+1))
{

	if (g<=NbStrata)  # only studies in the stratum g
	{
		IndicStrata<-strata==levels(strata)[g]
		StudyTemp<-study[IndicStrata]
	}
	if (g==(NbStrata+1))  ## all studies
	{
		IndicStrata<-rep(TRUE,length(strata))
		StudyTemp<-study[IndicStrata]
	}

#### Mise en forme des data

	# Arc sine tranformation of the conditional survival (Equations page 7)
	ArcSineCondSurv<-asin(sqrt((n.risk[IndicStrata]*CondSurv[IndicStrata]+0.25)/(n.risk[IndicStrata]+0.5)))
	VarArcSineCondSurv<-1/(4*(n.risk[IndicStrata]+0.5))

	# Number of studies
	IndiceStudies<-sort(unique(study[IndicStrata]))
	NbStudies<-length(IndiceStudies)

	# Number of points in time (number of pooled conditional survival to be estimated)
	IndiceTimes<-sort(unique(time[IndicStrata]))
	NbTimes<-length(IndiceTimes)
	TimeMax<-max(IndiceTimes)

	
	# ys est la matrice des proba conditionnelles arcsine transformées (1 ligne par étude, 1 colonne par temps, valeurs manquantes remplacées par la moyenne des autres)
	ys<-matrix(nrow=NbStudies,ncol=NbTimes)
	for (i in 1:NbStudies) ys[i,1:sum(StudyTemp==IndiceStudies[i])]<-ArcSineCondSurv[StudyTemp==IndiceStudies[i]]
	for (i in 1:NbTimes)  ys[is.na(ys[,i]),i]<-mean(ys[!is.na(ys[,i]),i])

	# vars est la matrice des variances des proba conditionnelles arcsine transformées (1 ligne par étude, 1 colonne par temps, valeurs manquantes remplacées par 10^12)
	vars<-matrix(nrow=NbStudies,ncol=NbTimes)
	for (i in 1:NbStudies) vars[i,1:sum(StudyTemp==IndiceStudies[i])]<-VarArcSineCondSurv[StudyTemp==IndiceStudies[i]]
	for (i in 1:NbTimes)  vars[is.na(vars[,i]),i]<-10^12



#### Application of the multivariate DerSimonian and Laird methodology (R code from Dan Jackson)
#### Equations page 7
#### Estimation of the pooled arc sine transformed conditional survival in each stratum


	mat<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	maty<-matrix(0, nrow=NbTimes, ncol=1)
	matfixedeffect<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	matyfixedeffect<-matrix(0, nrow=NbTimes, ncol=1)
	esttrun<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	estmat<-matrix(0, nrow=NbTimes, ncol=NbTimes)
	count<-1
	
	for(i in 1:(NbTimes-1))
	{
		for(j in (i+1):NbTimes) 
		{
			crossses<-(vars[,i]*vars[,j])^0.5
			mean1a<-sum(ys[,i]/crossses)/sum(1/crossses)
			mean2a<-sum(ys[,j]/crossses)/sum(1/crossses)
			Q<-sum((ys[,i]-mean1a)*(ys[,j]-mean2a)/crossses)
			bb<-sum(1/crossses)-sum(1/(crossses^2))/sum(1/crossses)
			est<-Q/bb
			estmat[i,j]<-est
		}
	}
	estmat<-estmat+t(estmat)
	for(i in 1:NbTimes)
	{
		cross_ses<-vars[,i]
		mean1a<-sum(ys[,i]/vars[,i])/sum(1/vars[,i])
		Q<-sum((ys[,i]-mean1a)*(ys[,i]-mean1a)/vars[,i])
		rhos<-rep(1,NbStudies)
		check<-(1-1*(vars[,i]==10^12))
		rhos<-rhos*check
		aa<-sum(rhos)-sum(rhos/vars[,i])/sum(1/vars[,i])
		bb<-sum(1/vars[,i])-sum(1/(vars[,i]^2))/sum(1/vars[,i])
		est<-(Q-aa)/bb
		estmat[i,i]<-est
	}
	eig<-eigen(estmat)
	for(i in 1:NbTimes) esttrun<-esttrun+max(0, eig$values[i])*eig$vectors[,i]%*%t(eig$vectors[,i])
	for(k in 1:NbStudies)
	{
		within<-diag(vars[k,])
		invwithin<-diag(1/vars[k,])
		matfixedeffect<-matfixedeffect+invwithin 
		mat1<-within+esttrun
		mat1<-qr.solve(mat1,diag(rep(1,length(mat1[,1]))))
		mat<-mat+mat1
		y<-matrix(ys[k,],nrow=NbTimes,ncol=1)
		mat2<-mat1%*%y
		maty<-maty+mat2
		matyfixedeffect<-matyfixedeffect+invwithin %*%y
	}
	covrandomeffect<-solve(mat)					## Variance matrix of the pooled arcsine transformed cond. surv. with random effects
	betahatrandomeffect<-covrandomeffect%*%maty     	## Pooled estimates of arcsine transformed cond. surv. with random effects
	covfixedeffect<-solve(matfixedeffect)			## Variance matrix of the pooled arcsine transformed cond. surv. with fixed effects
	betahatfixedeffect<-covfixedeffect%*%matyfixedeffect  ## Pooled estimates of arcsine transformed cond. surv. with fixed effects


#### Heterogeneity (section 2.3. of the paper)

	HeterogeneityQ<-0
	for (i in 1:NbStudies) 
	{
		IndicHeterogeneity<-vars[i,]!=10^12
		HeterogeneityQ<-HeterogeneityQ+t(ys[i,IndicHeterogeneity]-betahatfixedeffect[IndicHeterogeneity])%*%(diag(1/vars[i,]))[IndicHeterogeneity,IndicHeterogeneity]%*%(ys[i,IndicHeterogeneity]-betahatfixedeffect[IndicHeterogeneity])
	}
	HSquared<-HeterogeneityQ/(sum(as.vector(vars)!=10^12)-NbStudies)
	ISquared<-max(0,100*(HSquared-1)/HSquared)  ## in percentage
	HeterogeneityResults<-c(HeterogeneityQ,HSquared,ISquared)



#### Summarized survival probabilities

	PooledSurvivalFE<-cumprod((sin(betahatfixedeffect))^2)   	## Fixed effects
	PooledSurvivalRE<-cumprod((sin(betahatrandomeffect))^2)   	## Random effects


#### 95% confidence intervals around the summarized survival probabilities


	if (confidence=="Greenwood")
	{
		# Fixed effects
		VarLogPooledSurvivalFE<-cumsum(4*(cos(betahatfixedeffect)/sin(betahatfixedeffect))^2*diag(covfixedeffect))
		PooledSurvivalICinfFE<-exp(log(PooledSurvivalFE)-1.96*sqrt(VarLogPooledSurvivalFE))
		PooledSurvivalICsupFE<-exp(log(PooledSurvivalFE)+1.96*sqrt(VarLogPooledSurvivalFE))
		CIPooledSurvivalFE<-cbind(PooledSurvivalICinfFE,PooledSurvivalICsupFE)

		# Random effects
		DerivativeOfTransformation<-cos(betahatrandomeffect)/sin(betahatrandomeffect)
		VarLogPooledConditionalSurvivalProbaRE<-(DerivativeOfTransformation%*%t(DerivativeOfTransformation))*covrandomeffect
		MatriceTriangular<-diag(rep(1,length(betahatrandomeffect)))
		MatriceTriangular[upper.tri(MatriceTriangular)]<-1
		VarLogPooledSurvivalRE<-4*(t(MatriceTriangular)%*%VarLogPooledConditionalSurvivalProbaRE%*%MatriceTriangular)
		PooledSurvivalICinfRE<-exp(log(PooledSurvivalRE)-1.96*sqrt(diag(VarLogPooledSurvivalRE)))
		PooledSurvivalICsupRE<-exp(log(PooledSurvivalRE)+1.96*sqrt(diag(VarLogPooledSurvivalRE)))
		CIPooledSurvivalRE<-cbind(PooledSurvivalICinfRE,PooledSurvivalICsupRE)
}


	SimulatedArcSinCondSurvFE<-rmvnorm(n=10000,mean=betahatfixedeffect,sigma = covfixedeffect)
	SimulatedArcSinCondSurvRE<-rmvnorm(n=10000,mean=betahatrandomeffect,sigma = covrandomeffect)
	SimulatedPooledSurvFE<-apply((sin(SimulatedArcSinCondSurvFE))^2,1,cumprod)
	SimulatedPooledSurvRE<-apply((sin(SimulatedArcSinCondSurvRE))^2,1,cumprod)

	if (confidence=="MonteCarlo")
	{
		# Fixed effects

		CIPooledSurvivalFE<-t(apply(SimulatedPooledSurvFE,1,quantile,probs=c(0.025,0.975)))

		# Random effects

		CIPooledSurvivalRE<-t(apply(SimulatedPooledSurvRE,1,quantile,probs=c(0.025,0.975)))
	}




#### Median survival time from the summarized survival curve


	# Fixed effects
	if (PooledSurvivalFE[length(PooledSurvivalFE)]>0.50) MedianTimeFE<-NA
	if (PooledSurvivalFE[length(PooledSurvivalFE)]==0.50) MedianTimeFE<-IndiceTimes[length(IndiceTimes)]
	if (PooledSurvivalFE[length(PooledSurvivalFE)]<0.50)
	{
		IndexMin<-max((1:length(PooledSurvivalFE))[PooledSurvivalFE>0.5])
		IndexMax<-min((1:length(PooledSurvivalFE))[PooledSurvivalFE<0.5])
		MedianTimeFE<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(PooledSurvivalFE[IndexMax]-0.5)/(PooledSurvivalFE[IndexMax]-PooledSurvivalFE[IndexMin])
	}
	MedianTimeFEbtp<-rep(NA,10000)
	for (i in 1:10000)
	{
		if (SimulatedPooledSurvFE[length(PooledSurvivalFE),i]<0.50)
		{
			IndexMin<-max((1:length(SimulatedPooledSurvFE[,i]))[SimulatedPooledSurvFE[,i]>0.5])
			IndexMax<-min((1:length(SimulatedPooledSurvFE[,i]))[SimulatedPooledSurvFE[,i]<0.5])
			MedianTimeFEbtp[i]<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(SimulatedPooledSurvFE[IndexMax,i]-0.5)/(SimulatedPooledSurvFE[IndexMax,i]-SimulatedPooledSurvFE[IndexMin,i])
		}
	}
	if (sum(is.na(MedianTimeFEbtp))==0) CIMedianTimeFE<-quantile(MedianTimeFEbtp,probs=c(0.025,0.975))
	if (sum(is.na(MedianTimeFEbtp))>0) CIMedianTimeFE<-c(NA,NA)



	# Random effects
	if (PooledSurvivalRE[length(PooledSurvivalRE)]>0.50) MedianTimeRE<-NA
	if (PooledSurvivalRE[length(PooledSurvivalRE)]==0.50) MedianTimeRE<-IndiceTimes[length(IndiceTimes)]
	if (PooledSurvivalRE[length(PooledSurvivalRE)]<0.50)
	{
		IndexMin<-max((1:length(PooledSurvivalRE))[PooledSurvivalRE>0.5])
		IndexMax<-min((1:length(PooledSurvivalRE))[PooledSurvivalRE<0.5])
		MedianTimeRE<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(PooledSurvivalRE[IndexMax]-0.5)/(PooledSurvivalRE[IndexMax]-PooledSurvivalRE[IndexMin])
	}
	MedianTimeREbtp<-rep(NA,10000)
	for (i in 1:10000)
	{
		if (SimulatedPooledSurvRE[length(PooledSurvivalRE),i]<0.50)
		{
			IndexMin<-max((1:length(SimulatedPooledSurvRE[,i]))[SimulatedPooledSurvRE[,i]>0.5])
			IndexMax<-min((1:length(SimulatedPooledSurvRE[,i]))[SimulatedPooledSurvRE[,i]<0.5])
			MedianTimeREbtp[i]<-IndiceTimes[IndexMax]-(IndiceTimes[IndexMax]-IndiceTimes[IndexMin])*(SimulatedPooledSurvRE[IndexMax,i]-0.5)/(SimulatedPooledSurvRE[IndexMax,i]-SimulatedPooledSurvRE[IndexMin,i])
		}
	}
	if (sum(is.na(MedianTimeREbtp))==0) CIMedianTimeRE<-quantile(MedianTimeREbtp,probs=c(0.025,0.975))
	if (sum(is.na(MedianTimeREbtp))>0) CIMedianTimeRE<-c(NA,NA)




#### Mean survival time from the summarized survival curve

	DureeIntervalle<-IndiceTimes[-1]-IndiceTimes[-length(IndiceTimes)]

	#MeanTimeFE<-sum(DureeIntervalle*(PooledSurvivalFE[-length(PooledSurvivalFE)]+PooledSurvivalFE[-1])/2)
	
	MeanTimeFE <- IndiceTimes[1]*(1+PooledSurvivalFE[1])/2+sum(DureeIntervalle * (PooledSurvivalFE[-length(PooledSurvivalFE)] +PooledSurvivalFE[-1])/2)
	
	Temp<-(SimulatedPooledSurvFE[-length(PooledSurvivalFE),]+SimulatedPooledSurvFE[-1,])/2
	for (i in 1:(length(PooledSurvivalFE)-1)) Temp[i,]<-Temp[i,]*DureeIntervalle[i]
	#CIMeanTimeFE<-quantile(apply(Temp,2,sum),probs=c(0.025,0.975))
    CIMeanTimeFE <- quantile(apply(Temp, 2, sum)+IndiceTimes[1]*(SimulatedPooledSurvFE[1,]+rep(1,10000))/2, probs = c(0.025,0.975))
	
	#MeanTimeRE<-sum(DureeIntervalle*(PooledSurvivalRE[-length(PooledSurvivalRE)]+PooledSurvivalRE[-1])/2)
	
	MeanTimeRE <- IndiceTimes[1]*(1+PooledSurvivalRE[1])/2+sum(DureeIntervalle * (PooledSurvivalRE[-length(PooledSurvivalRE)] +PooledSurvivalRE[-1])/2)
	
	MeanTimeRE <- IndiceTimes[1]*(1+PooledSurvivalRE[1])/2+sum(DureeIntervalle * (PooledSurvivalRE[-length(PooledSurvivalRE)] +PooledSurvivalRE[-1])/2)
	
	Temp<-(SimulatedPooledSurvRE[-length(PooledSurvivalRE),]+SimulatedPooledSurvRE[-1,])/2
	for (i in 1:(length(PooledSurvivalRE)-1)) Temp[i,]<-Temp[i,]*DureeIntervalle[i]
	# CIMeanTimeRE<-quantile(apply(Temp,2,sum),probs=c(0.025,0.975))
    CIMeanTimeRE <- quantile(apply(Temp, 2, sum)+IndiceTimes[1]*(SimulatedPooledSurvRE[1,]+rep(1,10000))/2, probs = c(0.025,0.975))

#### Results per stratum

	SummarySurvivalFE<-data.frame(cbind(IndiceTimes,round(PooledSurvivalFE,3),round(CIPooledSurvivalFE,3)))
	names(SummarySurvivalFE)<-c("time","SurvivalFE","Lower CI95%","Upper CI95%")

	MedianSurvivalTimeFE<-round(c(MedianTimeFE,CIMedianTimeFE),2)
	MeanSurvivalTimeFE<-round(c(MeanTimeFE,CIMeanTimeFE),2)

	SummarySurvivalRE<-data.frame(cbind(IndiceTimes,round(PooledSurvivalRE,3),round(CIPooledSurvivalRE,3)))
	names(SummarySurvivalRE)<-c("time","SurvivalRE","Lower CI95%","Upper CI95%")

	MedianSurvivalTimeRE<-round(c(MedianTimeRE,CIMedianTimeRE),2)
	MeanSurvivalTimeRE<-round(c(MeanTimeRE,CIMeanTimeRE),2)


#### Put results of each stratum in the global results

	SurvivalFE[[g]]<-SummarySurvivalFE
	SurvivalRE[[g]]<-SummarySurvivalRE
	if (g<=NbStrata)
	{
		MedianTime[g,]<-c(levels(strata)[g],MedianSurvivalTimeFE,MedianSurvivalTimeRE)
		MeanTime[g,]<-c(levels(strata)[g],MeanSurvivalTimeFE,MeanSurvivalTimeRE)
		Heterogeneity[g,]<-c(levels(strata)[g],round(HeterogeneityResults,2))
	}
	if (g==(NbStrata+1))
	{
		MedianTime[g,]<-c("All strata",MedianSurvivalTimeFE,MedianSurvivalTimeRE)
		MeanTime[g,]<-c("All strata",MeanSurvivalTimeFE,MeanSurvivalTimeRE)
		Heterogeneity[g,]<-c("All strata",round(HeterogeneityResults,2))
	}
}

names(SurvivalFE)<-c(levels(strata),"All strata")
names(SurvivalRE)<-c(levels(strata),"All strata")


BetweenStudyQStatistic<-as.numeric(Heterogeneity$Qstatistic[NbStrata+1])-sum(as.numeric(Heterogeneity$Qstatistic[1:NbStrata]))
DegreeOfFreedom<-NbTimes*(NbStrata-1)
Pvalue<-1-pchisq(BetweenStudyQStatistic,df=DegreeOfFreedom)

return(list(
verif.data=verif.data,
summary.fixed=SurvivalFE,
summary.random=SurvivalRE,
median.fixed=MedianTime,
mean.fixed=MeanTime,
heterogeneity=Heterogeneity,
p.value=Pvalue))

}
}