Utilización de medidas parciales de la variable de salida en la modelización

Question

Mi pregunta es: ¿Cómo podemos utilizar datos de salida parcialmente medidos en un conjunto de entrenamiento? Esto es vago, así que lo concreto en un cuento caprichoso.

Las ardillas tienen nueces, pero ¿cuántas?

Configurar

Hay un conjunto $S$ de ardillas y un conjunto $T$ de los árboles del bosque. Ardilla $s_i \in S$ tiene $n_i \in N$ nueces. Las ardillas almacenan sus nueces en un árbol o en muchos árboles diferentes del bosque. Queremos predecir cuántas nueces ha recogido una ardilla determinada a partir de las características de entrada: el peso de la ardilla y la capacidad de sus mejillas.

Fase de investigación 1:

Hemos monitorizado un subconjunto de ardillas. Tomamos sus medidas (peso y capacidad de las mejillas), y contamos

1. cuántas nueces recogieron y
2. en cuántos árboles los almacenaron

(por ejemplo, la ardilla #55 almacenó 5 nueces en un árbol, 10 nueces en otro y 500 nueces en otro [se desconoce en qué 3 árboles])

Fase de investigación 2:

Controlamos un subconjunto de árboles. Cuando una ardilla llegaba a uno de nuestros árboles, la identificábamos (para saber si iba a otro árbol de nuestro subconjunto de estudio más tarde), la medíamos (peso y capacidad de las mejillas) y contábamos cuántas nueces soltaba. Esto nos dio una parcial medición de la recogida de frutos secos para un subconjunto de ardillas.

(por ejemplo, en el árbol nº 23 recogimos 10 nueces de la ardilla nº 99 y 50 de la ardilla nº 88, en el árbol nº 24 recogimos...)

(Nota importante: las identificaciones de las ardillas no persisten a través de las fases).

Pregunta

Supongamos que deseamos modelar el número de nueces totales de una ardilla de la fase 2. ¿Cómo podríamos utilizar sus mediciones parciales de nueces para aumentar los resultados de la modelización?

Además, ¿cómo podemos introducir resultados parcialmente medidos en el conjunto de entrenamiento?

Answer 1

Pregunta divertida. El problema clave, como señala @MartijnWeterings, es que el número de árboles en la fase 2 es sólo una medida parcial del número total de árboles. Sin embargo, si conociéramos el número total de árboles, podríamos resolver el problema construyendo un modelo del número de nueces observado en la fase 1 dado el número de árboles en la fase 1, y luego predecir el número de nueces en la fase 2 utilizando el número de árboles en la fase 2. Nuestra estrategia en esta respuesta es, por tanto, estimar primero el número de árboles en la etapa 2 y luego construir un modelo de frutos secos dados los árboles en la etapa 1.

Notación y suposición

A continuación, asumo que el muestreo de árboles y ardillas es aleatorio en todas las etapas. Sea $n_{1i}$ denotan la suma de todas las nueces recogidas por la ardilla $i$ en la fase 1. Sea $t_{1i}$ denotan el número total de árboles ardilla $i$ tuercas almacenadas en la fase 1. Dejemos que $n_{2j}$ denotan la suma no observada de nueces recogidas por la ardilla $j$ en la fase 2 y que $t_{2j}$ denotan el número de árboles ardilla $j$ tuercas almacenadas en la fase 2. Por último, dejemos $k_{2j}$ denotan el número parcial de árboles observados, donde $k_{2j} \le t_{2j}$ ,

Número de árboles en la fase 2

Como señala @MartijnWeterings $k_{2j}$ es siempre menor o igual que el número total de árboles $t_{2j}$ en la fase 2, que se desconoce. Nuestro objetivo es, por tanto, estimar $t_{2j}$ lo más cerca posible. Afortunadamente, tenemos alguna información sobre $t_{2j}$ . Dependiendo de su diseño de muestreo en la fase 2, hay una probabilidad $p$ que una ardilla sea capturada en una de las $t_{2j}$ árboles que visita. La probabilidad de $k_{2j}$ es por tanto binomial con parámetros $t_{2j}$ y $p$ . Sin embargo, nosotros observe binomio $k_{2j}$ el número de árboles $t_{2j}$ sin embargo, no tiene una distribución binomial dado $k_{2j}$ . No estaba seguro de su distribución exacta y por eso hice una pregunta al respecto en Matemáticas-Intercambio de pila . Recibí la útil respuesta de que la función de masa de probabilidad de $t=t_{2j}$ con $k=k_{2j}$ y $p$ viene dada por $$P(t; k ,p) = \binom{t-1}{k} p^t (1-p)^{(t-k)}, \quad t \in \{k,...,\infty\}.$$ para todos $j$ que tiene la expectativa $E(t)=k/p$ . Por lo tanto, si sabemos $k_{2j}$ y $p$ podríamos estimar $\hat{t}_{2j}=k_{2j}/p$ . Como se ha dicho, la probabilidad $p$ depende de su diseño de muestreo, pero afortunadamente podemos estimarlo a partir de los datos como $$\hat{p}=\frac{\sum_{j} k_{2j}}{\sum_{i} t_{1i}}$$ para que $\hat{t}_{2j}=k_{2j}/\hat{p}$ .

Estimación bajo el supuesto de proporcionalidad

Dejemos que

$$ \pi = \frac{1}{\#S_1} \sum_{i} \frac{n_{1i}}{t_{1i}}$$

sea la proporción media de nueces que deja una ardilla en un árbol. Una primera estimación del número total de nueces de la ardilla $j$ es

$$ \hat{n}_{2j} = \pi \hat{t}_{2j}.$$

Estimación mediante la relación entre los frutos secos y los árboles en la fase 1

Esto puede parecer insatisfactorio, porque no tiene en cuenta que puede haber una relación entre $n$ y $t$ que no sea una simple proporcionalidad. Por ejemplo, podemos imaginar que las ardillas tienen el extraño comportamiento de dejar menos nueces por árbol cuanto más nueces tienen a su disposición. Entonces el número total de nueces $n$ no aumentaría proporcionalmente con $t$ y en su lugar se aplana. Por lo tanto, podríamos decidir modelar

$$ n_{1i}= f(t_{1i},\theta) + \epsilon_i$$

donde $f$ es una función no lineal con parámetros theta y $\epsilon_i$ es un término de error de medición. Una opción obvia podría ser

$$ n_{1i} = \theta_0 + \theta_1 \log(t_{1i}) + \epsilon_i$$

con $\epsilon_i$ iid normal con expectativa 0. El modelo puede ajustarse por mínimos cuadrados no lineales o por máxima verosimilitud. Un estimador sería entonces

$$ \hat{n}_{2j} = \hat{\theta_0} + \hat{\theta_1} \log(\hat{t}_{2j})$$

Por supuesto, se podrían utilizar otras formas funcionales o se podrían utilizar técnicas de modelización flexibles para aproximar la relación funcional, como por ejemplo bosques aleatorios (juego de palabras).

Simulaciones

¿Funciona esto? Vamos a probarlo. Simulo los datos en R según las siguientes ideas. La probabilidad de que una ardilla recoja $n+1$ tuercas viene dado por $n \sim \text{Poisson}(20)$ . Una ardilla llega entonces al primer árbol y se esconde $h_1+1$ tuercas donde $h_1 \sim \text{Poisson}(\lambda)$ y $\lambda \sim \Gamma(60/n,1)$ . Sigue escondido en $1 + (h_2,...,h_t)$ nueces hasta que llegue al árbol $t$ y no tiene nueces. Lo hace independientemente de si lo observa en la fase 1 o en la 2; sin embargo, en la fase 1 observa todo $h_t$ mientras que en la fase 2 se observa una muestra de $\{h_1,...,h_t\}$ . Como se ha dicho asumo que se tiene una muestra aleatoria simple de árboles en la fase 2 y así se observa $h_{kj}$ (el k-ésimo árbol visitado por la ardilla j) con probabilidad $p$ (abajo en el código lo llamo truncamiento).

Ahora muestro 1000 ardillas en la fase 1. El gráfico siguiente ilustra la relación entre el número total de árboles y el número total de nueces recogidas. Se puede ver que hay una decadencia en esa relación a través de $t$ .

Ahora muestro en la etapa 2 con $p=0.5$ y considerar tres estimadores. Primero el estimador bajo proporcionalidad. En segundo lugar, creo un estimador que utiliza la media condicional de $n_1$ en cada nivel observado de $t_1$ como una estimación para $n_2$ en $\hat{t}_2$ . Para la evaluación comparativa utilizo de nuevo la media condicional de $n_1$ en cada nivel observado de $t_1$ como una estimación para $n_2$ pero ahora con el verdadero número de árboles $t_2$ en la fase 2. Por supuesto, este estimador no está disponible en la práctica.

Para dos muestras, una de cada fase 1 y 2, respectivamente, y los tres estimadores llego a los siguientes sesgos, respectivamente 5,61, -0,19 y 0,24. Además, observamos los siguientes errores cuadráticos medios 15.3, 4.07, 3.32. Vemos que el estimador de la media condicional con una estimación ajustada del número de árboles en la fase 2 tiene un rendimiento casi tan bueno como el estimador que utiliza el número verdadero desconocido de árboles en la fase 2. El error restante es la varianza, que quizá pueda reducirse un poco más utilizando un modelo mejor para $n_1$ dado $t_1$ que el modelo de media condicional no paramétrico.

---

Aquí hay una función que crea los datos para la simulación que hice.

# A squirrel collects nuts
squirrel_set = function(n, trunc= FALSE){
  nuts = rpois(n, 20) + 1
  nut_seq = list()
  for(i in 1:n){
    j = 1
    nuts_left = nuts[i]
    nuts_hidden = numeric()
    # squirrel hides nuts at tree j
    while(nuts_left>0){
      if(j  == 1) {lambda = rgamma(1,60/nuts_left,1) } 
      if(lambda < 1){ lambda = 1}
      nuts_hidden[j]  = rpois(1, lambda) + 1
      if(nuts_left - nuts_hidden[j] <0){
        nuts_hidden[j] = nuts_left
        nuts_left = 0
      }
      else{ nuts_left =  nuts_left - nuts_hidden[j] }
      j = j+1
    }
    nut_seq[[i]] = nuts_hidden
  } 
  # Truncated sample
  # A squirrel is caught with probability .5 at a tree 
  # (or a random half of the trees are observed and a squirrel is always caught)
  if(trunc == TRUE){
    trees = sapply(nut_seq , length)
    nut_seq_obs = list()
    for(i in 1:length(nut_seq)){
      caught = rbinom(trees[i],1,.5)
      nut_seq_obs[[i]] = nut_seq[[i]][as.logical(caught)]
    }
    return( list(nut_seq,nut_seq_obs) )
  }
  else{
    return(nut_seq)
  }
}

Y aquí el código utilizado en el análisis:

set.seed(12345)
n = 1000
# Phase 1
nut_seq1 = squirrel_set(n)

# Phase 2
nut_seq2 = squirrel_set(n, trunc = TRUE)
nut_seq2_true = nut_seq2[[1]]
nut_seq2_trunc = nut_seq2[[2]]

# Trees and nuts at phases 1 and 2
t1  = sapply(nut_seq1,length, simplify = TRUE) # Trees seen at phase 1

k   = sapply(nut_seq2_trunc , length) # Trees seen at phase 2
nut_seq2_trunc = nut_seq2_trunc[k>0] # Squirrels with k=0 have avtually not been observed
nut_seq2_true = nut_seq2_true[k>0]
k   = k[k>0]
n1  = sapply(nut_seq1,sum, simplify = TRUE) # Trees seen at phase 1
n2  = sapply(nut_seq2_true,sum, simplify = TRUE) # Trees at phase 2 (unobserved)
t2  = sapply(nut_seq2_true,length, simplify = TRUE) # Trees at phase 2 (unobserved)

# Plot
plot( t1, n1, xlab='Trees at phase 1', ylab='Total number of nuts at phase 1')
mnuts = numeric()
for(i in 1:max(t1)){
  mnuts[i] = mean(n1[t1 == i])
}
lines( 1:max(t1), mnuts, col=2)
legend("bottomright",lty=1, lwd=2, col='2', legend='Conditional mean')

# Estimators
p           = sum(k) / sum(t1) # Estimate of observational probability at phase 2
t2_est      = k/p  # Estimate of total number of trees for each squirrel at phase 2

n2_prop_est = t2_est * mean(sapply(n1,sum, simplify = TRUE)/t1 )  # proportionality

mnuts = numeric()
for(i in 1:max(t1)){
  mnuts[i] = mean(n1[t1 == i]) # Conditional mean at each level of trees at phase 1
}
n2_regest = mnuts[round(t2_est)] # Non-parametric regression estimate of n2
n2_regest_truet2 = mnuts[t2] 

# Bias and Variance
mean( n2_prop_est - n2)
sqrt(mean( (n2_prop_est - n2)^2 ))

mean( n2_regest - n2 , na.rm=TRUE)
sqrt(mean( (n2_regest - n2)^2 , na.rm=TRUE))

mean( n2_regest_truet2 - n2 , na.rm=TRUE)
sqrt(mean( (n2_regest_truet2 - n2)^2 , na.rm=TRUE))

Answer 2

En la fase 1 podrías hacer un modelo que relacione el comportamiento de una ardilla con el número total de nueces.

En la fase 2 no se observa exactamente la misma información completa que en la fase 1. Por ejemplo, no se sabe cuántos árboles en total utiliza una ardilla concreta. Pero sí se observa parte del comportamiento de la ardilla, a saber muestra de la distribución del número de nueces por árbol. A partir de esto se puede estimar la distribución y los parámetros que describen la distribución pueden ser introducidos para el modelo.

• Así que en la fase 1 se hace un modelo que relaciona el número total de nueces que almacena una ardilla con la distribución de las nueces por árbol para la ardilla. Cómo modelar esto exactamente es difícil de decir.

  Si crees que puedes hacer algún modelo mecanicista entonces podrías empezar con algún análisis exploratorio y conocimientos previos sobre el comportamiento de las ardillas para tener una idea sobre un modelo útil. Me faltan los datos y el conocimiento biológico para hacer esto en esta respuesta (Una dirección obvia podría ser ver si más nueces por árbol también se relacionan con más nueces en total, y posiblemente esto tendrá alguna relación más compleja con el peso y la mejilla de la ardilla, y otros factores, como la alta variación en las nueces por árbol puede ayudar a obtener también una indicación sobre el número total de árboles utilizados por una ardilla)

• En la fase 2 se hará una estimación de los parámetros que se necesitan para hacer predicciones con el modelo que se ha creado en la fase 1. Los parámetros que describen la distribución para el número de nueces por árbol pueden estimarse a partir de la muestra medida en el subconjunto de árboles.

  Una forma sencilla sería ignorar las identificaciones de los árboles y sólo utilizar los datos por ardilla para estimar los parámetros de distribución y ponerlos en el modelo de la fase 1.

  Un modelo más preciso modelo trataría las identificaciones de los árboles como un factor aleatorio, de manera que el comportamiento específicamente atribuido a la ardilla puede ser mejor estimado. Para tratar los árboles como un factor aleatorio habrá que saber cómo los árboles pueden ser un factor aleatorio. Se puede hacer una conjetura para esto, pero también podrías intentar aprenderlo de los datos (yo diría que con algún análisis exploratorio primero, comprobando la correlación entre un árbol y cuántas nueces se almacenan en él por ardilla, y si este efecto es independiente o quizás algunos árboles atraen a un tipo específico de ardillas, antes de llegar a algo cuantitativo). En la fase 1 no se observa información relacionada con las identificaciones de los árboles, pero en la fase 2 lo haces y puedes usar esos datos.

Así que en pocas palabras. Creo que se necesita un análisis exploratorio antes de poder hacer algo cuantitativo que vaya más allá del enfoque simple (simple es ignorar las identificaciones de los árboles en la fase 2 y usar sólo parámetros de distribución simples como entrada para el modelo cuyos coeficientes se aprenden en la fase 1).

---

Además, ¿cómo podemos introducir resultados parcialmente medidos en el conjunto de entrenamiento?

Cuando se hace el modelo en la fase 1 utilizando parámetros que describen la distribución de las nueces por árbol para una ardilla específica hay que tener en cuenta que debe ser posible estimar razonablemente esos números en la fase 2 y que los errores no afectarán demasiado al modelo. Por ejemplo, la media y la varianza (u otros estadísticos simples) podrían estimarse razonablemente a partir de las muestras de la fase 2 (suponiendo que la muestra no sea demasiado pequeña), pero los momentos de orden superior no.