# PDItest.R
# 2008-Ene-05 Viper@upch.edu.pe
# 2008-Ene-15 ponderaciones
# 2008-Ene-25 estimaciones
# 2008-Feb-14 viviendas
# 2008-Mar-13 panel
# 2008-Mar-24 N5 y N3, glm descriptivos, ajuste por MaEdu
# 2008-Mar-31 Crecer
# 2008-Mar-05 gráficos para informe final
# 2008-Mar-09 panel y gráficos
# 2008-Mar-11 solo Crecer
# 2008-Abr-18 Cobertura VL, Beneficiarios individuales


# ======================================================================================================================== #
# call: source("T:/PDItest.R",echo=TRUE,print.eval=TRUE)

#par(ask = dev.interactive(orNone = TRUE))
sink("G:/Temp/DRAFT/MONIN/p.txt")
pdf("G:/Temp/DRAFT/MONIN/p.pdf")

date()

#library(foreign) 
#library(stats)
#library(boot)
library(survey)
library(MASS) 
library(epicalc) 

#setwd
WkOrDay<- function ( wkd ) { 
	y <- c("Lun","Mar","Mié","Jue","Vie","Sáb","Dom")
	w <- 0
	for (i in 1:7) { if (y[i]==wkd) w<-i }
	w }
FSumX<- function ( x, g, xlab ) {
	x<-by(x, g, sum, na.rm=TRUE)
	paste(	xlab, format(x[1]+x[2]+x[3]+x[4]+x[5],trim=TRUE,width=8),
		format(x[1],trim=TRUE,width=8),
		format(x[2],trim=TRUE,width=8),
		format(x[3],trim=TRUE,width=8),
		format(x[4],trim=TRUE,width=8),
		format(x[5],trim=TRUE,width=8) ) }
FValN<- function ( n, k ) {  }
FPct <- function ( x ) { paste(ifelse(x==1,"",ifelse(x<0.1,"  "," ")),format(round(100*x,2),nsmall=2),"%",sep="") }
FPctP<- function ( p, v ) { paste ( FPct(p), " [",FPct(sqrt(v)),"]", sep="" ) }
FPctD<- function ( x, d ) { paste ( " ",substr(d,1,1),FPct(x),substr(d,2,2), sep="" ) }
FPRowR<- function ( xnum, xden, xdes, xlab ) { 
	x=svyratio(xnum, xden, design=xdes, na.rm=TRUE);
	y=svyby(xnum,by=~Estrato,denominator=xden,design=xdes,svyratio,na.rm=TRUE);
	paste(xlab,
		"  ", FPctP(x$ratio[1],x$var[1]),
		"  ", FPctP(y$statistics.ratio[[1]],y$statistics.var[[1]]),
		"  ", FPctP(y$statistics.ratio[[2]],y$statistics.var[[2]]),
		"  ", FPctP(y$statistics.ratio[[3]],y$statistics.var[[3]]),
		"  ", FPctP(y$statistics.ratio[[4]],y$statistics.var[[4]]),
		"  ", FPctP(y$statistics.ratio[[5]],y$statistics.var[[5]]), "" ) }
FPRowM<- function ( xnum, xdes, xlab ) { 
	x=svymean(xnum,xdes,na.rm=TRUE,deff=TRUE);
	y=svyby(xnum,by=~Estrato,design=xdes,svymean,na.rm=TRUE,deff=TRUE);
	paste(xlab,
		"  ", FPctP(as.numeric(x[1]),as.numeric(attributes(x)[2])),
		"  ", FPctP(y[1,2],y[1,3]^2),
		"  ", FPctP(y[2,2],y[2,3]^2),
		"  ", FPctP(y[3,2],y[3,3]^2),
		"  ", FPctP(y[4,2],y[4,3]^2),
		"  ", FPctP(y[5,2],y[5,3]^2), "" ) }
FPRwR2<- function ( xnum, xden, xdes, xlab, xby=~Estrato, fplot=FALSE ) { 
	x=svyratio(xnum, xden, design=xdes, na.rm=TRUE);
	y=svyby(xnum,by=xby,denominator=xden,design=xdes,svyratio,na.rm=TRUE);
	s1="";for (i in 1:length(y$statistics.ratio)) s1=paste (s1,"  ", FPctD(y$statistics.ratio[[i]],"  "),sep="");
	s2="";for (i in 1:length(y$statistics.var)) s2=paste (s2,"  ", FPctD(sqrt(y$statistics.var[[i]]),"[]"),sep="");
	if ( fplot ) {
		a=c(x$var[1]); for (i in 1:length(y$statistics.var)) a=c(a,y$statistics.var[[i]])
		b=c(x$ratio[1]); for (i in 1:length(y$statistics.ratio)) b=c(b,y$statistics.ratio[[i]])
		plot(c(0,as.numeric(y[,1])), c(x$ratio[1],y$statistics.ratio),
			ylim=c(0,1), type = "n", bty="n", xaxt="n", las=1,
			xlab="INS/CENAN MONIN 2007-Nov-19 a 2008-Mar-03", ylab="prevalencia e IC95", main=xlab )
		grid(nx=NA, ny=NULL, col="lightgray", lty="solid", lwd=1, equilogs = TRUE)
		arrows(	c(0,as.numeric(y[,1])), (b-1.96*(sqrt(a))), 
			c(0,as.numeric(y[,1])), (b+1.96*(sqrt(a))), 
			code=3, length=0.1, angle=90, col=c("red",rep("blue",times=length(y$statistics.ratio))) )
		points(c(0,as.numeric(y[,1])), c(x$ratio[1],y$statistics.ratio), col=c("red",rep("blue",times=length(y$statistics.ratio))), pch=19, cex=2)
		text(c(0,as.numeric(y[,1])),rep(0,times=1+length(y$statistics.ratio)), c("Total",levels(y[,1])), col=c("red",rep("blue",times=length(y$statistics.ratio)))) }
	cat(	paste(xlab,"  ", FPctD(x$ratio[1],"  "),s1, sep="" ),
		paste(format(" ",trim=TRUE,width=nchar(xlab)),"  ", FPctD(sqrt(x$var[1]),"[]"),s2, sep="" ),
		sep="\n") }
FPRwM2<- function ( xnum, xdes, xlab, xby=~Estrato, fplot=FALSE ) { 
	x=svymean(xnum,xdes,na.rm=TRUE,deff=FALSE);
	y=svyby(xnum,by=xby,design=xdes,svymean,na.rm=TRUE,deff=FALSE);
	s1="";for (i in 1:length(y[,2])) s1=paste (s1,"  ", FPctD(y[i,2],"  "),sep="");
	s2="";for (i in 1:length(y[,3])) s2=paste (s2,"  ", FPctD(y[i,3],"[]"),sep="");
	if ( fplot ) {
		plot(c(0,as.numeric(y[,1])), c(x[1],y[,2]),
			ylim=c(0,1), type = "n", bty="n", xaxt="n", las=1,
			xlab="INS/CENAN MONIN 2007-Nov-19 a 2008-Mar-03", ylab="prevalencia e IC95", main=xlab )
		grid(nx=NA, ny=NULL, col="lightgray", lty="solid", lwd=1, equilogs = TRUE)
		arrows(	c(0,as.numeric(y[,1])), (c(x[1],y[,2])-1.96*c(sqrt(as.numeric(attributes(x)[2])),y[,3])), 
			c(0,as.numeric(y[,1])), (c(x[1],y[,2])+1.96*c(sqrt(as.numeric(attributes(x)[2])),y[,3])), 
			code=3, length=0.1, angle=90,
			col=c("red",rep("blue",times=length(y[,2]))))
		points(c(0,as.numeric(y[,1])), c(x[1],y[,2]), col=c("red",rep("blue",times=length(y[,2]))), pch=19, cex=2)
		text(c(0,as.numeric(y[,1])),rep(0,times=1+length(y[,2])), c("Total",levels(y[,1])), col=c("red",rep("blue",times=length(y[,2])))) }
	cat(	paste(xlab,"  ", FPctD(as.numeric(x[1]),"  "),s1, sep="" ),
		paste(format(" ",trim=TRUE,width=nchar(xlab)),"  ", FPctD(sqrt(as.numeric(attributes(x)[2])),"[]"),s2, sep="" ),
		sep="\n" ) }
FPTbMR<- function ( xvec, xdes, xlab, xden=NA, xby=~Estrato, fratio=FALSE, fplot=FALSE ) { 
	nv = length(xvec)
	nc = max(nchar(xlab))
	t = svytable(xby,xdes[[1]],round=TRUE)
	nf = length(t)	#length(levels(xby))
	xp=rep(NA,times=nv*(nf+2));dim(xp)=c(nf+2,nv)
	se=rep(NA,times=nv*(nf+2));dim(se)=c(nf+2,nv)
	for (i in 1:nv) {
		if ( fratio ) {
			x=svyratio(xvec[i][[1]], xden[i][[1]], design=xdes[[i]], na.rm=TRUE);
			y=svyby(xvec[i][[1]],by=xby,denominator=xden[i][[1]],design=xdes[[i]],svyratio,na.rm=TRUE);
			xp[1,i]=x$ratio[1]
			se[1,i]=sqrt(x$var[1])
			for (j in 1:nf) {
				xp[1+j,i]=y$statistics.ratio[[j]]
				se[1+j,i]=sqrt(y$statistics.var[[j]]) }
			xp[nf+2,i]=NA }
		else {
			x=svymean(xvec[i][[1]],design=xdes[[i]],na.rm=TRUE,deff=FALSE)
			y=svyby(xvec[i][[1]],by=xby,design=xdes[[i]],svymean,na.rm=TRUE,deff=FALSE);
			xp[1,i]=as.numeric(x[1])
			se[1,i]=sqrt(as.numeric(attributes(x)[2]))
			for (j in 1:nf) {
				xp[1+j,i]=y[j,2]
				se[1+j,i]=y[j,3] }
			if ( TRUE ) {
				# length(levels(xvec[i][[1]]))>1
				z = svychisq(as.formula(paste("~",all.vars(xby),"+",all.vars(xvec[i][[1]]))),xdes[[1]],round=TRUE)
				xp[nf+2,i]=round(z$p.value,6) }
			else { xp[nf+2,i]=NA }  }
		}
	if ( fplot ) {
		# proporción versus indicador
		plot(1:nv, xp[1,],
			ylim=c(0,1), type = "n", bty="n", xaxt="n", las=1,
			xlab="INS/CENAN MONIN 2007-Nov-19 a 2008-Mar-03", ylab="proporción e IC95", main="Nacional" )
		grid(nx=NA, ny=NULL, col="lightgray", lty="solid", lwd=1, equilogs = TRUE)
		arrows(1:nv, xp[1,]-1.96*se[1,], 1:nv, xp[1,]+1.96*se[1,], code=3, length=0.1, angle=90, col=1:nv)
		points(1:nv, xp[1,], col=1:nv, pch=19, cex=2)
		text(1:nv,rep(0,times=nv), xlab, col=1:nv,cex=0.75)
		# proporción versus indicador y grupo
		plot(1:nv, xp[1,],
			ylim=c(0,1), type = "n", bty="n", xaxt="n", las=1,
			xlab="INS/CENAN MONIN 2007-Nov-19 a 2008-Mar-03", ylab="proporción e IC95", main="Nacional y por Grupos" )
		grid(nx=NA, ny=NULL, col="lightgray", lty="solid", lwd=1, equilogs = TRUE)
		for (j in 1:(nf+1)) { 
			lines(1:nv, xp[j,], col=j, lty=1, pch=j)
			points(1:nv, xp[j,], col=j, pch=j ) }
		text(1:nv,rep(0,times=nv), xlab, col=1:nv,cex=0.75)
		legend("topright", bg="white", legend=c("Total",labels(t)[[1]]),
			col=1:(nf+1), pch=1:(nf+1), lty=1:(nf+1), cex=0.75)
		# proporción versus grupo, para cada indicador
		for (i in 1:nv ) {
			plot(0:nf, c(xp[1:(nf+1),i]),
				ylim=c(0,1), type = "n", bty="n", xaxt="n", las=1,
				xlab="INS/CENAN MONIN 2007-Nov-19 a 2008-Mar-03", ylab="proporción e IC95", main=paste("Indicador ",xlab[[i]]) )
			grid(nx=NA, ny=NULL, col="lightgray", lty="solid", lwd=1, equilogs = TRUE)
			arrows(0:nf, (xp[1:(nf+1),i]-1.96*se[1:(nf+1),i]),0:nf, (xp[1:(nf+1),i]+1.96*se[1:(nf+1),i]), 
				code=3, length=0.1, angle=90,col=c("red",rep("blue",times=nf)))
			points(0:nf, xp[1:(nf+1),i], col=c("red",rep("blue",times=nf)), pch=19, cex=2)
			text(0:nf,rep(0,times=1+nf), c("Total",labels(t)[[1]]), col=c("red",rep("blue",times=nf+1)),cex=0.75) }
		}
	# tabla indicador versus grupo ( p y ES )
	s1="Total";for (j in 1:nf) s1=paste (s1,"  ", labels(t)[[1]][j],"  ",sep="");
	cat(	paste(format(" ",trim=TRUE,width=nc),"  ", s1, sep="" ),
		format(" ",trim=TRUE,width=nc),
		sep="\n" )
	s1=paste("n ponderado ",format(sum(t),trim=TRUE,width=nc),sep="")
	for (j in 1:nf) s1=paste (s1,"  ", format(t[j],trim=TRUE,width=10),"  ",sep="")
	cat(	paste(format(" ",trim=TRUE,width=nc),"  ", s1, sep="" ),
		format(" ",trim=TRUE,width=nc),
		sep="\n" )
	for (i in 1:nv ) {
		s1="";for (j in 1:(nf+1)) s1=paste (s1,"  ", FPctD(xp[j,i],"  "),sep="");s1=paste (s1,"  ", xp[nf+2,i], sep="")
		s2="";for (j in 1:(nf+1)) s2=paste (s2,"  ", FPctD(se[j,i],"[]"),sep="");
		cat(	paste(format(xlab[[i]],trim=TRUE,width=nc),"  ", s1, sep="" ),
			paste(format(" ",trim=TRUE,width=nc),"  ", s2, sep="" ),
			sep="\n" ) }
	}
	#FPRwR2 ( xvec[i][[1]], xden[i][[1]], xdes[[i]], format(xlab[i],trim=TRUE,width=nc+2), xby=xby, fplot=fplot )
	#FPRwM2 ( xvec[i][[1]], xdes[[i]], format(xlab[i],trim=TRUE,width=nc+2), xby=xby, fplot=fplot )
FPack1<- function ( xby=~Estrato, fplot=FALSE ) { 
	t<-"===== Niños 0-59 meses, Resumen";t
	FPTbMR ( xvec=c(~N5Nut,~N5Ali,~N5SaY,~N5ATE,~N5CPR,~N5VIV,~N5COM),
	 xdes=list(p5,p5,p5,p5,p5,p5,p5),
	 xlab=c("Nutr","Alim","Salud","Ate","Prep","Viv","Com"), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Nutrición";t
	FPTbMR ( xvec=c(~N5DTE,~N5ADE,~N5OBE,~N5Ane),
	 xdes=list(p5,p5,p5,p5),
	 xlab=c("TE<-2z","PT<-2z","Obesidad","Anemia"), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Alimentación";t
	FPTbMR ( xvec=c(~IPConLM,~N5LME,~IPVComio,~MasGra,~IPVitA,~IPVFe,~IPSalRx,~N5Ali),
	 xdes=list(subset(p3,EdadM<24),subset(p5,EdadM<6),subset(p3,EdadM<24),subset(p3,EdadM<24),subset(p3,EdadM<24),subset(p3,EdadM<24),p5,subset(p5,N5SaY==0)),
	 xlab=c("LM[<2]","LME[<6m]","3+SS[<2]","Grasa[<2]","Vit.A[<2]","S.Fe[<2]","Sal I","Alim.Enf."), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Salud";t
	FPTbMR ( xvec=c(~IPPFie2w,~IPPTos2w,~NER7,~NumeIR,~NumeRR,~IPPDia2w,~NED7,~NumeID,~NumeRD),
	 xdes=list(p5,p5,p5,p5,p5,p5,p5,p5,p5),
	 xlab=c("Fieb 2s","IRA 2s","pIRA","iIRA","rIRA","EDA 2s","pEDA","iEDA","rEDA"), 
	 xden=c(~one,~one,~DenomP,~DenomIR,~DenomRR,~one,~DenomP,~DenomID,~DenomRD), 
	 xby=xby, fratio=TRUE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Atención";t
	FPTbMR ( xvec=c(~N5CC0,~N5PIP,~IPNVBCG,~IPNVPOL,~IPNVDPT,~IPNVHEP,~IPNVHIB,~IPNVSAR,~IPNVPaR,~IPNVAMA),
	 xdes=list(p5,p5,p5,p5,p5,p5,p5,p5,p5,p5),
	 xlab=c("CRED<1a","PAI","BCG","Polio","DPT","Hep.","HiB","Sar.","Rub.Pap.","F.Am."), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Cuidado";t
	FPTbMR ( xvec=c(~MaPte,~IPMaEdu1,~IPMaEdu2,~IPNMSD,~IPNMCCD),
	 xdes=list(p5,p5,p5,p5,p5),
	 xlab=c("Presente","Alfabeta","Primaria","SD<6m","CRED c/Cont."), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Vivienda";t
	FPTbMR ( xvec=c(~IPAlumbra,~IPAgua,~IPAgRTC,~IPAgRTH,~IPExcre1,~IPExcre3,~IPTCoc1,~IPTCoc2),
	 xdes=list(p5,p5,p5,p5,p5,p5,p5,p5),
	 xlab=c("Elec.","Agua RD","Agua Cl","Agua H","Excr RD","Letr Ad.","Coc E/G","Coc As."), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	t<-"===== Niños 0-59 meses, Comunidad";t
	FPTbMR ( xvec=c(~IPCUI,~IPNBI,~IPNJtos,~IPSIS,~IPNBAA,~IPNPapB,~IPNBAA),
	 xdes=list(p5,p5,subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI==0)),
	 xlab=c("CUI","P=NBI>0","Juntos [P]","SIS [P]","AA [P]","Pap [P]","AA[noP]"), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	FPTbMR ( xvec=c(~IPNBPIN,~IPNBWW,~IPNBCOI,~IPNBFAR,~IPNBTBC,~IPNBCOP,~IPNBPRV,~IPNBVAL),
	 xdes=list(subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0)),
	 xlab=c("PIN[P]","WW[P]","CoI[P]","FAR[P]","TBC[P]","CoP[P]","Priv[P]","VL[P]"), 
	 xby=xby, fratio=FALSE, fplot=fplot ) 
	}
DescMod1<- function ( xdep, xdat, xset ) { 
	summary(stepAIC(glm(xdep~
	EdadM+Edad2+csexo+NBI+nHV109+cagua+nINEIE8+lINEIE8+ALMPM+FALT+FALT2+fINEIER+Estrato+WkDay,
	data=xdat,subset=xset,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE)) }
DescMod2<- function ( xdep, xdat, xset ) {
	summary(stepAIC(glm(xdep~
	EdadM+Edad2+csexo+NBI+cagua+nINEIE8+lINEIE8+ALMPM+FALT+FALT2+fINEIER+Estrato+WkDay,
	data=xdat,subset=xset,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE)) }
DescMod3<- function ( xdep, xdat, xset ) {
	summary(stepAIC(glm(xdep~
	EdadM+Edad2+csexo+nHV109+cagua+nINEIE8+lINEIE8+ALMPM+FALT+FALT2+fINEIER+Estrato+WkDay,
	data=xdat,subset=xset,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE)) }
DescMod4<- function ( xdep, xdat, xset ) {
	summary(stepAIC(glm(xdep~
	EdadM+Edad2+csexo+NBI+nHV109+nINEIE8+lINEIE8+ALMPM+FALT+FALT2+fINEIER+Estrato+WkDay,
	data=xdat,subset=xset,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE)) }
DescMod5<- function ( xdep, xdat, xset ) {
	summary(stepAIC(glm(xdep~
	EdadM*Edad2+nINEIE8*lINEIE8+FALT*FALT2*nVertOc,
	data=xdat,subset=xset,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE)) }


# ======================================================================================================================== #
t<-"===== MONIN 2008 0-59m Ronda 1";t;date()
# Archivo de conglomerados
tc<-read.csv("G:/Temp/Draft/MONIN/T.csv") 
tc$one = 1
# Archivo de niños
n5<-read.csv("G:/Temp/Draft/MONIN/P.csv") 
#n5$SWt=ifelse(n5$Cong1==35,1/((n5$P1*(1E-9))*0.125),n5$SWt)	# parche mientras se digitaba
n5$one = 1
# Variables antiguas
n5$P1=n5$P1*(1E-9)
n5$P2=n5$P2*(1E-6)
#n5$SWt=n5$SWt*(8/6)
#n5$sdom2=ifelse(n5$Estrato==4,3,n5$Estrato)
#n5$N5DTE=ifelse(n5$ZHAZ<(-2),1,0)
#n5$N5OBE=ifelse(n5$ZWHZ>(2),1,0)
#n5$N5Ali=ifelse(is.na(n5$ConAL),NA,
#	ifelse((n5$EdadM<6 & n5$ConLM==1 & n5$ConAL!=1) | (n5$EdadM>=6 & n5$EdadM<11 & n5$ConLM==1 & n5$AliOK==1 & n5$SupOK==1) | (n5$EdadM>=12 & n5$EdadM<24 & n5$ConLM==1 & n5$AliOK==1) | (n5$EdadM>=24 & n5$AliOK==1),1,0))
n5$pzho=ifelse(n5$HFAZ<(-2),1,0)
n5$HbCORR=ifelse(n5$HbCORR==-1,NA,n5$HbCORR)
n5$EstClu=as.numeric(n5$Estrato)*100+5*n5$Salida
n5$Estrato=factor(n5$Estrato,labels=c("Lima M","R Costa","Sierra U","Sierra R","Selva"))
# Denominadores para agregación por vivienda
ifelse(n5$EdadM<60 & !is.na(n5$ZHAZ),1,0) -> n5$N5DTED
ifelse(n5$EdadM<60 & !is.na(n5$ConAL),1,0) -> n5$N5CCDD
ifelse(n5$EdadM<12 & !is.na(n5$ConAL),1,0) -> n5$N5CC0D
ifelse(n5$EdadM<60 & !is.na(n5$ConAL),1,0) -> n5$N5AliD
ifelse(n5$EdadM<60 & !is.na(n5$VivOK),1,0) -> n5$N5VIVD
ifelse(n5$EdadM<60 & !is.na(n5$ZHAZ),1,0) -> n5$N5NutD
ifelse(n5$EdadM<60 & !is.na(n5$HbCORR),1,0) -> n5$N5AneD
# Indicadores de Estado Nutricional
ifelse(n5$EdadM<60 & !is.na(n5$ZHAZ),
	ifelse(n5$EdadM<60 & n5$ZHAZ < -2,1,0),NA) -> n5$N5DTE
ifelse(n5$EdadM<60 & !is.na(n5$ZWHZ),
	ifelse(n5$EdadM<60 & n5$ZWHZ < -2,1,0),NA) -> n5$N5ADE
ifelse(n5$EdadM<60 & !is.na(n5$ZWHZ),1,0) -> n5$N5OBED
ifelse(n5$EdadM<60 & !is.na(n5$ZWHZ),
	ifelse(n5$EdadM<60 & n5$ZWHZ > 2,1,0),NA) -> n5$N5OBE
ifelse(n5$EdadM<60 & !is.na(n5$HbCORR),
	ifelse(n5$EdadM<60 & n5$HbCORR < 11,1,0),NA) -> n5$N5Ane
# no incluimos  & !is.na(n5$HbCORR) por ser un subconjunto
n5$N5Nut <- ifelse(n5$EdadM<60 & !is.na(n5$ZHAZ),
	ifelse(n5$ZHAZ>=-2 & n5$ZWHZ>=-2 & n5$ZWHZ<2,1,0),NA)
# Indicadores de Alimentación
ifelse(n5$EdadM<60 & !is.na(n5$ConAL) & (
	(n5$EdadM<6 & n5$ConLM==1 & n5$ConAL!=1) | 
	(n5$EdadM>=6 & n5$EdadM<11 & n5$ConLM==1 & n5$AliOK==1 & n5$SupOK==1 & n5$NVVTA>0) | 
	(n5$EdadM>=12 & n5$EdadM<18 & n5$ConLM==1 & n5$AliOK==1 & n5$NVVTA>1) | 
	(n5$EdadM>=18 & n5$EdadM<24 & n5$ConLM==1 & n5$AliOK==1 & n5$NVVTA>2) | 
	(n5$EdadM>=24 & n5$AliOK==1)),1,0) -> n5$N5Ali
ifelse(n5$EdadM<6 & !is.na(n5$ConAL),ifelse(n5$EdadM<6 & n5$ConLM==1 & n5$ConAL!=1,1,0),NA) -> n5$N5LME
# Indicadores de Estado de Salud
ifelse(n5$EdadM<60 & !is.na(n5$ConAL) & n5$PDia2w==2 & n5$PTos2w==2 & n5$PFie2w==2,1,0) -> n5$N5San	# quincena
ifelse(n5$EdadM<60 & !is.na(n5$ConAL), ifelse(n5$NEDY==0 & n5$NERY==0,1,0),NA) -> n5$N5SaY	# ayer
# Indicadores de Atención (Vacunación)
ifelse(n5$EdadM<60 & !is.na(n5$ConAL) & (
	(n5$EdadM<1 & n5$NVBCG>0 & n5$NVHEP>0) | 
	(n5$EdadM>=1 & n5$EdadM<2 & n5$NVBCG>0 & n5$NVHEP>1 & n5$NVDPT>0 & n5$NVPOL>0 & n5$NVHIB>0) | 
	(n5$EdadM>=2 & n5$EdadM<3 & n5$NVBCG>0 & n5$NVHEP>1 & n5$NVDPT>1 & n5$NVPOL>1 & n5$NVHIB>1) | 
	(n5$EdadM>=3 & n5$EdadM<12 & n5$NVBCG>0 & n5$NVHEP>2 & n5$NVDPT>2 & n5$NVPOL>2 & n5$NVHIB>2) | 
	(n5$EdadM>=12 & n5$NVBCG>0 & n5$NVHEP>2 & n5$NVDPT>2 & n5$NVPOL>2 & n5$NVHIB>2 & n5$NVSAR>0 & n5$NVPaR>0 & n5$NVAMA>0)
	),1,0) -> n5$N5PIP
# Indicadores de Atención (Crecimiento)
ifelse(!is.na(n5$EdadM) & !is.na(n5$ConAL) & (
	(n5$EdadM<1 & n5$NCCD0>(1+(n5$EdadM-0)/0.5)) | 
	(n5$EdadM>=1 & n5$EdadM<2 & n5$NCCD0>(3+(n5$EdadM-1)/1)) | 
	(n5$EdadM>=2 & n5$EdadM<6 & n5$NCCD0>(5+(n5$EdadM-2)/2)) | 
	(n5$EdadM>=6 & n5$EdadM<24 & n5$NCCD0>(7+(n5$EdadM-6)/3)) | 
	(n5$EdadM>=24 & n5$NCCD0>(12+(n5$EdadM-24)/6)) ),1,0) -> n5$N5CCD
ifelse(!is.na(n5$EdadM) & !is.na(n5$ConAL) & (
	(n5$EdadM<1 & n5$NCCD0>as.integer(1+(n5$EdadM-0)/0.5)) | 
	(n5$EdadM>=1 & n5$EdadM<2 & n5$NCCD0>3) | 
	(n5$EdadM>=2 & n5$EdadM<6 & n5$NCCD0>as.integer(4+(n5$EdadM-2)/2)) | 
	(n5$EdadM>=6 & n5$EdadM<7 & n5$NCCD0>6) | 
	(n5$EdadM>=7 & n5$EdadM<9 & n5$NCCD0>7) | 
	(n5$EdadM>=9 & n5$EdadM<12 & n5$NCCD0>8) | 
	(n5$EdadM>=12 & n5$NCCD0>9) ),1,0) -> n5$N5CC0
# Porcentaje de último periodo cubierto por CRED
n5$ECCD0=ifelse (n5$ECCD0<0,-1,n5$ECCD0/30.4375)	# edad a última visita, meses
n5$UCCD0=ifelse (n5$ECCD0<0,100,ifelse(n5$EdadM>n5$ECCD0,n5$EdadM-n5$ECCD0,0))	# distancia a última visita, meses
x=ifelse(n5$EdadM<1,0.5,ifelse(n5$EdadM<2,1,ifelse(n5$EdadM<6,2,ifelse(n5$EdadM<24,3,ifelse(n5$EdadM<60,6,1)))))	# "duración" de la cobertura, según edad
n5$PCCD0=ifelse(n5$EdadM>=12,NA,ifelse(n5$UCCD0>x,0,ifelse(n5$UCCD0>x,(n5$EdadM-n5$UCCD0)-(n5$EdadM-x),n5$UCCD0))/x)	# porcentaje del último periodo
n5$PCCD30=ifelse(n5$EdadM>=12,NA,ifelse((n5$UCCD0+x)>n5$EdadM,1,ifelse((n5$UCCD0+x)<(n5$EdadM-1),0,1-(n5$EdadM-(n5$UCCD0+x)))))	# porcentaje de los últimos 30 dias "cubierto"
n5$NCCD30=round(n5$PCCD30*30,0)
n5$Deno30=ifelse(is.na(n5$NCCD30),NA,30)
n5$N5ATE <- ifelse(n5$N5PIP==1 & n5$N5CC0==1,1,0)
# Indicadores de Cuidadores Preparados
ifelse(n5$EdadM<60 & !is.na(n5$ConAL) & !is.na(n5$NMSD) & !is.na(n5$NMCCD) &
	(n5$MaOK==1 & n5$NMSD==1 & n5$NMCCD>0),1,0) -> n5$N5CPR
# Indicadores de Vivienda
ifelse(n5$EdadM<60 & !is.na(n5$VivOK) & n5$VivOK==1,1,0) -> n5$N5VIV
# Indicadores de Comunidad
ifelse(n5$EdadM<60 & !is.na(n5$ConAL) & (
	n5$CUI==1 & n5$SIS==1 & (n5$NBI==0 | (n5$NBI>0 & n5$NBAA>0 & n5$NJtos>0 & n5$NPapB>0))),1,0) -> n5$N5COM
# Indicador de Protección
n5$N5PRO <- ifelse(n5$N5Nut==1 & n5$N5Ali==1 & n5$N5San==1 & n5$N5ATE==1 & n5$N5CPR==1 & n5$N5VIV==1 & n5$N5COM==1,1,0)
# Variables de estratificación
n5$INEIER<-as.numeric(n5$Estrato)
n5$INEIE4<-n5$INEIE8
recode(INEIE4,8,6,dataFrame=n5)
recode(INEIE4,7,6,dataFrame=n5)
recode(INEIE4,n5$INEIE4>=3 & n5$INEIE4<6,3,dataFrame=n5)
n5$INEIE2<-n5$INEIE4
recode(INEIE2,n5$INEIE4>=1 & n5$INEIE4<6,1,dataFrame=n5)
recode(INEIER,4,3,dataFrame=n5)
n5$Crecer=factor(ifelse(n5$EtapaC>0,1,2),labels=c("Crecer","Resto"))
# variables de fecha
n5$fecha=as.date("1/1/1960", order="dmy")+ifelse(nchar(as.character(n5$FEnc))==0,NA,as.date(substr(as.character(n5$FEnc),1,10), order="dmy"))
n5$WkDay=0
for (i in 1:length(n5$WkDay)) n5$WkDay[i]<-WkOrDay(weekdays(as.Date(n5$fecha[i]),abbreviate=TRUE))
# variables indicador adicionales
n5$IPConLM=ifelse(n5$ConLM==2,0,n5$ConLM)
n5$IPVComio=ifelse(n5$VComio>=3,1,0)
n5$IPVitA=ifelse(n5$VitA==1,1,0)
n5$IPVFe=ifelse(n5$VFe>5,1,0)
n5$IPSalRx=ifelse(n5$SalRx==1,1,0)
n5$IPPFie2w=ifelse(n5$PFie2w==1,1,0)
n5$IPPTos2w=ifelse(n5$PTos2w==1,1,0)
n5$IPPDia2w=ifelse(n5$PDia2w==1,1,0)
n5$IPNVBCG=ifelse(n5$NVBCG>0,1,0)
n5$IPNVPOL=ifelse(n5$NVPOL>2,1,0)
n5$IPNVDPT=ifelse(n5$NVDPT>2,1,0)
n5$IPNVHEP=ifelse(n5$NVHEP>2,1,0)
n5$IPNVHIB=ifelse(n5$NVHIB>2,1,0)
n5$IPNVSAR=ifelse(n5$NVSAR>0,1,0)
n5$IPNVPaR=ifelse(n5$NVPaR>0,1,0)
n5$IPNVAMA=ifelse(n5$NVAMA>0,1,0)
n5$IPMaEdu1=ifelse(n5$MaEdu>1,1,0)
n5$IPMaEdu2=ifelse(n5$MaEdu>2,1,0)
n5$IPNMSD=ifelse(n5$NMSD==1,1,0)
n5$IPNMCCD=ifelse(n5$NMCCD>0,1,0)
n5$IPAlumbra=ifelse(n5$Alumbra==1,1,0)
n5$IPAgua=ifelse(n5$Agua==1,1,0)
n5$IPAgRTC=ifelse(n5$AgRT==2 | n5$AgRT==3,1,0)
n5$IPAgRTH=ifelse(n5$AgRT==1 | n5$AgRT==3,1,0)
n5$IPExcre1=ifelse(n5$Excre==1,1,0)
n5$IPExcre3=ifelse(n5$Excre==3 & n5$LetOK==1,1,0)
n5$IPTCoc1=ifelse(n5$TCoc==1,1,0)
n5$IPTCoc2=ifelse(n5$TCoc==2,1,0)
n5$IPCUI=ifelse(n5$CUI==1,1,0)
n5$IPNBI=ifelse(n5$NBI>0,1,0)
n5$IPNJtos=ifelse(n5$NJtos>0,1,0)
n5$IPSIS=ifelse(n5$SIS==1,1,0)
n5$IPNBAA=ifelse(n5$NBAA>0,1,0)
n5$IPNPapB=ifelse(n5$NPapB>0,1,0)
n5$IPNBPIN=ifelse(n5$NBPIN>0,1,0)
n5$IPNBWW=ifelse(n5$NBWW>0,1,0)
n5$IPNBCOI=ifelse(n5$NBCOI>0,1,0)
n5$IPNBFAR=ifelse(n5$NBFAR>0,1,0)
n5$IPNBTBC=ifelse(n5$NBTBC>0,1,0)
n5$IPNBCOP=ifelse(n5$NBCOP>0,1,0)
n5$IPNBPRV=ifelse(n5$NBPRV>0,1,0)
n5$IPNBVAL=ifelse(n5$NBVAL>0,1,0)
# variables de evaluación
n5$EJtos = factor(ifelse(n5$NJtos>0,1,0) + 2*ifelse(n5$NBI>0,1,0),labels=c("E-J-","E-J+","E+J-","E+J+"))
n5$ESIS  = factor(ifelse(!is.na(n5$SIS) & !is.na(n5$NBI),ifelse(n5$SIS==1,1,0) + 2*ifelse(n5$NBI>0,1,0),0),labels=c("E-S-","E-S+","E+S-","E+S+"))
n5$EPIN  = factor(ifelse(n5$NBPIN>0,1,0) + 2*ifelse(n5$NBI>0,1,0),labels=c("E-P-","E-P+","E+P-","E+P+"))
n5$EVL   = factor(ifelse(n5$NBVAL>0,1,0) + 2*ifelse(n5$NBI>0,1,0),labels=c("E-V-","E-V+","E+V-","E+V+"))
# variables de modelamiento
n5$ESI8H9 = as.numeric(n5$Estrato)*100+n5$INEIE8*10+n5$HV109
n5$cagua = ifelse(is.na(n5$Agua),NA,ifelse(n5$Agua<=2,1,0))
n5$Edad2 = n5$EdadM^2
n5$FALT2 = n5$FALT^2
n5$RegNat = factor(ifelse(n5$Estrato==1&n5$Estrato==2,1,ifelse(n5$Estrato==5,3,2)))
n5$ALMPM = sqrt((n5$FLAT+(12+2/60+45.59/60/60))^2+(n5$FLON+(77+1/60+49.97/60/60))^2)
n5$csexo = 2 - n5$Sexo
n5$fHV109 = factor(ifelse(n5$HV109==0,7,n5$HV109))
n5$fINEIE8 = factor(n5$INEIE8)
n5$fINEIE4 = factor(n5$INEIE4)
n5$fINEIE2 = factor(n5$INEIE2)
n5$nHV109 = ifelse(n5$HV109==1,3,ifelse(n5$HV109==2,6,ifelse(n5$HV109==3,8.5,ifelse(n5$HV109==4,11,ifelse(n5$HV109==5,15,0)))))
n5$nINEIE8 = ifelse(n5$INEIE8==1,(5000-100)/2,ifelse(n5$INEIE8==2,(100-20)/2,ifelse(n5$INEIE8==3,(20-10)/2,ifelse(n5$INEIE8==4,(10-4)/2,ifelse(n5$INEIE8==5,(4-0.4)/2,ifelse(n5$INEIE8==6,(0.4-0.2)/2,ifelse(n5$INEIE8==7,(0.2-0.1)/2,ifelse(n5$INEIE8==8,(0.1-0)/2,0))))))))
n5$lINEIE8 = log(n5$nINEIE8)
n5$fINEIER = factor(n5$INEIER)
n5$MaEdA2 = n5$MaEdA^2
n5$gMaEdad = 5*as.integer(n5$MaEdA/5)
n5$nVertOc = 2 - n5$VertOc 	# 1 occidente, 0 oriente
n5$lCampoD=log(n5$CampoD)
# variables de conglomerados y ponderación
table(n5$N5AliD,n5$Estrato)
n5<-subset(n5,n5$Salida<14 & n5$StatusM!=15 & n5$EdadM>=0 & n5$EdadM<60)
n5$ExFNi=n5$SWt*2472049/sum(n5$SWt)
n5$SWtNi=n5$SWt*length(n5$SWt)/sum(n5$SWt)
n5$gedad=factor(ifelse(n5$EdadM<12,6*as.integer(n5$EdadM/6),12*as.integer(n5$EdadM/12)),
	labels=c("0-5m","6-11m","12-23m","24-35m","36-47m","48-59m"))
# archivo de visitas
#v<-read.csv("G:/Temp/Draft/MONIN/V.csv") 
x<-aggregate(x=n5[c("SWt","INEIE8","Salida")],
	by=list(n5$Estrato,n5$Cong1,n5$Cong2),FUN=max,na.rm=TRUE)
w<-aggregate(x=n5[c("N5Nut", "N5NutD", "N5DTED","N5DTE","N5OBED","N5OBE","N5AneD","N5Ane",
	"N5AliD","N5CCDD","N5CC0D","N5Ali","N5San","N5SaY","N5PIP","N5CCD","N5CC0","N5ATE","N5CPR","N5VIVD","N5VIV","N5COM","N5PRO")],
	by=list(n5$Estrato,n5$Cong1,n5$Cong2),FUN=sum,na.rm=TRUE)
v<-aggregate(x=n5[c("IdM")],
	by=list(n5$Estrato,n5$Cong1,n5$Cong2),FUN=length)
w<-merge(x=w,y=x,by=c("Group.1", "Group.2", "Group.3"),all.y=TRUE)
v5<-merge(x=v,y=w,by=c("Group.1", "Group.2", "Group.3"),all.y=TRUE)
v5$Estrato=v5$Group.1; v5$Cong1=v5$Group.2; v5$Cong2=v5$Group.3
v5$SWtVi=v5$SWt*length(v5$SWt)/sum(v5$SWt)
#s<-as.svrepdesign(s, type="bootstrap", replicates=1000)
# subgrupo 0-2
n3<-subset(n5,n5$EdadM<36)
n3$ExFNi=n3$SWt*2472049*(3/5)/sum(n3$SWt)
n3$SWtNi=n3$SWt*length(n3$SWt)/sum(n3$SWt)
# diseños muestrales
p5<-svydesign(id=~Cong1+Cong2, strata=~Estrato, weights=~SWtNi, data=n5)
p3<-svydesign(id=~Cong1+Cong2, strata=~Estrato, weights=~SWtNi, data=n3)
#p<-svydesign(id=~Cong1, strata=~Estrato, weights=~SWt2, data=n)
#b<-as.svrepdesign(p, type="bootstrap", replicates=1000)
s5<-svydesign(id=~Cong1+Cong2, strata=~Estrato, weights=~SWtVi, data=v5)
v=1; w=1; x=1; y=1
# almacenamiento temporal
save.image( "G:/Temp/DRAFT/MONIN/M.RData" )
load( "G:/Temp/DRAFT/MONIN/M.RData" )
summ(n5); summ(v5); summary(p5); summary(s5)
summ(n3); summary(p3)
date()


# ======================================================================================================================== #
# Tabulaciones y graficaciones principales


t<-"===== MONIN Nov 2007 - Feb 2008 (distribución de conglomerados)";t;date()
#by(tc[,c("NVivs","TotVivO","NClasV","NEleDes","NNoEle","NEleNo","NCoEViv","NNi04","NNi02","NMEF","NGes")], tc$MONIN5, sum,na.rm=TRUE)
#by(data.frame(tc$NVivs,tc$TotVivO,tc$NClasV,tc$NEleDes,tc$NNoEle,tc$NEleNo,tc$NCoEViv,tc$NNi04,tc$NNi02,tc$NMEF,tc$NGes), tc$MONIN5, sum,na.rm=TRUE)
table(tc$StatusM,tc$MONIN5)
table(tc$INEIE8,tc$MONIN5)
FSumX ( tc$one, tc$MONIN5,     "Conglomerados (Marco INEI)         " )
FSumX ( tc$NVivs, tc$MONIN5,   "Viviendas (Marco INEI)             " )
FSumX ( tc$NVivsL, tc$MONIN5,  "Viviendas (Marcos Locales)         " )
FSumX ( tc$TotVivO, tc$MONIN5, "Viviendas Enumeradas               " )
FSumX ( tc$NClasV, tc$MONIN5,  "Viviendas Clasificadas             " )
FSumX ( tc$NEleDes, tc$MONIN5, "Viviendas Desconocidas             " )
FSumX ( tc$NNoEle, tc$MONIN5,  "Hogares No Elegibles               " )
FSumX ( tc$NEleNo, tc$MONIN5,  "Hogares Elegibles Excluidos        " )
FSumX ( tc$NEleRe, tc$MONIN5,  "Hogares Elegibles que Rechazaron   " )
FSumX ( tc$NCoEViv, tc$MONIN5, "Hogares Entrevistados              " )
FSumX ( tc$NNi04, tc$MONIN5,   "Niños 0-4 años Censados            " )
FSumX ( tc$NNi02, tc$MONIN5,   "Niños 0-2 años Censados            " )
FSumX ( tc$NMEF, tc$MONIN5,    "Mujeres 15-49 años Censadas        " )
FSumX ( tc$NGes, tc$MONIN5,    "Mujeres Gestantes Censadas         " )
FSumX ( tc$N04HA, tc$MONIN5,   "Niños 0-59 meses con Antropometría " )
FSumX ( tc$N04HB, tc$MONIN5,   "Niños 0-59 meses con Hemoglobina   " )
FSumX ( n5$one, n5$Estrato,    "Niños 0-59 meses Entrevistados     " )
FSumX ( n5$N5DTED, n5$Estrato, "Niños 0-59 meses con Antropometría " )
FSumX ( n5$N5AneD, n5$Estrato, "Niños 0-59 meses con Hemoglobina   " )
#svymean(~ZHAZ,p5,na.rm=TRUE,deff=TRUE)
#svyby(~ZHAZ,~Estrato,design=p5,svymean,na.rm=TRUE,deff=TRUE)
#svyratio(~N5DTE, ~N5DTED, design=s5)
#svyby(~N5DTE,by=~Estrato,denominator=~N5DTED,design=s5,svyratio)


t<-"===== MONIN Nov 2007 - Feb 2008 (indicadores principales)";t;date()
FPack1 ( xby=~Estrato, fplot=TRUE )


t<-"===== MONIN Nov 2007 - Feb 2008 (indicadores de Crecer)";t;date()
table(n5$EtapaC)
table(n5$Crecer)
svymean(~ifelse(EtapaC>0,1,0),p5,na.rm=TRUE,deff=FALSE)
svyby(~ifelse(EtapaC>0,1,0),~Estrato,design=p5,svymean,na.rm=TRUE,deff=FALSE)
t<-"===== Niños 0-59 meses, Indicadores Crecer - Decreto Supremo 080";t
#incidencia de probreza y pobreza extrema
#prevalencia de desnutrición en menores de cinco años
#proporción de menores de 6 meses con lactancia materna exclusiva 
#prevalencia de anemia en menores de 36 meses
#incidencia de infección respiratoria aguda en menores de 36 meses
#incidencia de enfermedad diarréica aguda en menores de 36 meses
##incidencia de bajo peso al nacer
##tasa de mortalidad neonatal por mil nacidos vivos
#incremento de identificación (CUI)
FPRwM2 ( ~ifelse(NBI>0,1,0), p5,                           "NBI>0                ", xby=~Crecer )
FPRwM2 ( ~ifelse(ZHAZ<(-2),1,0), p5,                       "TE <-2zs OMS 2006    ", xby=~Crecer )
FPRwM2 ( ~N5LME, subset(p3,EdadM<6),                       "Lact. M. Excl. [<6m] ", xby=~Crecer )
FPRwM2 ( ~ifelse(HbCORR<11,1,0), p3,                       "Hbc<11 [<36m]        ", xby=~Crecer )
FPRwR2 ( ~NumeIR, ~DenomIR, p3,                            "Inc. Enf. Resp. Ag.  ", xby=~Crecer )
FPRwR2 ( ~NumeID, ~DenomID, p3,                            "Inc. Enf. Diarr. Ag. ", xby=~Crecer )
FPRwM2 ( ~ifelse(CUI==1,1,0), p5,                          "CUI                  ", xby=~Crecer )
t<-"===== Niños 0-59 meses, Indicadores Crecer - Presentación de Estrategia";t
#.desnutrición infantil
#..dieta pobre en cantidad y calidad
#...hábitos alimenticios adecuados
#....educación alimentaria nutricional
#....promoción de productos locales
#...complementación nutricional y alimentaria
#....provisión de alimentos fortificados
#....seguridad alimentaria
#....transferencias condicionadas
#...suplemento micronutrientes
#....provisión de micronutrientes
#..enfermedades diarreicas agudas +
#..infecciones respiratorias aguda +
#..enfermedades prevalentes regionales
#...agua segura
#....saneamiento básico
#...prácticas saludables
#....cocinas mejoradas
#....letrinas adecuadas
#....alfabetización
#...atención integral +
#..bajo peso al nacer
#...atención prenatal de calidad +
#...atención integral de la mujer
#....aseguramiento de la salud
#ciudadanía
#monitoreo y evaluación
FPRwM2 ( ~ifelse(ConLM==2,0,ConLM), subset(p3,EdadM<24),   "Lact. Materna   [<24]", xby=~Crecer )
FPRwM2 ( ~ifelse(VComio>=3,1,0), subset(p3,EdadM<24),      "Comió 3+ SS     [<24]", xby=~Crecer )
FPRwM2 ( ~MasGra, subset(p3,EdadM<24),                     "Grasa adicional [<24]", xby=~Crecer )
FPRwM2 ( ~ifelse(NMSD==1,1,0), p5,                         "Sesión demostrativa  ", xby=~Crecer )
FPRwM2 ( ~ifelse(NBAA>0,1,0), subset(p5,NBI>0),            "Asistencia Alim.[NBI]", xby=~Crecer )
FPRwM2 ( ~ifelse(NJtos>0,1,0), subset(p5,NBI>0),           "Juntos          [NBI]", xby=~Crecer )
FPRwM2 ( ~ifelse(VitA==1,1,0), subset(p3,EdadM<24),        "Supl. Vit. A    [<24]", xby=~Crecer )
FPRwM2 ( ~ifelse(VFe>5,1,0), subset(p3,EdadM<24),          "Supl. Fe        [<24]", xby=~Crecer )
FPRwM2 ( ~ifelse(Alumbra==1,1,0), p5,                      "Electricidad         ", xby=~Crecer )
FPRwM2 ( ~ifelse(Agua==1,1,0), p5,                         "Agua en Red dentro   ", xby=~Crecer )
FPRwM2 ( ~ifelse(AgRT==2|AgRT==3,1,0), p5,                 "Recibe Agua Clorada  ", xby=~Crecer )
FPRwM2 ( ~ifelse(AgRT==1|AgRT==3,1,0), p5,                 "Consume Agua Hervida ", xby=~Crecer )
FPRwM2 ( ~ifelse(Excre==1,1,0), p5,                        "Desagüe en Red dentro", xby=~Crecer )
FPRwM2 ( ~ifelse(Excre>1&Excre<5,1,0), p5,                 "Letrina/Pozo/Red f   ", xby=~Crecer )
FPRwM2 ( ~ifelse(Excre==3&LetOK==1,1,0), p5,               "Letrina adecuada      ", xby=~Crecer )
FPRwM2 ( ~ifelse(TCoc==1,1,0), p5,                         "Combustible Elec/Gas ", xby=~Crecer )
FPRwM2 ( ~ifelse(TCoc==2,1,0), p5,                         "Cocina Asesorada     ", xby=~Crecer )
FPRwM2 ( ~ifelse(MaEdu>1,1,0), p5,                         "Madre Alfabeta       ", xby=~Crecer )
FPRwM2 ( ~N5CC0, p5,                                       "CRED al dia 1er año  ", xby=~Crecer )
FPRwM2 ( ~N5PIP, p5,                                       "PAI al dia           ", xby=~Crecer )
FPRwM2 ( ~ifelse(SIS==1,1,0), subset(p5,NBI>0),            "SIS             [NBI]", xby=~Crecer )
# Crecer versus paquete estándar
FPack1 ( xby=~Crecer, fplot=TRUE )
FPTbMR ( xvec=c(~IPNBPIN,~IPNBVAL),
	 xdes=list(subset(p5,NBI>0),subset(p5,NBI>0)),
	 xlab=c("PIN[P]","VL[P]"), 
	 xby=~Crecer, fratio=FALSE, fplot=FALSE ) 


t<-"===== MONIN Nov 2007 - Feb 2008 (evaluación de intervenciones)";t;date()
FPack1 ( xby=~EJtos, fplot=FALSE )
FPack1 ( xby=~ESIS, fplot=FALSE )
FPTbMR ( xvec=c(~IPCUI,~IPNBI,~IPNJtos,~IPNBAA,~IPNPapB,~IPNBAA),
	 xdes=list(p5,p5,subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI>0),subset(p5,NBI==0)),
	 xlab=c("CUI","P=NBI>0","Juntos [P]","AA [P]","Pap [P]","AA[noP]"), 
	 xby=~ESIS, fratio=FALSE, fplot=FALSE ) 
FPack1 ( xby=~EPIN, fplot=FALSE )
FPack1 ( xby=~EVL, fplot=FALSE )


t<-"===== MONIN Nov 2007 - Feb 2008 (indicadores de DI)";t;date()
t<-"===== Niños 0-59 meses, Indicadores Principales - Prevalencia [e.s.], est. media en niños";t
FPack1 ( xby=~Estrato, fplot=TRUE )
t<-"===== Niños 0-35 meses, Indicadores Principales - Prevalencia [e.s.], est. razón en hogares";t
FPRwM2 ( ~N5Nut, p3,                 "Buena Nutrición      " )
FPRwM2 ( ~N5Ali, p3,                 "Buena Alimentación   " )
FPRwM2 ( ~N5SaY, p3,                 "Buena Salud          " )
FPRwM2 ( ~N5ATE, p3,                 "Buena Atención       " )
FPRwM2 ( ~N5CPR, p3,                 "Cuidadores Preparados" )
FPRwM2 ( ~N5VIV, p3,                 "Vivienda Saludable   " )
FPRwM2 ( ~N5COM, p3,                 "Comunidad Saludable  " )
t<-"===== Niños 0-35 meses, Indicadores Componentes - Prevalencia [e.s.], est. media en niños";t
FPRwM2 ( ~N5DTE, p3,                 "Talla-Edad<-2z 2006  " )
FPRwM2 ( ~N5ADE, p3,                 "Peso-Talla<-2z 2006  " )
FPRwM2 ( ~N5OBE, p3,                 "Obesidad>2z 2006     " )
FPRwM2 ( ~N5Ane, p3,                 "Anemia (Hb(c)<11)    " )
#svymean(~N5Ane,p3,na.rm=TRUE,deff=TRUE)
#svyby(~N5Ane,~Estrato,design=subset(p3,Estrato!=4),svymean,na.rm=TRUE,deff=TRUE)
FPRwM2 ( ~ifelse(SalRx==1,1,0), p3,  "Sal Iodada (1=15+)   " )
FPRwM2 ( ~ifelse(ConLM==2,0,ConLM), subset(p3,EdadM<24),         "Lact. Materna [<24m] " )
FPRwM2 ( ~N5LME, subset(p3,EdadM<6),              "L.M. Exclusiva [<6m] " )
FPRwM2 ( ~ifelse(AliOK==2,0,AliOK), subset(p,EdadM>=6),         "Alimentos No LM [6+m]" )
FPRwM2 ( ~N5Ali, subset(p3,N5SaY==0),"Alim. en Enfermedad  " )
FPRwM2 ( ~ifelse(SupOK==1,1,0), p3,  "Suplementos          " )
FPRwM2 ( ~ifelse(NVVTA>0,1,0), p3,   "Suplemento Vitamina A" )
FPRwM2 ( ~ifelse(PDia2w==1,1,0), p3, "Diarrea (2 sem)      " )
FPRwM2 ( ~ifelse(PTos2w==1,1,0), p3, "Tos (2 sem)          " )
FPRwM2 ( ~ifelse(PFie2w==1,1,0), p3, "Fiebre (2 sem)       " )
FPRwM2 ( ~NEDY, p3,                  "Enf. Diarreica ayer  " )
FPRwM2 ( ~NERY, p3,                  "Enf. Resp. ayer      " )
FPRwR2 ( ~NED7, ~DenomP, p3,         "Prev.Enf.Diar. (dia) " )
FPRwR2 ( ~NER7, ~DenomP, p3,         "Prev.Enf.Resp. (dia) " )
FPRwR2 ( ~NumeID, ~DenomID, p3,      "Inc.Enf.Diar. (dia)  " )
FPRwR2 ( ~NumeIR, ~DenomIR, p3,      "Inc.Enf.Resp. (dia)  " )
FPRwR2 ( ~NumeRD, ~DenomRD, p3,      "Rec.Enf.Diar. (dia)  " )
FPRwR2 ( ~NumeRR, ~DenomRR, p3,      "Rec.Enf.Resp. (dia)  " )
FPRwM2 ( ~N5PIP, p3,                 "Vacunado al dia      " )
FPRwM2 ( ~N5CC0, p3,                 "CRED 1er. Año al dia " )
FPRwM2 ( ~ifelse(NVBCG>0,1,0), p3,   "Vacuna BCG           " )
FPRwM2 ( ~ifelse(NVPOL>2,1,0), p3,   "Vacuna AntiPolio     " )
FPRwM2 ( ~ifelse(NVDPT>2,1,0), p3,   "Vacuna DPT           " )
FPRwM2 ( ~ifelse(NVHEP>2,1,0), p3,   "Vacuna Hepatitis     " )
FPRwM2 ( ~ifelse(NVHIB>2,1,0), p3,   "Vacuna Haemophilus   " )
FPRwM2 ( ~ifelse(NVSAR>0,1,0), p3,   "Vacuna Sarampión     " )
FPRwM2 ( ~ifelse(NVPaR>0,1,0), p3,   "Vacuna RubeolaPaperas" )
FPRwM2 ( ~ifelse(NVAMA>0,1,0), p3,   "Vacuna FiebreAmarilla" )
FPRwM2 ( ~MaPte, p3,                 "Madre Presente       " )
FPRwM2 ( ~ifelse(MaEdu>1,1,0), p3,   "Madre Lee y Escribe  " )
FPRwM2 ( ~ifelse(MaEdu>2,1,0), p3,   "Madre con Primaria   " )
FPRwM2 ( ~ifelse(NMSD==1,1,0), p3,   "Asistió a SD hace <6m" )
FPRwM2 ( ~ifelse(NMCCD>0,1,0), p3,   "CRED con 1+ contenido" )
FPRwM2 ( ~ifelse(Agua==1,1,0), p3,   "Agua en Red Domic.   " )
FPRwM2 ( ~ifelse(Excre==1,1,0), p3,  "Desagüe en Red Domic." )
FPRwM2 ( ~LetOK, p3,                 "Letrina Adecuada[L]  " )
svymean(~LetOK,p3,na.rm=TRUE,deff=TRUE)
svyby(~LetOK,~Estrato,design=subset(p3,Estrato!=1),svymean,na.rm=TRUE,deff=TRUE)
FPRwM2 ( ~ifelse(TCoc==1,1,0), p3,   "Cocina Elec. o Gas   " )
FPRwM2 ( ~ifelse(TCoc==2,1,0), subset(p3,NBI>0),  "Cocina Mejorada [NBI]" )
FPRwM2 ( ~ifelse(AgRT==2|AgRT==3,1,0), p3,        "Reciben Agua Clorada " )
FPRwM2 ( ~ifelse(AgRT==1|AgRT==3,1,0), p3,        "Beben Agua Hervida   " )
FPRwM2 ( ~ifelse(CUI==1,1,0), p3,    "Con CUI              " )
FPRwM2 ( ~ifelse(NBI>0,1,0), p3,     "NBI > 0 (pobreza)    " )
FPRwM2 ( ~ifelse(NBAA>0,1,0), subset(p3,NBI>0),   "A.Alimentaria [NBI]  " )
FPRwM2 ( ~ifelse(NJtos>0,1,0), subset(p3,NBI>0),  "Subs. Juntos [NBI]   " )
FPRwM2 ( ~ifelse(NPapB>0,1,0), subset(p3,NBI>0),  "Papilla [NBI]        " )
FPRwM2 ( ~ifelse(SIS==1,1,0), subset(p3,NBI>0),   "Seguro SIS [NBI]     " )


t<-"===== MONIN Nov 2007 - Feb 2008 (modelos descriptivos, sin ajustar diseño)";t;date()
DescMod1 ( n5$N5Nut, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5Ali, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5SaY, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5ATE, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5CPR, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5VIV, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5COM, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5DTE, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5ADE, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5OBE, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5Ane, n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$SalRx==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$ConLM==2,0,n5$ConLM), n5, n5$EdadM<24 )
DescMod1 ( n5$N5LME, n5, n5$EdadM<6 )
DescMod1 ( ifelse(n5$AliOK==2,0,n5$AliOK), n5, n5$EdadM>=6 )
DescMod1 ( n5$N5Ali, n5, n5$N5SaY==0 )
DescMod1 ( ifelse(n5$SupOK==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVVTA>0,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$PDia2w==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$PTos2w==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$PFie2w==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5PIP, n5, !is.na(n5$IdM) )
DescMod1 ( n5$N5CC0, n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVBCG>0,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVPOL>2,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVDPT>2,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVHEP>2,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVHIB>2,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVSAR>0,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVPaR>0,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NVAMA>0,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( n5$MaPte, n5, !is.na(n5$IdM) )
DescMod2 ( ifelse(n5$MaEdu>1,1,0), n5, !is.na(n5$IdM) )
DescMod2 ( ifelse(n5$MaEdu>2,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NMSD==1,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$NMCCD>0,1,0), n5, !is.na(n5$IdM) )
DescMod4 ( ifelse(n5$Agua==1,1,0), n5, !is.na(n5$IdM) )
DescMod4 ( ifelse(n5$Excre==1,1,0), n5, !is.na(n5$IdM) )
#LetOK
DescMod1 ( ifelse(n5$TCoc==1,1,0), n5, !is.na(n5$IdM) )
#ifelse(TCoc==2,1,0)
DescMod1 ( ifelse(n5$AgRT==2|n5$AgRT==3,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$AgRT==1|n5$AgRT==3,1,0), n5, !is.na(n5$IdM) )
DescMod1 ( ifelse(n5$CUI==1,1,0), n5, !is.na(n5$IdM) )
DescMod3 ( ifelse(n5$NBI>0,1,0), n5, !is.na(n5$IdM) )
DescMod3 ( ifelse(n5$NBAA>0,1,0), n5, n5$NBI>0 )
DescMod3 ( ifelse(n5$NJtos>0,1,0), n5, n5$NBI>0 )
DescMod3 ( ifelse(n5$NPapB>0,1,0), n5, n5$NBI>0 )
DescMod3 ( ifelse(n5$SIS==1,1,0), n5, n5$NBI>0 )


t<-"===== MONIN Nov 2007 - Feb 2008 (hipótesis adicionales)";t;date()
t<-"===== Niños 0-59 meses, Cobertura de Vaso de Leche, Pobres vs. No Pobres";t
n = n5
n$IPNBI=factor(n$IPNBI,labels=c("NBI=0","NBI>0"))
p<-svydesign(id=~Cong1+Cong2, strata=~Estrato, weights=~SWtNi, data=n)
# Nacional: Cobertura de Vaso de Leche según Pobreza
svyby(~IPNBVAL,~IPNBI,design=p,svymean,na.rm=TRUE,deff=FALSE)
svychisq(~IPNBI+IPNBVAL,p)
# Por dominios
svyby(~IPNBVAL,~IPNBI+Estrato,design=p,svymean,na.rm=TRUE,deff=FALSE)
svyby(~IPNBVAL,~IPNBI+Crecer,design=p,svymean,na.rm=TRUE,deff=FALSE)
# Análisis Combinado, dependiente: cobertura VL, predictores: NB, Estrato, Crecer
> # triaje aproximado
summary(stepAIC(glm(IPNBVAL~IPNBI*Estrato*Crecer,
	data=n,family=binomial(link="logit"),na.action=na.exclude),trace=FALSE))
# confirmación ajustando a diseño muestral
r<-svyglm(IPNBVAL~IPNBI+Estrato+Crecer,p,family=binomial(link="logit"))
summary(r)
regTermTest(r,~IPNBI)
regTermTest(r,~Estrato)
regTermTest(r,~Crecer)
# conclusión: cobertura de VL depende separadamente de NBI, Estrato y Crecer;
# las diferencias por NBI no son significativamente distintas por Estrato o Crecer;
# la diferencia por NBI es pequeña, OR=1.3;
# las mayores diferencias son por Sierra Rural, OR=5.1, y Crecer, OR=3.4;
# (NBI, Sierra Rural y Crecer tienen mayores coberturas de VL).


t<-"===== Niños 0-59 meses, Población con Superposición de Hierro";t
n = n5
n$EBFE=factor(8*ifelse(n$IPVFe==1,1,0)+4*ifelse(n$EBVAL>0,1,0)+2*ifelse(n$EBPIN>0,1,0)+ifelse(n$EBWW>0,1,0),
	labels=c("Nada","WW","PIN","VL","PIN+VL","SF","SF+WW","SF+PIN","SF+VL","SF+PIN+VL"))
# 0   1   2   4   6   8   9  10  12  14 
# 3: "PIN+WW", 5:"VL+WW", 7:"PIN+VL+WW", 12:"SF+PIN+WW", 13:"SF+VL+WW", 15:,"SF+VL+PIN+WW"
n$NUFE=ifelse(n$IPVFe==1,1,0)+ifelse(n$EBVAL>0,1,0)+ifelse(n$EBPIN>0,1,0)+ifelse(n$EBWW>0,1,0)
#n$EBPIN=factor(ifelse(n$EBPIN>0,1,0),labels=c("NoPIN","PIN"))
#n$EBWW =factor(ifelse(n$EBWW>0,1,0),labels=c("NoWW","WW"))
#n$EBVAL=factor(ifelse(n$EBVAL>0,1,0),labels=c("NoVL","VL"))
#n$IPVFe=factor(n$IPVFe,labels=c("NoSFe","SFe"))
p<-svydesign(id=~Cong1+Cong2, strata=~Estrato, weights=~SWtNi, data=n)
# Distribución de todos los niños 0-59, combinaciones de fuentes de hierro
svytable(~EBFE,p,Ntotal=1)
# Distribución de todos los niños 0-59, número de fuentes de hierro
svytable(~NUFE,p,Ntotal=1)
# svytable(~EBPIN+EBVAL+EBWW+IPVFe,p,Ntotal=1)
# Distribución solo de los niños 0-59 que reciben alguna de las fuentes
# svytable(~EBPIN+EBVAL+EBWW+IPVFe,subset(p,n$EBPIN=="PIN"|n$EBWW=="WW"|n$IPVFe=="SFe"),Ntotal=1)
# 10% reciben mas de una fuente (incluye 9% declarado como beneficiario VL siendo 0-59)


date()
sink()
q("no")


# ======================================================================================================================== #
