Nothing
F1addfg<-function(t,lam0=0.5,beta=c(-0.5,-0.005,-0.004),x=rep(0,3)) # FG
{ ## {{{
baset <- 0.13*pnorm((t-.70)/.13)
xm <- matrix(x,ncol=3)
return( 1 - exp(-baset*exp(xm %*% matrix(beta,3,1))))
} ## }}}
##' @export
corsim.prostate <- function(n,theta=1,thetaslope=0,crate=2,test=0,pcens=0,mt=1,same.cens=TRUE,country=TRUE)
{ ## {{{
###n <- 10; theta <- 1; thetaslope <- 0; mt <- 1
if (country==TRUE) xl <- sample(1:4,n,replace=TRUE) else xl <- rep(1,n)
x<-cbind(xl==2,xl==3,xl==4)*1
tt<-seq(0,mt,length=mt*100)
###
###n=100;theta=1;lam0=0.5;beta=0.3;crate=2
thetat <- exp(log(theta))
F11x<-F1addfg(mt,x=x)
F12x<-F1addfg(mt,x=x)
###
thetaslut <- exp(log(theta)+thetaslope*(mt-mt/2))
p11 <- thetaslut*F11x*F12x/((1-F11x)+thetaslut*F11x)
p12 <- F11x-p11
p21 <- F12x-p11
p22 <- 1-F12x-F11x+p11
###apply(cbind(p11,p12,p21),1,sum)
if (test==1) { ## {{{
for (i in 1:2) {
print(x[i,]);
F11xt<-F1addfg(tt,x=x[i,])
F12xt<-F1addfg(tt,x=x[i,])
p11t <- thetat* F11xt*F12xt/((1-F11xt)+thetat*F11xt)
cortt <- ((p11t)/(F12xt-p11t))/(F11xt/(1-F11xt))
###plot(tt,log(cortt))
if (i==1) {
plot(tt,p11t,type="l",ylim=c(0,0.1),xlim=c(0,mt))
###lines(tt,F11x[i]-p11t,col=2)
###lines(tt,F12x[i]-p11t,col=2)
} else lines(tt,p11t,col=2);
###if (sum(diff(p11t<0))>0) stop("dec\n");
###p11 <- max(p11t)
###p12 <- F11x[i]-p11
###p21 <- F12x[i]-p11
###p22 <- 1- F12x[i]-F11x[i]+p11
###pnn <- 1- F12x[i]-F11x[i]+p11
}
} ## }}}
###apply(cbind(p11,p12,p21,p22),1,sum)
###
types <- rep(0,n)
causes <- matrix(0,n,2)
stime<-matrix(mt+1,n,2);
for (i in 1:n)
{
ptype <- runif(1)
if (ptype<=p11[i]) {
types[i] <- 1
myhazx<-F1addfg(tt,x=x[i,])/F12x[i]
### if (abs(max(myhazx)-1)> 0.001) stop("not dist\n");
stime[i,2]<-Cpred(cbind(myhazx,tt),runif(1))[1,2]+runif(1,0,0.001)
f1<- F1addfg(tt,x=x[i,])
myhazx<- (F12x[i]/p11[i]) * (thetat*f1/((1-f1)+thetat*f1))
### if (abs(max(myhazx)-1)> 0.001) stop("not dist\n");
stime[i,1]<-Cpred(cbind(myhazx,tt),runif(1))[1,2]+runif(1,0,0.001)
causes[i,] <- c(1,1)
}
if ((ptype>p11[i]) & (ptype<=p12[i]+p11[i])) {
types[i] <- 2
f1 <- F1addfg(tt,x=x[i,])
myhazx<- ( f1 - thetat*F12x[i]*f1/((1-f1)+thetat*f1))/p12[i];
myhazx <- f1/F11x[i]
### if (abs(max(myhazx)-1)> 0.001) stop("not dist 2 \n");
stime[i,1]<-Cpred(cbind(myhazx,tt),runif(1))[1,2]+runif(1,0,0.001)
causes[i,] <- c(1,2)
stime[i,2] <- runif(1)*mt
}
if ((ptype>p11[i]+p12[i]) && (ptype<=p21[i]+p12[i]+p11[i])) {
types[i] <- 3
f2 <- F1addfg(tt,x=x[i,])
myhazx <- (f2 - (thetat*F11x[i]*f2/((1-F11x[i])+thetat*F11x[i])))/p21[i];
myhazx <- f2/F12x[i]
### if (abs(max(myhazx)-1)> 0.001) stop("not dist3 \n");
stime[i,2]<-Cpred(cbind(myhazx,tt),runif(1))[1,2]+runif(1,0,0.001)
causes[i,] <- c(2,1)
stime[i,1] <- runif(1)*mt
}
if (ptype>p11[i]+p12[i]+p21[i] ) {
types[i] <- 4
causes[i,] <- c(2,2)
stime[i,1:2] <- runif(2)*mt
}
}
###stime
###causes
stime <- c(t(stime))
cause <- c(t(causes))
###same.cens=TRUE
if (same.cens==TRUE) {
ctime <- rep(rbinom(n,1,pcens),each=2)
ctime[ctime==1] <- rep(runif(sum(ctime==1)/2),each=2)*crate
}
else {
ctime<- rbinom(n,1,pcens)
ctime[ctime==1] <- runif(sum(ctime==1))*crate
}
ctime[ctime==0] <- mt;
cens <- (ctime< stime)
time <- ifelse(cens,ctime,stime)
cause <- ifelse(cens,0,cause)
id <- rep(1:n,rep(2,n))
country <- c()
country[xl==1] <- "SWE"
country[xl==2] <- "DK"
country[xl==3] <- "FIN"
country[xl==4] <- "NOR"
data<-data.frame(time=time,cause=cause,xl=rep(xl,each=2),
country=rep(country,each=2),id=id,cens=cens,stime=stime,type=rep(types,each=2),
f1=rep(F11x,each=2),p11=rep(p11,each=2),p12=rep(p12,each=2),p21=rep(p21,each=2),
p22=rep(p22,each=2))
return(data)
} ## }}}
##' @export
simnordic <- function(n,cordz=2,cormz=3,cratemz=2,cratedz=2,pcensmz=0.8,pcensdz=0.8,country=TRUE,same.cens=TRUE)
{ ## {{{
outdz <- corsim.prostate(n,theta=cordz,crate=cratedz,pcens=pcensdz,mt=1,same.cens=same.cens,test=0,country=country)
outmz <- corsim.prostate(n,theta=cormz,crate=cratemz,pcens=pcensmz,mt=1,same.cens=same.cens,test=0,country=country)
outdz$zyg <- "DZ"
outmz$zyg <- "MZ"
outmz$id <- outmz$id+nrow(outdz)
###
out <- rbind(outdz,outmz)
out$time <- out$time*100
###table(out$type,out$country)
###table(out$type,out$cause)
if (country==TRUE) out$country <- relevel(factor(out$country),ref="SWE")
###table(out$country)
outk <- out[,c("country","cause","id","time","zyg","type")]
###table(outk$cause)
###table(outk$type,outk$country)
###table(outk$cause,outk$country)
###
return(outk)
} ## }}}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.