while(min(stickers)==0){
ico<-ico+1
ival<-sample(1:N,size=5,rep=T)
stickers[ival]<-1
}
ipack[irep]<-ico
if(irep %% 100 == 0) print(c(irep,ico))
}
mean(ipack)
N*log(N)
N*ln(N)
N*log(N)
N*log(N)/5
N<-639
nreps<-10000
ipack<-rep(0,nreps)
for(irep in 1:nreps){
stickers<-rep(0,N)
ico<-0
while(min(stickers)==0){
ico<-ico+1
ival<-sample(1:N,size=5,rep=T)
stickers[ival]<-1
}
ipack[irep]<-ico
if(irep %% 100 == 0) print(c(irep,ico))
}
hisr(ipack)
hist(ipack)
mean(ipack)
q()
library(DPpackage)
update.packages(ask='graphics',checkBuilt=TRUE)
library(DPpackage)
>PT
??PT
# Prostate cancer data example
data(nodal)
attach(nodal)
lacid<-log(acid)
# Initial state
state <- NULL
# MCMC parameters
nburn<-20000
nsave<-10000
nskip<-10
ndisplay<-100
mcmc <- list(nburn=nburn,nsave=nsave,
nskip=nskip,ndisplay=ndisplay,
tune1=1.1,tune2=1.1)
# Prior distribution
prior <- list(alpha=1, beta0=c(0,rep(0.75,5)),
Sbeta0=diag(c(100,rep(25,5)),6),M=5)
# Fitting the Finite Polya tree model
fit1 <- FPTbinary(ssln~age+lacid+xray+size+grade,
prior=prior,mcmc=mcmc,
state=state,status=TRUE)
fit1
# Summary with HPD and Credibility intervals
summary(fit1)
plot(fit1)
plot(fit1,ask=FALSE,nfigr=1,nfigc=2,param="xray")
plot(fit1,ask=FALSE,param="link",nfigc=1,nfigr=1)
anova(fit1)
??Bernstein
q()
Sigma<-rWishart(1,p,diag(rep(1,p))/p)
Sigma
PCF.mat<-array(0,c(nreps,p,p))
library(MASS)
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
nreps<-1000
p<-10
n<-100
Sigma<-rWishart(1,p,diag(rep(1,p))/p)[,,1]
PCF.mat<-array(0,c(nreps,p,p))
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
X
cor(X)
Cor.mat<-PCF.mat<-array(0,c(nreps,p,p))
Cor.mat[irep,,]<-cor(X)
??partial
?cor
install.packages('ggm')
pcor(X)
library(ggm)
pcor(X)
?pcor
parcor(X)
parcor(Cor.mat[irep,,])
nreps<-1000
p<-10
n<-100
Sigma<-rWishart(1,p,diag(rep(1,p))/p)[,,1]
Cor.mat<-PCF.mat<-array(0,c(nreps,p,p))
library(MASS)
library(ggm)
for(irep in 1:nreps){
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
Cor.mat[irep,,]<-cor(X)
PCF.mat[irep,,]<-parcor(Cor.mat[irep,,])
}
hist(PCF.mat[,1,1])
nreps
parcor(Cor.mat[irep,,])
hist(PCF.mat[,1,2])
diag(PCF.mat)<-NA
hist(PCF.mat[,1,3])
hist(PCF.mat[,1,2:3])
plot(PCF.mat[,1,2:3])
nreps<-10000
p<-10
n<-100
Sigma<-rWishart(1,p,diag(rep(1,p))/p)[,,1]
Cor.mat<-PCF.mat<-array(0,c(nreps,p,p))
library(MASS)
library(ggm)
for(irep in 1:nreps){
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
Cor.mat[irep,,]<-cor(X)
PCF.mat[irep,,]<-parcor(Cor.mat[irep,,])
}
apply(PCF.mat,2:3,range)
apply(PCF.mat,2:3,min)
apply(PCF.mat,2:3,max)
hist(PCF.mat[,8,5])
q()
1150.26+226.70-43.38
182.32- 119.47
21+31+3
21+31+30+8
(55/90)*86.52 + 0.21
22+30+31+8
38*(62.85/91)
62.85 + 95.43 + 53.08 + 26.25
62.85/91
1333.58 + 237.61
1571.19 - 1350
q()
install.packages('data.table');library(data.table)
install.packages('plyr');
library(plyr)
install.packages('data.table');library(data.table)
q()
t<-seq(-10,10,by=0.1)
rm(t)
tv<-seq(-10,10,by=0.1)
i
?Mod
cos(tv)+1i*sin(tv)
Mod(cos(tv)+1i*sin(tv))
q()
?sink
?capture.output
?file
q()
hist(x<-rgamma(10000,0.1,0.1))
hist(x<-rgamma(10000,0.1,0.1),nclass=50)
qgamma(c(1:100)/100,0.1,0.1)
cbind(c(1:100)/100,qgamma(c(1:100)/100,0.1,0.1)
)
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
quantile(x<-1/rgamma(10000,0.1,0.1),prob=c(0.25,.5,.75))
qgamma(c(0.25,.5,.75),0.1,0.1)
1/qgamma(c(0.25,.5,.75),0.1,0.1)
1/qgamma(c(0.25,.5,.75),0.5,0.5)
quantile(x<-1/rgamma(10000,0.5,0.5),prob=c(0.25,.5,.75))
hist(x<-1/rgamma(10000,0.5,0.5),nclass=50)
hist(1/rgamma(10000,0.5,0.5),nclass=50)
hist(1/rgamma(10000,1,1),nclass=50)
hist(1/rgamma(10000,10,10),nclass=50)
hist(1/rgamma(10000,0.1,0.1),nclass=50)
hist(1/rgamma(10000,0.5,0.5),nclass=50)
hist(1/rgamma(10000,1,1),nclass=50)
hist(1/rgamma(10000,2,2),nclass=50)
q()
library(survival)
options(na.action=na.exclude) # retain NA in predictions
fit <- coxph(Surv(time, status) ~ age + ph.ecog + strata(inst), lung)
#lung data set has status coded as 1/2
mresid <- (lung$status-1) - predict(fit, type='expected') #Martingale resid 
predict(fit,type="lp")
predict(fit,type="expected")
predict(fit,type="risk",se.fit=TRUE)
predict(fit,type="terms",se.fit=TRUE)
lung$time
predict(fit,type="expected")
plot(lung$time,predict(fit,type="expected"))
plot(fit)
fit
plot.survfit(fit)
?coxph
plot(survfit(fit))
survfit(fit)
q()
update.packages(ask='graphics',checkBuilt=TRUE)
q()
filecode<-'Set1'
fname<-paste('PS-DR-',filecode','.pdf',sep='')
fname<-paste('PS-DR-',filecode,'.pdf',sep='')
fname
fname<-paste(filecode,'-PS-DR.pdf',sep='')
fname
p1<-glm(Z~X,family=binomial)
ps.vals<-fitted(p1)
psr.fit1<-lm(Y~Z+p1)
#Predictors
X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)
#Data simulation
Xmat<-cbind(rep(1,n),X)
eta.Z<-Xmat %*% (c(1,al.select)*al)
p.Z<-expit(eta.Z)
Z<-rbinom(n,1,p.Z)
eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)
Y<-eta.Y+rnorm(n)*sig.Y
source('GenerateSetUp.R')
filecode<-'Set1'
source('GenerateSetUp.R')
filecode<-'Set1'
source("BasicSim.R") #Stratification
source('GenerateSetUp.R')
rm(list=ls())
source('GenerateSetUp.R')
filecode<-'Set1'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
psr.fit1<-lm(Y~Z+ps.vals)
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
ls()
?save
save(ests.mat,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set1.RData')
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set1.RData')
elapsed.time()
e.time()
??timin
?proc.time
proc.time()
t0<-proc.time()[3]
t0
t0<-proc.time()[3]
source('GenerateSetUp.R')
filecode<-'Set1'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set1.RData')
t1<-proc.time()[3]
print(t1-t0)
source('GenerateSetUp2.R')
filecode<-'Set2'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set2.RData')
t2<-proc.time()[3]
source('GenerateSetUp3.R')
filecode<-'Set3'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set3.RData')
t3<-proc.time()[3]
library(MASS)
set.seed(4323)
#Expit transform
expit<-function(x){return(1/(1+exp(-x)))}
n<-1000
p<-10
x.index<-1:p
#Two spurious
#Two instruments
#Three confounders
#Three predictors of outcome only 
al.select<-c(0,0,1,1,1,1,1,0,0,0)
be.select<-c(0,0,0,0,1,1,1,1,1,1)
#Confounders
confound.list<-x.index[al.select == 1 & be.select==1]
#Instruments
intrument.list<-x.index[al.select ==1 & be.select==0]
#Predictors of outcome only
outcomeonly.list<-x.index[al.select ==0  & be.select==1]
#Spurious variables
spurious.list<-x.index[al.select ==0  & be.select==0]
#Predictors of outcome
outcome.list<-x.index[be.select ==1]
#Predictors of exposure
treatment.list<-x.index[al.select ==1]
################################################
#Generate alpha parameters (for treatment allocation)
al<-rnorm(p+1,0,1.0)
#Generate beta parameters (for outcome regression)
be<-rnorm(p+1,0,3.0)
theta<-2.0
sig.Y<-1
p0<-2
p1<-2
p2<-3
p3<-3
Sig0<-rWishart(1,p0+1,diag(p0)/(p0+1))[,,1]
Sig1<-rWishart(1,p1+1,diag(p1)/(p1+1))[,,1]
Sig2<-rWishart(1,p2+1,diag(p2)/(p2+1))[,,1]
Sig3<-rWishart(1,p3+1,diag(p3)/(p3+1))[,,1]
Sig<-matrix(0,p,p)
Sig[1:p0,1:p0]<-Sig0
Sig[(p0+1):(p0+p1),(p0+1):(p0+p1)]<-Sig1
Sig[(p0+p1+1):(p0+p1+p2),(p0+p1+1):(p0+p1+p2)]<-Sig2
Sig[(p0+p1+p2+1):(p0+p1+p2+p3),(p0+p1+p2+1):(p0+p1+p2+p3)]<-Sig3
#Predictors
X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)
#Data simulation
Xmat<-cbind(rep(1,n),X)
eta.Z<-Xmat %*% (c(1,al.select)*al)
p.Z<-expit(eta.Z)
Z<-rbinom(n,1,p.Z)
eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)
hist(p.Z)
library(MASS)
set.seed(4323)
#Expit transform
expit<-function(x){return(1/(1+exp(-x)))}
n<-1000
p<-10
x.index<-1:p
#Two spurious
#Two instruments
#Three confounders
#Three predictors of outcome only 
al.select<-c(0,0,1,1,1,1,1,0,0,0)
be.select<-c(0,0,0,0,1,1,1,1,1,1)
#Confounders
confound.list<-x.index[al.select == 1 & be.select==1]
#Instruments
intrument.list<-x.index[al.select ==1 & be.select==0]
#Predictors of outcome only
outcomeonly.list<-x.index[al.select ==0  & be.select==1]
#Spurious variables
spurious.list<-x.index[al.select ==0  & be.select==0]
#Predictors of outcome
outcome.list<-x.index[be.select ==1]
#Predictors of exposure
treatment.list<-x.index[al.select ==1]
################################################
#Generate alpha parameters (for treatment allocation)
al<-rnorm(p+1,0,1.0);al[1]<-al[1]+1
#Generate beta parameters (for outcome regression)
be<-rnorm(p+1,0,3.0)
theta<-2.0
sig.Y<-1
p0<-2
p1<-2
p2<-3
p3<-3
Sig0<-rWishart(1,p0+1,diag(p0)/(p0+1))[,,1]
Sig1<-rWishart(1,p1+1,diag(p1)/(p1+1))[,,1]
Sig2<-rWishart(1,p2+1,diag(p2)/(p2+1))[,,1]
Sig3<-rWishart(1,p3+1,diag(p3)/(p3+1))[,,1]
Sig<-matrix(0,p,p)
Sig[1:p0,1:p0]<-Sig0
Sig[(p0+1):(p0+p1),(p0+1):(p0+p1)]<-Sig1
Sig[(p0+p1+1):(p0+p1+p2),(p0+p1+1):(p0+p1+p2)]<-Sig2
Sig[(p0+p1+p2+1):(p0+p1+p2+p3),(p0+p1+p2+1):(p0+p1+p2+p3)]<-Sig3
#Predictors
X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)
#Data simulation
Xmat<-cbind(rep(1,n),X)
eta.Z<-Xmat %*% (c(1,al.select)*al)
p.Z<-expit(eta.Z)
Z<-rbinom(n,1,p.Z)
eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)
hist(p.Z)
library(MASS)
set.seed(4323)
#Expit transform
expit<-function(x){return(1/(1+exp(-x)))}
n<-1000
p<-10
x.index<-1:p
#Two spurious
#Two instruments
#Three confounders
#Three predictors of outcome only 
al.select<-c(0,0,1,1,1,1,1,0,0,0)
be.select<-c(0,0,0,0,1,1,1,1,1,1)
#Confounders
confound.list<-x.index[al.select == 1 & be.select==1]
#Instruments
intrument.list<-x.index[al.select ==1 & be.select==0]
#Predictors of outcome only
outcomeonly.list<-x.index[al.select ==0  & be.select==1]
#Spurious variables
spurious.list<-x.index[al.select ==0  & be.select==0]
#Predictors of outcome
outcome.list<-x.index[be.select ==1]
#Predictors of exposure
treatment.list<-x.index[al.select ==1]
################################################
#Generate alpha parameters (for treatment allocation)
al<-rnorm(p+1,0,1.0);al[1]<-al[1]+1
#Generate beta parameters (for outcome regression)
be<-rnorm(p+1,0,4.0)
theta<-2.0
sig.Y<-1
p0<-2
p1<-2
p2<-3
p3<-3
Sig0<-rWishart(1,p0+1,diag(p0)/(p0+1))[,,1]
Sig1<-rWishart(1,p1+1,diag(p1)/(p1+1))[,,1]
Sig2<-rWishart(1,p2+1,diag(p2)/(p2+1))[,,1]
Sig3<-rWishart(1,p3+1,diag(p3)/(p3+1))[,,1]
Sig<-matrix(0,p,p)
Sig[1:p0,1:p0]<-Sig0
Sig[(p0+1):(p0+p1),(p0+1):(p0+p1)]<-Sig1
Sig[(p0+p1+1):(p0+p1+p2),(p0+p1+1):(p0+p1+p2)]<-Sig2
Sig[(p0+p1+p2+1):(p0+p1+p2+p3),(p0+p1+p2+1):(p0+p1+p2+p3)]<-Sig3
#Predictors
X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)
#Data simulation
Xmat<-cbind(rep(1,n),X)
eta.Z<-Xmat %*% (c(1,al.select)*al)
p.Z<-expit(eta.Z)
Z<-rbinom(n,1,p.Z)
eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)
source('GenerateSetUp2.R')
filecode<-'Set2'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set2.RData')
t2<-proc.time()[3]
source('GenerateSetUp3.R')
filecode<-'Set3'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set3.RData')
t3<-proc.time()[3]
source('GenerateSetUp4.R')
filecode<-'Set4'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set4.RData')
t4<-proc.time()[3]
q()
t0<-proc.time()[3]
source('GenerateSetUp.R')
filecode<-'Set1'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set1.RData')
t1<-proc.time()[3]
print(t1-t0)
source('GenerateSetUp2.R')
filecode<-'Set2'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set2.RData')
t2<-proc.time()[3]
print(t2-t1)
source('GenerateSetUp3.R')
filecode<-'Set3'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set3.RData')
t3<-proc.time()[3]
print(t3-t2)
source('GenerateSetUp4.R')
filecode<-'Set4'
source("BasicSim.R") #Stratification
source("BasicSim-Reg.R") #PS Regression
source("BasicSim-IPW.R") #IPW
source("BasicSim-DR.R")  #AIPW
source("BasicSim-Match.R") #Matching
save(ests.mat,ests.mat.psr,ests.mat.ipw,ests.mat.dr,ests.mat.match,file='Set4.RData')
t4<-proc.time()[3]
print(t4-t3)
q()
