# Using partial measurements of output variable in modeling

My question is: How can we use partially measured output data in a training set? This is vague, so I concretize it in a whimsical tale.

# Squirrels Have Nuts, But How Many?

## Setup

There is a set $$S$$ of squirrels and a set $$T$$ of trees in the forest. Squirrel $$s_i \in S$$ has $$n_i \in N$$ nuts. Squirrels store their nuts in one tree or many different trees in the forest. We wish to predict how many nuts a given squirrel has collected from input features: squirrel weight, and cheek capacity.

### Research phase 1:

We monitored a subset of squirrels. Took their measurements (weight and cheek capacity), and counted

1. how many nuts they collected and
2. how many trees they stored them in

(e.g. Squirrel #55 stored 5 nuts in one tree, 10 nuts in another, and 500 nuts in another [which 3 trees is unknown])

### Research phase 2:

We monitored a subset of trees. When a squirrel arrived at one of our trees, we I.D.'ed them (to know if they went to another tree in our study subset later), measured them (weight and cheek capacity), and counted how many nuts they dropped off. This gave us a partial nut collection measurement for a subset of squirrels.

(e.g. in Tree #23 we collected 10 nuts from Squirrel #99 and 50 from Squirrel #88, in Tree #24 we collected...)

(important note: Squirrel IDs don't persist across phases.)

## Question

Suppose we wish to model the number of total nuts of a squirrel from phase 2. How could we use their partial nut measurements to augment the modelling results?

Furthermore, how can we introduce partially measured outputs into the training set?

• Do you know the number of trees they were captured at at phase 2? – tomka Oct 1 '18 at 14:10
• And at phase 1: do you only know the total count of nuts they collected or do you know this per tree? (E.g. Squirrel 1 has collected 5 nuts from tree 1, 3 nuts from tree 2, and 10 nutes from tree 3; squrrel 2 has collected 7 nuts from (a different) tree 1, etc.. – tomka Oct 1 '18 at 14:19
• What fraction of the total is the subset of trees from phase 2?(if this fraction is large then you could reasonably estimate the total number of nuts collected by a squirrel based on the subset of nuts divided by that fraction) – Sextus Empiricus Oct 1 '18 at 22:46
• Please be more specific, do you know how many nuts they stored in each tree at phase 1 or 2 or both? – tomka Oct 2 '18 at 20:35
• Interesting question. Coincidentally, at the very moment I clicked this link, a squirrel came to eat some nuts scattered for him in my yard! Phase 1 is underway... – user20160 Oct 7 '18 at 16:21

Fun question. The key problem as noted by @MartijnWeterings is that the number of trees at phase 2 is only a partial measurement of the total number of trees. If we knew the total number of trees, however, we could solve the problem by building a model of the number of nuts observed at stage 1 given the number of trees at stage 1, and then predict the number of nuts at stage 2 using the number of trees at stage 2. Our strategy in this answer is therefore to first estimate the number of trees at stage 2 and then build a model of nuts given trees at stage 1.

## Notation and assumption

In the following, I assume that the sampling of trees and squirrels is random at all stages. Let $$n_{1i}$$ denote the sum of all nuts collected by squirrel $$i$$ in phase 1. Let $$t_{1i}$$ denote the total number of trees squirrel $$i$$ stored nuts at in phase 1. Let $$n_{2j}$$ denote the unobserved sum of nuts collected by squirrel $$j$$ in phase 2 and let $$t_{2j}$$ denote the number of trees squirrel $$j$$ stored nuts at in phase 2. Finally let $$k_{2j}$$ denote the partial number of trees observed, where $$k_{2j} \le t_{2j}$$,

## Number of trees at stage 2

As noted by @MartijnWeterings $$k_{2j}$$ is always smaller or equal to the total number of trees $$t_{2j}$$ at phase 2, which is unknown. Our goal thus becomes that of estimating $$t_{2j}$$ as closely as possible. Fortunately, we have some information on $$t_{2j}$$. Depending on your sampling design in phase 2, there is a probability $$p$$ that a squirrel is captured at one of the total $$t_{2j}$$ trees that it visits. The probability of $$k_{2j}$$ is thus binomial with parameters $$t_{2j}$$ and $$p$$. However, we observe binomial $$k_{2j}$$; the number of trees $$t_{2j}$$, however, is not binomial distributed given $$k_{2j}$$. I was not sure about its exact distribution and therefore I asked a question about it on Mathematics-StackExchange. I received the useful reply that the probability mass function of $$t=t_{2j}$$ with $$k=k_{2j}$$ and $$p$$ is given by $$P(t; k ,p) = \binom{t-1}{k} p^t (1-p)^{(t-k)}, \quad t \in \{k,...,\infty\}.$$ for all $$j$$ which has expectation $$E(t)=k/p$$. Hence if we know $$k_{2j}$$ and $$p$$ we could estimate $$\hat{t}_{2j}=k_{2j}/p$$. As said, probability $$p$$ depends on your sampling design, but fortunately we can estimate it from the data as $$\hat{p}=\frac{\sum_{j} k_{2j}}{\sum_{i} t_{1i}}$$ so that $$\hat{t}_{2j}=k_{2j}/\hat{p}$$.

## Estimation under proportionality assumption

Let

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

be the average proportion of nuts left by a squirrel at a tree. A first estimate of the total number of nuts of squirrel $$j$$ is

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

## Estimation using relationship between nuts and trees at phase 1

This may seem unsatisfactory, because it does not take into account that there may be a relationship between $$n$$ and $$t$$ other than a simple proportional one. For example we may imagine squirrels having the strange behavior of leaving less nuts per tree the more nuts they have at their disposal. Then the total number of nuts $$n$$ would not proportionately increase with $$t$$ and instead flatten off. Hence we could decide to model

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

where $$f$$ is a non-linear function with parameters theta and $$\epsilon_i$$ is a measurement error term. An obvious choice might be

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

with $$\epsilon_i$$ iid normal with 0 expectation. The model could be fit by non-linear least squares or maximum likelihood. An estimator would then be

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

Of course other functional forms could be used or you could use flexible modeling techniques to approximate the functional relationship, such as random forests (pun intended).

## Simulations

Does this work? Let's try it. I simulate data in R according to the following ideas. The probability that a squirrel collects $$n+1$$ nuts is given by $$n \sim \text{Poisson}(20)$$. A squirrel then arrives at the first tree and hides $$h_1+1$$ nuts where $$h_1 \sim \text{Poisson}(\lambda)$$ and $$\lambda \sim \Gamma(60/n,1)$$. It continues hiding at $$1 + (h_2,...,h_t)$$ nuts until it arrives at tree $$t$$ and is out of nuts. It does so regardless of whether you observe it in phase 1 or 2; however in phase 1 you observe all $$h_t$$, whereas in phase 2 you observe a sample from $$\{h_1,...,h_t\}$$. As said I assume you have a simple random sample of trees at phase 2 and so you observe $$h_{kj}$$ (the k-th tree visited by squirrel j) with probability $$p$$ (below in the code I call this truncation).

I now sample 1000 squirrels at phase 1. The plot below illustrates the relationship of the total number of trees and total number of nuts collected. It can be seen that there is a decay in that relationship across $$t$$. I now sample at stage 2 with $$p=0.5$$ and consider three estimators. First the estimator under proportionality. Second, I create an estimator which uses the conditional mean of $$n_1$$ at each observed level of $$t_1$$ as an estimate for $$n_2$$ at $$\hat{t}_2$$. For benchmarking I use again the conditional mean of $$n_1$$ at each observed level of $$t_1$$ as an estimate for $$n_2$$, but now at the true number of trees $$t_2$$ at phase 2. This estimator is of course not available in practice.

For two samples, one from each of phase 1 and 2, respectively, and the three estimators I arrive at the following biases, respectively: 5.61, -0.19, and 0.24. Furthermore we observe the following root mean square errors: 15.3, 4.07, 3.32. We see that the conditional mean estimator with an adjusted estimate for the number of trees at phase 2 has almost as good performance as the estimator using the unknown true number of trees at phase 2. The remaining error is variance which can perhaps be reduced a bit further by using a better model for $$n_1$$ given $$t_1$$ than the non-parametric conditional mean model.

Here is a function creating the data for the simulation I made.

# 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)
}
}


And here the code used in analysis:

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[]
nut_seq2_trunc = nut_seq2[]

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

• This is a very nicely written answer and using phase 1 to get some relationship between $n_{1i}$ and $t_{1i}$ seems the way to go. But, the problem is not cracked that easily. The $t_{2i}$ in phase 2 is a different thing. It is based on measurements from a subset. If one is lucky then there is a clear relationship between the number of total nuts and the number of nuts per tree. Since this 'number of nuts per tree' can be reasonably estimated in phase 2. – Sextus Empiricus Oct 1 '18 at 22:36
• Your thoughtfulness and puns are appreciated, and I also like the idea of finding relationships between number of nuts and number of trees. @MartijnWeterings is correct though. Squirrel j stored nuts at at least t2j trees. You seem to be addressing this in your deprecated draft with "we could now look at phase 1 only at those squirrels that stored nuts at the same number of trees or more". I think this idea is really cool! It looks like you're simulating partial nut measurements with the total nut measurements. This is where my thoughts were going. – timwiz Oct 2 '18 at 20:57
• @timwiz Thanks. You and MartijnWeterings are correct, this is still a problem in the regression / t-n modeling approach. I have an idea how to fix it though, because I believe we can estimate the probability that the squirrel has visited $T_2$ trees in total if we've seen it at $t_{2j}$ trees. As for the deprecated idea: yes indeed that's what I had in mind. I will first try to fix the regression approach and then continue to think about this idea. – tomka Oct 2 '18 at 21:06
• @timwiz I fully revised my answer to take into account the unknown number of trees at phase 2. I added a simulation which demonstrates its performance. I also looked at the bootstrap estimator I disucssed in an earlier version of this answer in work not shown here. This estimator had however worse performance than the estimator I report on and I therefore do not discuss it here (the answer is already quite long and I wanted to report the relevant finding). – tomka Oct 7 '18 at 13:03
• @timwiz Is there anything missing in the answer or open questions? – tomka Oct 8 '18 at 19:52

In phase 1 you could make a model that relates the behavior of a squirrel to the total number of nuts.

In phase 2 you do not observe the exact same full information as in phase 1. For instance, you do not know how many trees in total are used by a specific squirrel. But you do observe some of the squirel's behavior, namely a sample from the distribution of the number of nuts per tree. From this you can estimate the distribution and the parameters that describe the distribution can be input for the model.

• So in phase 1 you make a model that relates the total number of nuts that a squirrel stores with the distribution of the nuts per tree for the squirrel. How exactly to model this is difficult to say.

If you think that you could make some mechanistic model then you could start with some exploratory analysis and prior insights about squirrel behavior to get an idea about a useful model. I lack the data and the biological knowledge to do this in this answer (One obvious direction might be to see whether more nuts per tree will also relate to more nuts in total, and possibly this will have some more complex relationship with the squirrel weight and cheek, and other factors, like high variation in nuts per tree may help to get also an indication about the total number of trees used by a squirrel)

• In phase 2 you will make an estimate of the parameters that are needed to make predictions with the model that has been created in phase 1. The parameters that describe the distribution for the number of nuts per tree can be estimated from the sample measured at the subset of trees.

A simple way would be to ignore the tree id's and just use the data per squirrel to estimate the distribution parameters and put them into the model from phase 1.

A more precise model would treat the tree id's as a random factor such that the behavior specifically attributed to the squirrel can be better estimated. To treat the trees as a random factor you will have to know how the trees can be a random factor. You can make an educated guess for this, but you could also try to learn this from the data (I'd say with some exploratory analysis first, checking out the correlation between a tree and how many nuts get stored in it per squirrel, and whether this effect is independent or maybe some trees attract specific type of squirrels, before coming up with something quantitative.). In the phase 1 you do not observe information related to the tree ids but in phase 2 you do and you can use that data.

So in a nutshell. I think you need some exploratory analysis before you can actually do something quantitative that is more than the simple approach (simple is ignoring the tree ids in phase 2 and using just simple distribution parameters as input for the model whose coefficients are learned in phase 1).

Furthermore, how can we introduce partially measured outputs into the training set?

When you make the model in 1 by using parameters that describe the distribution of the nuts per tree for a specific squirrel you need to take into account that it must be possible to reasonably estimate those numbers in phase 2 and that errors will not effect the model too much. For instance mean and variance (or other simple statistics) could be reasonably estimated from the samples in phase 2 (assuming your sample is not too small) but higher order moments may not.