Código em construção, ainda não é a versão final!
# Dados para teste
x <- c(0,0,0,1,0,0,NA,NA,1,0,0,1,NA,1,NA,NA,0,NA,NA,0,0,0,0,NA,0,0,0,1,0,NA,1,NA,0,1,0,1,NA,0,NA,0,0,0,1,0,0,0,NA,1,0,NA,0,0,NA,1,0,1,1,0,0,0,1,1,1,0,NA,NA,NA,1,NA,NA,NA,NA,NA,NA,0,0,0,0,NA,0,1,0,0,1,0,0,0,1,1,1,1,1,NA,NA,0,0,0,NA,0,NA,0,1,0,1,0,0,1,0,0,NA,NA,0,1,NA,NA,0,0,1,0,1,1,1,1,1,NA,0,1,NA,0,1,1,0,1,0)
y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0)
prev <- function(x,y, datatype=c("rrt", "qd", "RRT.qd"), probs=c(0.25, 0.5833), nsim=1000) # função prev com argumentos: x, y, probs (probabilidades de "sim" forçado e "diga a verdade", ver Hox & Lensvelt-Mulders 2005), teste (se usuário quer realizar teste de permutação), nsim (número de simulações)
{
if(datatype=="rrt")
{
rrt <- function(x) ((sum(na.omit(x))/length(na.omit(x)))-probs[1])/probs[2] # fórmula de Hox & Lensvelt-Mulders (2005) para estimar prevalência de comportamento com base em dados RRT
a <- rrt(x) # criando objeto a com estimativa de prevalência para os dados de respostas randomizadas
rrt.boot <- rep(NA, nsim) # criando objeto rrt.boot para receber valores de simulação
for(i in 1:nsim) # criando contador
{
sample.x <- sample(x, replace=TRUE) # reamostrando vetor de dados de respostas randomizadas com substituição
rrt.boot[i] <- rrt(sample.x) # preenchendo objeto rrt.boot
}
b <- quantile(rrt.boot, probs=c(0.025, 0.975), na.rm=T) # criando objeto b com quantis referentes ao intervalo de confiança de 95%
n.NAx <- (length(x))-(length(na.omit(x))) # calcula numero de NAs no vetor de dados x
cat("\n\t", n.NAx," valores NA omitidos de x\n") # insere texto informando numero de valores de NA omitidos de x
ret.rrt <- list("Prevalencia RRT" = a, "intervalo de confiança de 95% RRT" = b) # cria objeto com lista de itens a serem exibidos no return
return(ret.rrt)
}
if(datatype=="qd")
{
qd <- function(y) sum(na.omit(y))/length(na.omit(y)) # fórmula para estimar prevalência de comportamento com base em dados de questionamento direto
a2 <- qd(y) # criando objeto a2 com estimativa de prevalência para os dados de questionamento direto
qd.boot <- rep(NA, nsim) # criando objeto qd.boot para receber valores de simulação
for(i in 1:nsim) # criando contador
{
sample.y <- sample(y, replace=TRUE) # reamostrando vetor de dados de questionamento direto com substituição
qd.boot[i] <- qd(sample.y) # preenchendo objeto qd.boot
}
b2 <- quantile(qd.boot, probs=c(0.025, 0.975), na.rm=T) # criando objeto b2com quantis referentes ao intervalo de confiança de 95%
n.NAy <- (length(y))-(length(na.omit(y))) # calcula numero de NAs no vetor de dados y
cat("\t", n.NAy," valores NA omitidos de y\n\n")
ret.qd <- list("Prevalencia QD" = a2, "intervalo de confiança de 95% QD" = b2) # cria objeto com lista de itens a serem exibidos no return
return(ret.qd)
}
if(datatype=="rrt.qd")
{
rrt <- function(x) ((sum(na.omit(x))/length(na.omit(x)))-probs[1])/probs[2] # fórmula de Hox & Lensvelt-Mulders (2005) para estimar prevalência de comportamento com base em dados RRT
a <- rrt(x) # criando objeto a com estimativa de prevalência para os dados de respostas randomizadas
rrt.boot <- rep(NA, nsim) # criando objeto rrt.boot para receber valores de simulação
for(i in 1:nsim) # criando contador
{
sample.x <- sample(x, replace=TRUE) # reamostrando vetor de dados de respostas randomizadas com substituição
rrt.boot[i] <- rrt(sample.x) # preenchendo objeto rrt.boot
}
b <- quantile(rrt.boot, probs=c(0.025, 0.975), na.rm=T) # criando objeto b com quantis referentes ao intervalo de confiança de 95%
qd <- function(y) sum(na.omit(y))/length(na.omit(y)) # fórmula para estimar prevalência de comportamento com base em dados de questionamento direto
a2 <- qd(y) # criando objeto a2 com estimativa de prevalência para os dados de questionamento direto
qd.boot <- rep(NA, nsim) # criando objeto qd.boot para receber valores de simulação
for(i in 1:nsim) # criando contador
{
sample.y <- sample(y, replace=TRUE) # reamostrando vetor de dados de questionamento direto com substituição
qd.boot[i] <- qd(sample.y) # preenchendo objeto qd.boot
}
b2 <- quantile(qd.boot, probs=c(0.025, 0.975), na.rm=T) # criando objeto b2com quantis referentes ao intervalo de confiança de 95%
boxplot(rrt.boot, qd.boot, xlab="Metodo de amostragem", ylab="Prevalência", xaxt="n") # boxplot com estimativas de prevalência obtidas pelos dois métodos
axis(side=1, at=(1:2), labels=c("RRT", "QD")) # adicionando labels ao boxplot
dif <- rrt.boot-qd.boot # calculando a diferença entre as estimativas de prevalência obtidas pelos dois métodos
n <- sum(abs(dif) >= a-a2) # calculando quantas vezes o módulo do valor observado foi obtido ao acaso
prob <- n/length(dif) # calculando a probabilidade de o valor observado se dever ao acaso
x11() # abrindo novo graphic device
hist(dif, ylab="Frequência", xlab="RRT-QD", main=NULL) # plota histograma com diferenças simuladas entre as estimativas de prevalência obtidas pelos dois métodos
abline(v=a-a2, col="red") # indica a posição do valor observado no histograma
n.NAx <- (length(x))-(length(na.omit(x))) # calcula numero de NAs no vetor de dados x
cat("\n\t", n.NAx," valores NA omitidos de x\n") # insere texto informando numero de valores de NA omitidos de x
n.NAy <- (length(y))-(length(na.omit(y))) # calcula numero de NAs no vetor de dados y
cat("\t", n.NAy," valores NA omitidos de y\n\n")
ret.rrt.qd <- list("Prevalencia RRT" = a, "intervalo de confiança de 95% RRT" = b, "Prevalencia QD" = a2, "intervalo de confiança de 95% QD" = b2, "Probabilidade de diferença entre estimativas RRT e QD ter sido obtida ao acaso"=prob) # cria objeto com lista de itens a serem exibidos no return
return(ret.rrt.qd)
}
}
# testando:
prev(x, datatype="rrt")
prev(,y, datatype="qd")
prev(x, y, datatype="rrt.qd")