Hi, May be this helps you.
#Using set.seed(12345) S=10 simdata <- replicate(S, generate(250)) lstpshat<-lapply(seq_len(ncol(simdata)),function(i) {glm.t<-glm(t~x1+x2+x3+x4+x5+x6+x7+I(x2^2)+I(x4^2)+I(x7^2)+x1:x3+x2:x4+x3:x5+x4:x6+x5:x7+x1:x6+x2:x3+x3:x4+x4:x5+x5:x6,family=binomial,data=simdata[,i]); pshat<- predict(glm.t,type="response")}) simdata1<-rbind(simdata,pshat=lstpshat) pdf("hist1.pdf") lapply(seq_len(ncol(simdata1)),function(i){ x1<- simdata1[,i]; pshat0<-x1$pshat[x1$t==0];pshat1<- x1$pshat[x1$t==1]; hist(pshat1,xlim=c(0,1),col=rgb(0.7,0,0,0.5)); hist(pshat0,add=T,col=rgb(0,0,1,0.3))}) dev.off() (You need to change the colors as per your requirements) A.K. >thanks! it really helps! > >anyway, how will i have a histogram of the lstpshat basing on the value of t, that is, fot t=1, the color is red, for t=0, the color is blue and for their >overlap, the color is green? thanks a lot! > >set.seed(12345) >S=1000 >generate <- function(size) { >x1 <- rnorm(size, mean=0, sd=1) >x2 <- rnorm(size, mean=0, sd=1) >x3 <- rnorm(size, mean=0, sd=1) >x4 <- rnorm(size, mean=0, sd=1) >x5 <- rnorm(size, mean=0, sd=1) >x6 <- rnorm(size, mean=0, sd=1) >x7 <- rnorm(size, mean=0, sd=1) >x8 <- rnorm(size, mean=0, sd=1) >x9 <- rnorm(size, mean=0, sd=1) >x10 <- rnorm(size, mean=0, sd=1) >e<-rnorm(size, mean=0, sd=1) >t_trueps <- (1 + exp( -(b0 + b1*x1 + b2*x2 + b3*x3 + b4*x4 + b5*x5 + b6*x6 + >b7*x7 >+ b2*x2*x2 + b4*x4*x4 + b7*x7*x7 + b1*0.5*x1*x3 + b2*0.7*x2*x4 +b3*0.5*x3*x5 >+ b4*0.7*x4*x6 + b5*0.5*x5*x7 + b1*0.5*x1*x6 + b2*0.7*x2*x3 + b3*0.5*x3*x4 >+ b4*0.5*x4*x5 + b5*0.5*x5*x6) ) )^-1 >prob.exposure <- runif(size) >t <- ifelse(t_trueps > prob.exposure, 1, 0) >y <- a0 + a1*x1 + a2*x2 + a3*x3 + a4*x4 +a5*x8 + a6*x9 + a7*x10 + g1*t + e >sim <- as.data.frame(cbind(x1, x2, x3 ,x4, x5, x6, x7, x8, x9, x10, t, y)) >return(sim) >} >b0 <- 0.05 >b1 <- 0.95 >b2 <- -0.25 >b3 <- 0.6 >b4 <- -0.4 >b5 <- -0.8 >b6 <- -0.5 >b7 <- 0.7 >a0 <- -3.85 >a1 <- 0.3 >a2 <- -0.36 >a3 <- -0.73 >a4 <- -0.2 >a5 <- 0.71 >a6 <- -0.19 >a7 <- 0.26 >g1 <- -0.4 >simdata <- replicate(S, generate(3000)) > >lstpshat<-lapply(seq_len(ncol(simdata)),function(i) >{glm.t<-glm(t~x1+x2+x3+x4+x5+x6+x7+I(x2^2)+I(x4^2)+I(x7^2)+x1:x3+x2:x4+x3:x5+x4:x6+x5:x7+x1:x6+x2:x3+x3:x4+x4:x5+x5:x6,family=binomial,data=simdata[,i]); > > >pshat<- predict(glm.t,type="response")}) >simdata1<-rbind(simdata,pshat=lstpshat) > >simdata.ps1<- simdata1 >simdata.ps1[]<-do.call(c,lapply(seq_len(ncol(simdata1)),function(i) lapply(simdata1[,i],function(x) x[simdata1[,i]$t==1]))) >lstm1<- lapply(seq_len(ncol(simdata.ps1)),function(i) {dat<-do.call(data.frame,lapply(simdata.ps1[,i],function(x) x));if(nrow(dat)!=0) >{glm.1<-glm(y~x1+x2+x3+x4+x8+x9+x10,data=dat)} else NULL; glm.1; m1<- predict(glm.1)}) > >simdata.ps0<- simdata1 >simdata.ps0[]<-do.call(c,lapply(seq_len(ncol(simdata1)),function(i) lapply(simdata1[,i],function(x) x[simdata1[,i]$t==0]))) >lstm0<-lapply(seq_len(ncol(simdata.ps0)),function(i) {dat<-do.call(data.frame,lapply(simdata.ps0[,i],function(x) x));if(nrow(dat)!=0) >{glm.0<-glm(y~x1+x2+x3+x4+x8+x9+x10,data=dat)} else NULL; glm.0; m0<- predict(glm.0)}) > >simdata.psm1<- rbind(simdata.ps1,m1=lstm1) >simdata.psm0<- rbind(simdata.ps0,m0=lstm0) >
hist1.pdf
Description: Adobe PDF document
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.