1 Muestreo aleatorio simple con reemplazo

Considere una población de tamaño \(N\) donde se realizan \(m=3\) experimentos independientes. En cada uno, cada unidad de las \(N\) tiene probabilidad de ser seleccionada igual a \(1/N\). Dé la función de probabilidad de \(n_S\), es decir, de la variable asociada al número de elementos distintos en la muestra.

Denotemos a nuestra población con la letra \(U\). Para este ejercicio haré uso de la forma de un “Muestreo aleatorio simple con reemplazo”, el cual se ve de la siguiente forma:
Sea \(a_k\) el número de veces en que el elemento \(k\) (elemento de la población U) aparece en la muestra ordenada \(os=\{k_1,k_2,k_3\}\), donde sabemos por la forma del “Muestreo aleatorio simple con reemplazo” que la variable \(a_k\) tiene distribución \(Bin(m=3,\frac{1}{N})\). Una forma de notar el número de elementos distintos (o sea \(n_S\)) es viendo que esto se puede escribir de la forma: \(n_S=\sum_{k\in U}a_k\), o bien, también podemos definir una función indicadora como se ve adelante: \[\mathbb{I}(a_k):=\left\{ 1 \text{, si, } a_k\geq 0 \atop 0 \text{, si, } a_k= 0\right.\] Entonces denotamos a \(n_S\) como \(n_S = \sum_{k\in U}\mathbb{I}(a_k)\), de aquí la variable aleatoria \(n_S\) toma valores entre 1 y 3. Ahora definamos \(n_S(m)\) al número de elementos distintos obtenidos en el “Muestreo aleatorio simple con reemplazo” con m experimentos independientes de la población U, con \(m\in \{1,2,3\}\).
\(P(n_S(3)=n)=\\P(n_S(2)=n-1)\cdot P(\text{Seleccionar en el tercer experimento un elemento no seleccionado}|n_S(2)=n-1)+P(n_S(2)=n)\cdot P(\text{Seleccionar en el tercer experimento un elemento ya seleccionado} | n_S(2)=n)=\\ P(n_S(2)=n-1)\frac{N-n+1}{N} + P(n_S(2)=n)\frac{n}{N} .........(I)\)
Ahora notemos el hecho de que un “Número de Stirling de segundo tipo” \(\left\{n\atop 2 \right\}\) representa el número de formas de partir un conjunto de 2 elementos en n conjuntos no vacíos, además de que hay que notar que hay \(\frac{N!}{(N-n+1)!}\) maneras de elegir \(n-1\) elementos de la población U (aquí recordemos que el orden importa, dado que estamos en experimentos independientes y que estos mismos dependen del orden en que se hagan) y además donde en total tenemos \(N^{2}\) maneras de seleccionar 2 elementos cualesquiera de U; por lo que la probabilidad \(P(n_S(2)=n-1)\) es igual a \(\frac{N!}{(N-n+1)!N^{2}}\left\{n\atop 2 \right\}\). Este mismo análisis se hace con \(P(n_S(2)=n)\).
Dicho de esta manera, nuesta expresión (I) la podemos llevar como:
\(P(n_S(2)=n-1)\frac{N-n+1}{N} + P(n_S(2)=n)\frac{n}{N}=\frac{N!}{(N-n+1)!N^{2}}\left\{n-1\atop 2 \right\}\frac{N-n+1}{N} + \frac{N!}{(N-n)!N^{2}}\left\{n\atop 2 \right\}\frac{n}{N}\)
Ahora simplifiquemos esta expresión:
\(\frac{N!}{(N-n+1)!N^{2}}\left\{n-1\atop 2 \right\}\frac{N-n+1}{N} + \frac{N!}{(N-n)!N^{2}}\left\{n\atop 2 \right\}\frac{n}{N}=\frac{N!}{(N-n)!N^3} (\left\{n-1\atop 2 \right\}+u\left\{n\atop 2 \right\})=\frac{N!}{(N-n)!N^3}\left\{n\atop 3 \right\}\)
Por lo tanto tenemos que \[p(n_s(3)=n):=\left\{\frac{N!}{(N-n)!N^3}\left\{n\atop 3 \right\} \quad si \quad n\in\{1,2,3\} \atop 0 \quad e.o.c \right.\]

\[ p(n_s(3)=n) = \left\{\begin{array}{lr} \frac{1}{N^2}*1, & \text{for } n = 1\\ \frac{N-1}{N^2}*3, & \text{for } n=2\\ \frac{(N-1)(N-2)}{N^2}*1, & \text{for } n=3\\ 0, & \text{e.o.c } \end{array}\right\} \] \(\blacksquare\)

2 Esquemas de selección sin reemplazo

Considere el siguiente esquema de selección para obtener una muestra de tamaño \(n=2\).

  1. Del conjunto de \(N\) unidades se seleciona de forma aleatoria una unidad, en esta primera selección cada unidad tiene la misma probabilidad de ser seleccionada.
  2. Una vez seleccionada una unidad, se selecciona otra unidad del conjunto de las \(N-1\) unidades restantes (sin contar la seleccionada en la primera selección). En esta selección cada unidad de las \(N-1\) restantes tiene la misma probabilidad de ser seleccionada.

Indique cuáles son las probabilidades de primer y segundo orden si se usa un diseño de muestreo con este esquema de selección.

Notemos que el esquema de selección tenemos un orden en la selcción, ya que se hace la distinción de probabilidades en cada elección de los elementos de nuestra población \(U=\{u_1,...,u_N\}\). Al importar el orden, tenemos el siguiente número de posibles elecciones: \[\frac{N!}{(N-2)!}\] Es decir, tenemos \(N(N-1)\) posibles parejas de elementos.
En este sentido, también notemos que por hipótesis, tenemos que la probabilidad de elegir el primer elemento es de \(\frac{1}{N}\), de aquí tenemos el caso en que si agarramos una muestra \(s_i\), con \(i\in \{1,...,N(N-1)\}\), tenemos que elegir un elemento \(k\) tiene probabilidad de \(\frac{1}{N}\) de ser seleccionada, ahora veamos el hecho de que en nuestro conjunto de todas las muestras posibles con nuestro esquema de selección, tenemos las siguientes dos muestras: \((u_i,u_j)\) y \((u_j, u_i)\), por lo tanto, tenemos que: \[\pi_k = \frac{2}{N}\] Del mismo modo, ahora obtenemos que la selección del primer elemento es independiente a la selección del segundo, de este modo la propabilidad de elegir al elemento \(k\) y \(j\) es \(\frac{1}{N(N-1)}\), siendo de la misma forma, tenemos a las muestras \((u_i,u_j)\) y \((u_j, u_i)\). Por lo tanto la probabilidad de segundo orden es: \[\pi_{kj}=\frac{2}{N(N-1)}\] \(\blacksquare\)

3 Muestreo con probabilidad proporcional al tamaño sin reemplazo

Considerando un diseño de muestreo con probabilidad proporcional al tamaño sin reemplazo y tamaño de muestra igual a 12, se seleccionó una muestra de una población de tamaño \(N=284\) con el objetivo de estimar el total de la variable \(REV84\). La variable tamaño con la que se construyeron las probabilidades de inclusión es la variable \(P75\) (el total de la variable \(P75\) en la población es 8,182). La información recolectada de la variable \(REV84\) en la muestra seleccionada es la siguiente:

  1. Dé un intervalo de confianza al 95% para el total de la variable \(REV84\) usando las fórmulas vistas en clase.

  2. Obtenga el intervalo de confianza usando el paquete survey en R.

Para este ejercicio Consideramos un diseño de muestreo conprobabilidad proporcional al tamaño SIN reemplazo (PPT sin reemplazo) pero una vez observa de la muestra podemos suponer que fue seleccionada CON reemplazo debido a que no tenemos manera de guardar las probabilidades de inclusion de segundo orden.

Y realizamos los calculos necesarios para estimar nuestros intervalos. En este caso consideramos la expresion (120) de las notas para calcular nuestras probabilidades de inclusión de primer orden \[\pi_k = min\{\frac{nx_k}{t_x}, 1\}\] Y la expresión conocida para calcular los pesos \[w_k = \frac{1}{\pi_k}\] y el calculo del factor de corrección por población finita (fpc) \[fpc = n/N\] con el objetivo de ayudar un poco a la función a calcular un mejor diseño y que sea un poco más semejante a un diseño SIN reemplazo.

De manera que calculando estos parámetros y agregandolos a nuestro dataframe obtenemos lo siguiente

N <- 284 # El tamaño de la poblacion es N = 284.
n <- 12 # El tamaño de la muestra es n = 12 denotados por la variable REV84.
#### Cargamos los datos
REV84 <- c(5246, 59877, 2208, 2546, 2903, 6850,
           3773, 4055, 4014, 38945, 1162, 4852)
P75 <- c(54, 671, 28, 27, 29, 62,
         42, 48, 33, 446, 12, 46)
datos <- data.frame(REV84, P75)

# Se recolectaron los siguientes datos de la variable REV84.
# P75 es la variable tamaño con la que se contruyen
# las probabilidades de inclusion. 
# El total de P75 en la poblacion es 8182.

###Calculo de las probabilidades de inclusión de primer orden:
datos$pik <- pmin(n * datos$P75 /8182 , 1) # utilizamos la expresion (120) de las notas
#sum(datos$pik) # NO es algo muy cercano a n = 12 (mala noticia)
#hist(datos$pik) # bastante dispersión en las probabilidades para solo tener 12 obs.

### Calculo de los pesos
datos$pesos <- 1/datos$pik # calculamos los pesos dado que ya tenemos las probas de inclusion

### Calculo del factor de correccion por poblacion finita
datos$fpc <- n/N  # de esta manera la estimacion puede asemejarse más a un
        # PPT SIN reemplazo aunque en realidad se sigue calculando CON reemplazo

head(datos) # Visualizamos los datos

Ahora que nuestro dataframe esta completo con los parámetros y valores necesarios podemos dedicarnos a aplicar las ecuaciones desarrolladas en clase.

En primer lugar estamos buscando estimar el total, por lo que decidimos utilizar 2 y compararlos. En primer lugar tenemos el estimar HT para totales muy conocido. En segundo lugar tenemos el estimador HT alternativo. Comparando los totales estimados concluimos que HT para totales dio mejor rendimiento, es por eso que decidimos continuar con este.

Las varianzas se calculan para este ultimo estimador en sus 2 opciones que vemos en las expresiones (59) y (60) de las notas.

Finalmente el paso importante fue elegir qué aproximación de intervalos nos convenía y nuestros datos cumplian con los requisitos. Notemos que las expresiones (82) y (83) nos brindan dos maneras de calcular estos intervalos, sin embargo, por tener una muestra pequeña para el tamaño de nuestra población en teoría el intervalo “correcto” estaria dado por la expresion (83): \[[\widehat{\theta}-t_{n-1, 1-\alpha/2}\sqrt{\widehat{V}(\widehat{\theta})}, \widehat{\theta}+t_{n-1, 1-\alpha/2}\sqrt{\widehat{V}(\widehat{\theta})}]\] Entonces veamos cómo quedaron nuestros calculos:

################################################################################
### Calculamos el intervalo y el total de REV84 usando las formulas de la clase.
################################################################################
### 
# Estimador HT para el total 
(TotalHTest <- sum(datos$REV84*datos$pesos)) # 787998.5
## [1] 787998.5
# Estimador ALTERNATIVO para el total 
TotalALTest <- N * (sum(datos$REV84*datos$pesos) / sum(datos$pesos)) # 1007947
# estimación de la varianza del total asumiendo muestreo con reemplazo (expresiones 59 y 60)
VarEst.a.TotalHTest <- (n/(n-1))*sum( (datos$REV84*datos$pesos - sum(datos$REV84*datos$pesos)/n)^2 )  #792781286
VarEst.b.TotalHTest <- (1-n/N)*(n/(n-1))*sum( (datos$REV84*datos$pesos - sum(datos$REV84*datos$pesos)/n)^2 )  #759283485


valor_critico_norm <- qnorm (0.975) # expresion 82 usa cuantil de una normal
valor_critico_t <- qt (0.975, df = n-1) # expreion 83 usa cuantil de una t-studen n-1 grados
(intervalo_norm <- c(TotalHTest - valor_critico_norm*sqrt(VarEst.b.TotalHTest),
                     TotalHTest + valor_critico_norm * sqrt(VarEst.b.TotalHTest) ) ) # 733991.5, 842005.5
## [1] 733991.5 842005.5
(intervalo_t <- c(TotalHTest - valor_critico_t * sqrt(VarEst.b.TotalHTest),
                  TotalHTest + valor_critico_t * sqrt(VarEst.b.TotalHTest) ) ) # 727350.1, 848646.9
## [1] 727350.1 848646.9

Asi concluimos que

  • Usando la expresion (82) nuestro intervalo es (733991.5, 842005.5)

  • Usando la expresion (83) nuestro intervalo es (727350.1, 848646.9)

Analogamente, tomando en cuenta lo explicado en el inciso anterior. Ahora nos interesa calcular estos intervalos con la paqueteria survey la cual nos ahorra tener que escribir las expresiones y simplemente pasarle nuestros parametros calculados al inicio del ejercicio.

Aquí cabe destacar que brindaremos el parámetro fpc a la hora de crear nuestro diseño para asemejarlo lo más que se pueda a un muestreo proporcional al tamaño sin reemplazo.

################################################################################
################# Calculamos el intervalo de confianza usando el paquete survey.
################################################################################
PPT <- svydesign(id=~1, weight=~pesos, fpc=~fpc, data=datos)
summary(PPT) # Independent Sampling design
## Independent Sampling design
## svydesign(id = ~1, weight = ~pesos, fpc = ~fpc, data = datos)
## Probabilities:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01760 0.04217 0.06453 0.18308 0.08213 0.98411 
## Population size (PSUs): 284 
## Data variables:
## [1] "REV84" "P75"   "pik"   "pesos" "fpc"
### Estimador HT para el total                  #        total     SE   DEff
(TotalHTest <- svytotal(~REV84, PPT, deff = T)) # REV84 787999  27555 0.0069
##        total     SE   DEff
## REV84 787999  27555 0.0069
#   Varianza del total
#(vcov(TotalHTest)) # 759283485

#   Intervalo de confianza para el total de REV84. Por defecto al 95%
                       #          2.5 %   97.5 %
(confint(TotalHTest) ) # REV84 733991.5 842005.5
##          2.5 %   97.5 %
## REV84 733991.5 842005.5
                                       #          2.5 %   97.5 %
(confint(TotalHTest, df = degf(PPT)) ) # REV84 727350.1 848646.9
##          2.5 %   97.5 %
## REV84 727350.1 848646.9

Observamos que afortunadamente los intervalos calculados con survey coinciden completamente con los intervalos calculados “manualmente” por lo que nos queda decir que en lo personal nos convence más el intervalo_norm (733991.5, 842005.5) porque es mas pequeño y contiene al total estimado 787998.5, pero la teoria dice que el correcto es el otro intervalo.

4. Comparación de diseños unietápicos

Suponga que cada mes se levanta un censo en N=502 puntos de venta para obtener los precios del jitomate y calcular el precio promedio. Además de este cálculo, lo importante de este seguimiento es analizar la diferencia de los precios promedio asociada a los dos últimos meses.

Dado que el movimiento de los precios ha sido muy variable, se plantea realizar este seguimiento con más oportunidad, de ser posible cada semana. Sin embargo, para hacer la comparación cada semana no existe mucho dinero y se plantea la opción de levantar una muestra del 10% de los puntos de venta. Para esto se analizan dos características para definir la selección de la muestra.

  1. Considerar un m.a.s. del 10% de los puntos de venta o un diseño Poisson con el mismo tamaño esperado de muestra, donde la variable tamaño se define por los precios de enero. Para el caso del m.a.s no se estimaría \(N\), pero para el diseño Poisson, dado que se usan probabilidades diferentes, se considera estimar \(N\).

  2. Considerar una sóla muestra que será visitada cada semana o bien seleccionar una muestra independiente cada semana.

Para definir la mejor estrategia se acercan a usted y la información que le proporcionan son los precios del jitomate en enero y febrero recolectada en los 502 puntos de venta (DatosJitomatePobEneFeb.csv).

Con esta información de apoyo, ayude a definir la mejor estrategia considerando como base la estimación de la diferencia de los precios promedio entre enero y febrero.

Para la solución de este ejercicio dividiremos nuestros análisis de 2 partes.

  • En la primera parte analizaremos el uso de un diseño m.a.s realizando un estudio longitudinal y un estudio transversal (muestras independientes en el tiempo).

  • En la segunda parte analizaremos el uso de un diseño Poisson realizando un estudio longitudinal y un estudio transversal al igual que en el anterior.

Finalmente en cada una de las fases elegiremos al mejor diseño-estudio y compararemos los 2 mejores para definir cuál es el diseño y estudio que recomendariamos seguir para un estudio semanal a posteriori.

Para comenzar nuestro desarrollo es importante realizar el siguiente análisis general.

Dado que nuestro objetivo es comparar la diferencia de los promedios de los precios en cada mes, entonces nuestro parámetro a estimar será: \[\theta = \bar{y}_1 - \bar{y}_2\] donde \(\bar{y}_1\) es el precio promedio en Enero y \(\bar{y}_2\) es el precio promedio en Febrero.

Notemos que si tomamos un estudio longitudinal, entonces nuestro estimador se puede ver como: \[ \widehat{\theta} = \frac{\widehat{t}_1-\widehat{t}_2}{\widehat{N_1}} = \frac{\sum_{k \in S} \frac{y_{1k}}{\pi_k}-\sum_{k \in S} \frac{y_{2k}}{\pi_k}}{\widehat{N_1}}= \frac{\sum_{k \in S} \frac{z_{k}}{\pi_k}}{\widehat{N_1}} \] donde \(z_k = y_{1k}-y_{2k}\) y el denominador es un total que se puede tratar de manera común.

Y su varianza dada por: \[ V(\widehat{\theta}) = V(\frac{\sum_{k \in S} \frac{z_{k}}{\pi_k}}{\widehat{N_1}}) \] > Nota: Para el diseño m.a.s. la \(N\) es conocida por lo que no haria falta estimarla por lo que la formula anterior no es exactamente la que se usa en el m.a.s. En el caso del diseño Poisson la \(N\) sí se debe de estimar y se ve como en la formula expresada anteriormente.

Ahora veamos que si el estudio es transversal entonces: \[ \widehat{\theta} = \frac{\widehat{t}_1}{\widehat{N_1}} -\frac{\widehat{t}_2}{\widehat{N_2}} \] En tal caso la varianza estará dada por \[ V(\widehat{\theta}) = V(\frac{\widehat{t}_1}{\widehat{N_1}} -\frac{\widehat{t}_2}{\widehat{N_2}}) =V(\frac{\widehat{t}_1}{\widehat{N_1}})+V(\frac{\widehat{t}_2}{\widehat{N_2}}) \] Esta propiedad se debe a que ambos estimadores son independientes por lo que, aplicando propiedades de la varianza obtenemos que \(Var(X-Y) = Var(X)+Var(Y)\) cuando \(X,Y\) son independientes.

Nota: aquí es lo mismo explicado en la nota anterior. La N es estima solo en el muestreo Poisson.

Ahora sí comencemos con el m.a.s.:

# Cargamos el archivo DatosJitomatePobEneFeb.csv
jitomate <- read.csv("C:/Users/noeam/Documents/git/github/Muestreo/data/DatosJitomatePobEneFeb.csv", header = TRUE, sep=",", fileEncoding="latin1")
N <- 502
n <- 50 # Tomamos nuestras muestras del 10% (aprox 50 datos)

# precio promedio "real" del jitomate en Enero
(meanJitomateEne2022 <- mean(jitomate$JitomateEne2022)) # 20.14016
## [1] 20.14016
# precio promedio "real" del jitomate en Febrero
(meanJitomateFeb2022 <- mean(jitomate$JitomateFeb2022)) # 16.28371
## [1] 16.28371
# diferencia de los precios promedio "real"
(meanDiferencia <- meanJitomateEne2022 - meanJitomateFeb2022 )# 3.856454
## [1] 3.856454
# Definimos zk = y2-y1, es decir la diferencia de los precios para las mismas tiendas
# entre febrero y enero
jitomate$zk <- jitomate$JitomateEne2022 - jitomate$JitomateFeb2022
(mean_zk <- mean(jitomate$zk)) # 3.856454 es lo mismo que la "real": meanDiferencia
## [1] 3.856454

Para nuestro diseño m.a.s debemos considerar 3 expresiones importantes:

  1. La probabilidades de inclusion de primer orden \[ \pi_k = \frac{n}{N} \]

  2. El estimador de la media HT para m.a.s. con \(N\) conocida expresion (99) de las notas \[ \widehat{Y_\pi} = \bar{y}_S \]

  3. La varianza poblacional (porque tenemos los datos poblacionales) para nuestro estimador expresión 100 \[ V(\widehat{Y_\pi}) = \frac{N^2}{n}(1-\frac{n}{N})S^2_{y_U}\frac{1}{N^2} \] > Nota: El estimador y la varianza ya estan expresados para las medias.

De esta manera, utilizando estas expresiones obtenemos lo siguiente

Para un estudio longitudinal M.A.S:

################################################################################
################################################################################
# Estudio LONGITUDINAL M.A.S
# Consideramos la diferencia de los precios promedio para dos muestras con los
# mismos indices en Enero y Febrero
#    Dado que se usan probabilidades diferentes, consideramos estimar N. Es decir
#    usaremos el estimador alternativo (109) para la media y la varianza de este
#    estimador (110) para cada una de las muestras
#    Como guia sabemos que N = 502
################################################################################
################################################################################
pk <- n/N  # Probabilidades de primer orden M.A.S.
pkl <- (n*(n-1))/(N*(N-1)) # Probabilidades de segundo orden M.A.S.

B <- 10000 #numero de muestras que se tomaran
set.seed(1234) 
muestra1 <- replicate(B, sample(1:N, n))
# Estimador HT mas con N conocida expresion (99)
MeanDifference.HT.L <- sapply(1:B,function(x, Indices, Pob){mean(Pob[Indices[,x],"zk"])}, Indices=muestra1, Pob=jitomate)
h1 <- hist(MeanDifference.HT.L, plot=FALSE, breaks=100)
plot(h1,
     freq=FALSE,
     main=TeX("$\\widehat{\\bar{Y}}_{\\pi}=\\widehat{t}_{\\pi}/N$"),
     xlab=paste("\n Estimación M.A.S Longitudinal \n E(ns)=",n, ", ECM=", round(mean((MeanDifference.HT.L-meanDiferencia)^2), 2) ,  sep=""),
     ylab="",
     cex.lab=.9,
     yaxt='n',
     xlim=c(0,8) )
abline(v=meanDiferencia, col="red")

#Expresion 100 varianza poblacional de un mas
(VarPob.MeanDifference.HT.L <- ((N^2/n)*(1-(n/N))*sum((jitomate$zk-meanDiferencia)^2)/N-1)/N^2)
## [1] 0.3027038

Para un estudio transversal M.A.S:

################################################################################
################################################################################
# Estudio TRANSVERSAL M.A.S
# Consideramos la diferencia de los precios promedio para dos muestras con 
# diferentes indices en Enero y Febrero
#    Como guia sabemos que N = 502
################################################################################
################################################################################
set.seed(1234) 
muestra2 <- replicate(B, sample(1:N, n))
muestra3 <- replicate(B, sample(1:N, n))
# Estimador HT mas con N conocida expresion (99)
MeanDifference.HT.T <- sapply(1:B,function(x, Indices1, Indices2, Pob){mean(Pob[Indices1[,x],"JitomateEne2022"]) - mean(Pob[Indices2[,x],"JitomateFeb2022"])},
                            Indices1=muestra2,
                            Indices2=muestra3,
                            Pob = jitomate)
h2 <- hist(MeanDifference.HT.T, plot=FALSE, breaks=100)
plot(h2,
     freq=FALSE,
     main=TeX("$\\widehat{\\bar{Y}}_{\\pi}=\\widehat{t}_{\\pi}/N$"),
     xlab=paste("\n Estimación M.A.S Transversal\n E(ns)=",n, ", ECM=", round(mean((MeanDifference.HT.T-meanDiferencia)^2), 2) ,  sep=""),
     ylab="",
     cex.lab=.9,
     yaxt='n',
     xlim=c(0,8) )
abline(v=meanDiferencia, col="red")

#Expresion 100 varianza poblacional de un mas
VarPob.MeanEne.HT <- ((N^2/n)*(1-(n/N))*sum((jitomate$JitomateEne2022-meanJitomateEne2022)^2)/N-1)/N^2
VarPob.MeanFeb.HT <- ((N^2/n)*(1-(n/N))*sum((jitomate$JitomateFeb2022-meanJitomateFeb2022)^2)/N-1)/N^2

# Var(X-Y) = Var(X) + Var(Y) cuando X, Y son independientes
(VarPob.MeanDifference.HT.T <- VarPob.MeanEne.HT + VarPob.MeanFeb.HT)
## [1] 1.948116

Y de aquí concluimos que la estimación usando un estudio longitudinal es mucho mejor debido a que la dispersión, varianza y ECM son menores, lo cual se refleja claramente en la gráfica que nos muestra una concentrácion de los datos más cercanda a la media.

Por lo tanto si usamos un diseño m.a.s lo conveniente es realizar un estudio longitudinal.

Para el diseño Poisson también hay 3 expresiones a considerar:

  1. La probabilidades de inclusion de primer orden para un diseño Poisson que cumple con tener tamaños \(x_k\) positivos. Expresion (“108.5”) de las notas \[ \pi_k = \min\{\frac{nx_k}{t_x}, 1 \} \]

  2. El estimador alternativo para tamaños de muestra variables (diseño Poisson) con probabilidades de inclusion diferentes y \(N\) a estimar \[ \widehat{Y_\pi}_{alt} = \frac{\widehat{t}_{y\pi}}{\widehat{N}_\pi} \]

  3. La varianza poblacional (porque tenemos los datos poblacionales) para nuestro estimador expresión 100 \[ V(\widehat{Y_\pi}_{alt}) = \frac{1}{N^2}\sum_{k\in U}\frac{(y_k - \bar{y}_U)^2}{\pi_k}(1-\pi_k) \] > Nota: Nuevamente estas expresiones ya estan calculadas para estimar medias y calcular la varianza de estimadores de medias.

Así que con esto en mente desarrollamos nuestro problema utilizando estas expresiones encontradas:

################################################################################
################################################################################
# Ahoras consideremos un diseño Poisson con un tamaño de muestra n = 50
################################################################################
################################################################################

# Primero calculamos las probabilidades de primer orden usando expresion (108.5)
# Esto debido a que nuestra variable tamaño X_k es positiva siempre y 
# cumple los supuestos requeridos para este calculo.
jitomate$pik <- pmin(n*jitomate$JitomateEne2022/sum(jitomate$JitomateEne2022), 1)
describe(jitomate$pik) # descripcion de las pi_k's
min(jitomate$pik) # el valor minimo de las pi_k's
## [1] 0.04638806
sum(jitomate$pik) # cumple con la teoria, es decir suma n = 50.
## [1] 50
# Utilizamos el muestreo Poisson utilizado en clase
samplePoisson <- function(pks){
  Nd <- length(pks)
  index <- 1:Nd
  unif <- runif(Nd)
  Unos <- (unif < pks) * 1
  indexs <- index[Unos==1]
  return(indexs)
}

Para un estudio longitudinal Poisson:

################################################################################
################################################################################
# Estudio LONGITUDINAL Poisson
# Consideramos la diferencia de los precios promedio para dos muestras con los
# mismos indices en Enero y Febrero
#    Dado que se usan probabilidades diferentes, consideramos estimar N. Es decir
#    usaremos el estimador alternativo (109) para la media y la varianza de este
#    estimador (110) para cada una de las muestras
#    Como guia sabemos que N = 502
#    Definimos z_k = y1_k - y2_k (febrero menos enero) debido al algebrita
#    Mostrada anteriormente
################################################################################
################################################################################
# Tomamos B muestras con diseño Poisson
set.seed(1234)
muestra4 <- replicate(B, samplePoisson(jitomate$pik)) # indices de las B muestras

# Ahora calculamos utilizando estimador alternativo ecuacion 109 es mejor cuando las probas de inclusion son totalmente diferentes
# Estimador de la diferencia de medias con diseño Poisson y estudio Longitudinal
estDiferenciaPL <- sapply(1:B, function(x, Indices, Pob){sum(Pob[Indices[[x]],"zk"]/Pob[Indices[[x]],"pik"])/sum(1/Pob[Indices[[x]],"pik"]) },
                                    Indices=muestra4,
                                    Pob=jitomate)
h3 <- hist(estDiferenciaPL, plot=FALSE, breaks=100)
plot(h3,
     freq=FALSE,
     main=TeX("$\\widehat{\\bar{Y}}_{alt}=\\widehat{t}_{\\pi}/\\widehat{N}_{\\pi}$"),
     xlab=paste("\n Estimación Poisson Longitudinal\n E(ns)=",n, ", ECM=", round(mean((estDiferenciaPL-meanDiferencia)^2), 2) ,  sep=""),
     ylab="",
     cex.lab=.9,
     yaxt='n',
     xlim=c(0,8) )
abline(v=meanDiferencia, col="red")

# Expresion 110 varianza poblacional del estimador alternativo para Poisson
(VarPob.estDiferenciaPL <- sum(((jitomate$zk-meanDiferencia)^2)*(1-jitomate$pik)/jitomate$pik)/N^2 )
## [1] 0.2313726

Para un estudio transversal Poisson:

################################################################################
################################################################################
# Estudio TRANSVERSAL Poisson
# Consideramos la diferencia de los precios promedio para dos muestras con 
# diferentes indices en Enero y Febrero
#    Dado que se usan probabilidades diferentes, consideramos estimar N.
#    Como guia sabemos que N = 502
################################################################################
################################################################################
# Ahora calculamos utilizando estimador alternativo ecuacion 109 es mejor cuando las probas de inclusion son totalmente diferentes
set.seed(1234) 
muestra5 <- replicate(B, samplePoisson(jitomate$pik))
muestra6 <- replicate(B, samplePoisson(jitomate$pik))
# Ahora calculamos utilizando estimador alternativo ecuacion 109 es mejor cuando las probas de inclusion son totalmente diferentes
# Para Enero
MeanJitomateEne2022.alt.P <- sapply(1:B, function(x, Indices, Pob){sum(Pob[Indices[[x]],"JitomateEne2022"]/Pob[Indices[[x]],"pik"])/sum(1/Pob[Indices[[x]],"pik"]) },
                               Indices=muestra5,
                               Pob=jitomate)

# Ahora sacamos el estimador de la media para Febrero con indices diferentes pues es el caso de muestras independientes
MeanJitomateFeb2022.alt.P <- sapply(1:B, function(x, Indices, Pob){sum(Pob[Indices[[x]],"JitomateFeb2022"]/Pob[Indices[[x]],"pik"])/sum(1/Pob[Indices[[x]],"pik"]) },
                               Indices=muestra6,
                               Pob=jitomate)

# Estimador de la diferencia de medias con diseño Poisson y estudio Transversal
estDiferenciaPT <- MeanJitomateEne2022.alt.P - MeanJitomateFeb2022.alt.P
h4 <- hist(estDiferenciaPT, plot=FALSE, breaks=100)
plot(h4,
     freq=FALSE,
     main=TeX("$\\widehat{\\bar{Y}}_{alt}=\\widehat{t}_{\\pi}/\\widehat{N}_{\\pi}$"),
     xlab=paste("\n Estimación Poisson Transversal\n E(ns)=",n, ", ECM=", round(mean((estDiferenciaPT-meanDiferencia)^2), 2) ,  sep=""),
     ylab="",
     cex.lab=.9,
     yaxt='n',
     xlim=c(0,8) )
abline(v=meanDiferencia, col="red")

# Expresion 110 varianza poblacional del estimador alternativo para Poisson
VarPob.MeanJitomateEne2022.alt.P <- sum(((jitomate$JitomateEne2022-meanJitomateEne2022)^2)*(1-jitomate$pik)/jitomate$pik)/N^2 
VarPob.MeanJitomateFeb2022.alt.P <- sum(((jitomate$JitomateFeb2022-meanJitomateFeb2022)^2)*(1-jitomate$pik)/jitomate$pik)/N^2

# Var(X-Y) = Var(X) + Var(Y) cuando X, Y son independientes
(VarPob.estDiferenciaPT <- VarPob.MeanJitomateEne2022.alt.P + VarPob.MeanJitomateFeb2022.alt.P)
## [1] 1.505065

CONCLUSIÓN

Finalmente, después de observar los resultados y las gráficas, podemos concluir que las estimaciones se dan mejor cuando realizamos estudios longitudinales pues estas simulaciones suelen tener una menor dispersión y menor ECM, por lo que de ser posible lo mejor será realizar un estudio longitudinal, es decir, recabar datos de las mismas tiendas cada semana e ir dandoles seguimiento para tener la ventaja de analizar el comportamiento de los precios en cada una de las tiendas seleccionadas desde el inicio del estudio.

Por otro lado, ya que notamos que el estudio longitudinal arroja mejores resultados queda elegir entre el diseño M.A.S y el Poisson.

En este caso, si nos guiamos por los ECM, podemos notar que el diseño Poisson tiene un ECM ligeramente menor al M.A.S por lo que podría darnos aún mejores resultados. Pero si tomamos en cuenta que el M.A.S necesita que conozcamos los datos poblacionales entonces el Poisson vuelve a ganar, debido a que en la práctica es difícil obtener datos poblacionales por lo que estimar estos parámetros suele ser algo más común.

Así que, dado el análisis descrito, nosotros consideramos que la mejor estrategia a seguir sería implementar un diseño de muestreo Poisson con, idealmente, un estudio longitudinal para dar el seguimiento semanal que se plantea y obtener los mejores resultados posibles.