################
### Preparation
################
rm(list = ls())
require(foreign)
require(Hmisc)
require(combinat)
require(simex)
EULIS_EUP <- read.dta("EULIS_EUP.dta", convert.factors=FALSE)
nman <- nrow(EULIS_EUP)
nrepl <- 2000
actorBSn <- actorBSnRand <- array(as.matrix(EULIS_EUP[,17:56]),
c(nman, 40, nrepl+1),
dimnames=list(1:nman, names(EULIS_EUP[,17:56]), 0:nrepl))
###############################
## Bootstrap (multinomial draw)
###############################
require(combinat)
n <- apply(actorBSn[1:nrow(actorBSn),,1],1,sum)
p <- actorBSn[,,1]/n
cat("Resample number: ")
for (i in 1:nrepl) {
cat(i," ")
actorBSn[,,i] <- rmultinomial(n, p)
}
cat("\n")
############
## Positions
############
pnsBS <-
((log(actorBSn[,"p101s_107",]+0.5)-log(actorBSn[,"p101n_109",]+0.5)) +
(log(actorBSn[,"p102s_108",]+0.5)-log(actorBSn[,"p102n_110",]+0.5)) +
(log(actorBSn[,"p103s_602",]+0.5)-log(actorBSn[,"p103n_601",]+0.5))) /3
pnsMean_rep <- apply(pnsBS,1,mean)
pnsSE_rep <- apply(pnsBS,1,sd)
plrBS <-
((log(actorBSn[,"p201r_104",]+0.5)-log(actorBSn[,"p201l_105",]+actorBSn[,"p201l_106",]+0.5)) +
(log(actorBSn[,"p202r_605",]+0.5)-log(actorBSn[,"p202l_201",]+actorBSn[,"p202l_202",]+0.5)) +
(log(actorBSn[,"p203r_303",]+0.5)-log(actorBSn[,"p203l_404",]+actorBSn[,"p203l_405",]+0.5)) +
(log(actorBSn[,"p204r_401",]+0.5)-log(actorBSn[,"p204l_412",]+actorBSn[,"p204l_413",]+0.5)) +
(log(actorBSn[,"p205r_402",]+0.5)-log(actorBSn[,"p205l_403",]+0.5)) +
(log(actorBSn[,"p206r_407",]+0.5)-log(actorBSn[,"p206l_406",]+0.5)) +
(log(actorBSn[,"p207r_414",]+0.5)-log(actorBSn[,"p207l_409",]+0.5)) +
(log(actorBSn[,"p208r_410",]+0.5)-log(actorBSn[,"p208l_416",]+actorBSn[,"p208l_501",]+0.5)) +
(log(actorBSn[,"p209r_505",]+0.5)-log(actorBSn[,"p209l_503",]+actorBSn[,"p209l_504",]+0.5)) +
(log(actorBSn[,"p210r_603",]+0.5)-log(actorBSn[,"p210l_604",]+0.5)) +
(log(actorBSn[,"p211r_608",]+0.5)-log(actorBSn[,"p211l_607",]+0.5)) +
(log(actorBSn[,"p212r_702",]+0.5)-log(actorBSn[,"p212l_701",]+0.5)) +
(log(actorBSn[,"p213r_704",]+0.5)-log(actorBSn[,"p213l_705",]+0.5))) /13
plrMean_rep <- apply(plrBS,1,mean)
plrSE_rep <- apply(plrBS,1,sd)
eulisnsBS <-
(((log(actorBSn[,"p101s_107",]+0.5)-log(actorBSn[,"p101n_109",]+0.5))*EULIS_EUP$s101) +
((log(actorBSn[,"p102s_108",]+0.5)-log(actorBSn[,"p102n_110",]+0.5))*EULIS_EUP$s102) +
((log(actorBSn[,"p103s_602",]+0.5)-log(actorBSn[,"p103n_601",]+0.5))*EULIS_EUP$s103)) /100
eulisnsMean_rep <- apply(eulisnsBS,1,mean)
eulisnsSE_rep <- apply(eulisnsBS,1,sd)
eulislrBS <-
(((log(actorBSn[,"p201r_104",]+0.5)-log(actorBSn[,"p201l_105",]+actorBSn[,"p201l_106",]+0.5))*EULIS_EUP$s201) +
((log(actorBSn[,"p202r_605",]+0.5)-log(actorBSn[,"p202l_201",]+actorBSn[,"p202l_202",]+0.5))*EULIS_EUP$s202) +
((log(actorBSn[,"p203r_303",]+0.5)-log(actorBSn[,"p203l_404",]+actorBSn[,"p203l_405",]+0.5))*EULIS_EUP$s203) +
((log(actorBSn[,"p204r_401",]+0.5)-log(actorBSn[,"p204l_412",]+actorBSn[,"p204l_413",]+0.5))*EULIS_EUP$s204) +
((log(actorBSn[,"p205r_402",]+0.5)-log(actorBSn[,"p205l_403",]+0.5))*EULIS_EUP$s205) +
((log(actorBSn[,"p206r_407",]+0.5)-log(actorBSn[,"p206l_406",]+0.5))*EULIS_EUP$s206) +
((log(actorBSn[,"p207r_414",]+0.5)-log(actorBSn[,"p207l_409",]+0.5))*EULIS_EUP$s207) +
((log(actorBSn[,"p208r_410",]+0.5)-log(actorBSn[,"p208l_416",]+actorBSn[,"p208l_501",]+0.5))*EULIS_EUP$s208) +
((log(actorBSn[,"p209r_505",]+0.5)-log(actorBSn[,"p209l_503",]+actorBSn[,"p209l_504",]+0.5))*EULIS_EUP$s209) +
((log(actorBSn[,"p210r_603",]+0.5)-log(actorBSn[,"p210l_604",]+0.5))*EULIS_EUP$s210) +
((log(actorBSn[,"p211r_608",]+0.5)-log(actorBSn[,"p211l_607",]+0.5))*EULIS_EUP$s211) +
((log(actorBSn[,"p212r_702",]+0.5)-log(actorBSn[,"p212l_701",]+0.5))*EULIS_EUP$s212) +
((log(actorBSn[,"p213r_704",]+0.5)-log(actorBSn[,"p213l_705",]+0.5))*EULIS_EUP$s213)) /100
eulislrMean_rep <- apply(eulislrBS,1,mean)
eulislrSE_rep <- apply(eulislrBS,1,sd)
##########################
## Write out Stata dataset
##########################
save.image(file="EULIS_EUP_rep.Rdata", safe=FALSE)
EULIS_EUPbs <- data.frame(cbind(EULIS_EUP[,1:134],
pnsMean_rep,
pnsSE_rep,
plrMean_rep,
plrSE_rep,
eulisnsMean_rep,
eulisnsSE_rep,
eulislrMean_rep,
eulislrSE_rep))
write.dta(EULIS_EUPbs, "EULIS_EUP_rep.dta", version=7)