Pareto.TypeI.simu <- function(G,N,xi1,xi2,shape1=1,shape2=1,beta1,beta2, eta=0.5,theta=2,alpha=1,beta12=0, C.min=0,C.max=5,Z.dist=runif,...){ X.vec=D.vec=C.vec=t.event=t.death=event=death=Z=group=NULL ij=0 for(i in 1:G){ u=rgamma(1,shape=1/eta,scale=eta) for(j in 1:N){ ij=ij+1 group[ij]=i Z[ij]=Z.dist(1,...) theta12=theta*exp(beta12*Z[ij]) r1=(-shape1*u*exp(beta1*Z[ij]))^-1 r2=(theta12*shape2*(u^alpha)*exp(beta2*Z[ij]))^-1 V1=runif(1) V2=runif(1) X=xi1*(1-V1)^r1 W=(1-V1)^(-theta12) D=xi2*(1-W+W*(1-V2)^(-theta12/(theta12+1)))^r2 C=runif(1,min=C.min,max=C.max) X.vec[ij]=X D.vec[ij]=D C.vec[ij]=C t.event[ij]=min(X,D,C) t.death[ij]=min(D,C) event[ij]=as.numeric( t.event[ij]==X ) death[ij]=as.numeric( t.death[ij]==D ) } } data.frame(X=X.vec,D=D.vec,C=C.vec, t.event=t.event,event=event, t.death=t.death,death=death, group=group,Z=Z) }