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")