I'm trying to run a R code made by another researcher, but it always gives me error when it comes to return risk ratios from posterior simulations based on a zelig function modelling. The replication's material can be downloaded here. More specifically, I have trouble running the "_replic.ind2006.AJPS" code, which makes use of "individualdataset2006" data (it's all inside the zip archive).
The following code is written from line 74 to 92.
a matched.data <- match.data(matching)
matched.data$lula <- as.numeric(as.character(matched.data$lula))
matched.data$renda <- factor(matched.data$renda)#refactor because highest category was dropped
matched.data$renda2 <- factor(matched.data$renda2)#re
fit.ate <- lm(lula ~ bf + sexo + escola + renda + idade + region, data=matched.data, weights=matched.data$weights)
fit1 <- zelig(lula ~ bf + sexo + escola + renda + idade + region, model="logit",data=tmp)
fit1m <- zelig(lula ~ bf + sexo + escola + renda + idade + region + metropolitan, model="logit",weights=matched.data$weights,data=matched.data)
#FOR ALL OBSERVATIONS
inc <- levels(new.data$renda)
rr <- pred.bf0 <- pred.bf1 <- sig <-vector(length=length(levels(new.data$renda)))
for(i in 1:length(inc)){
bf0 <- setx(fit1 , bf=0, renda=inc[i])
bf1 <- setx(fit1 ,bf=1, renda=inc[i])
rr[i] <-summary(sim(fit1,x=bf0,x1=bf1))$qi.stats$rr[1]
sig[i] <- summary(sim(fit1,x=bf0,x1=bf1))$qi.stats$rr[3]>1
pred.bf0[i] <- mean(sim(fit1,x=bf0)$qi$ev)
pred.bf1[i] <- mean(sim(fit1,x=bf1)$qi$ev)
}
This part starts by using matching procedures, and it's followed by estimating models through logit function. So where does the problem begin? It begins precisely when running the loop
for(i in 1:length(inc))
where length(inc) represents the number of a categorical variable's levels, which means, in practice, 1:4.
The code fails to run
rr[i] <-summary(sim(fit1,x=bf0,x1=bf1))$qi.stats$rr
which should store risk ratios from the estimation. It gives the following error:
Error in rr[i] <- summary(sim(fit1, x = bf0, x1 = bf1))$qi.stats$rr[1] : replacement has length zero
I have first thought it was a looping problem that could be resolved if I stored "manually" the risk ratios, but running
summary(sim(fit1, x = bf0, x1 = bf1))$qi.stats$rr[1]
alone returns NULL. Indeed, if I just return the summary without $qi.stats$rr[1], there is no risk ratios in the output, it gives me only expected values, predicted values and first differences.
The problem of "replacement has zero length" was also reported when running
pred.bf0[i] <- mean(sim(fit1,x=bf0)$qi$ev)
but the difference is that I can get expected values through get_qi, so it's not really a big problem right now, although it's bothersome.
Does anyone have any clues about it?
Thank you in advance.
