Identifiability notebook

Let's start by just fitting all the subjects with model.diff to get a sense of what reasonable parameter values are.

fits.diff = ordf(mse, `_data` = ddply(obsts, .(s), function(slice)
   {f = fit.model(slice, model.diff)
    as.data.frame(t(c(f$par, mse = f$mse)))}))
row.names(fits.diff) = char(fits.diff$s)
fits.diff = ss(fits.diff, select = -s)
rd(rbind(fits.diff,
    mapcols(fits.diff, function(v) quantile(v, c(.025, .975)))))

All right, now, which portions of the parameter space are "interesting", in the sense of not resulting in degenerate choice patterns for the Builder test set?

llsgrid.diff = cached(local(
   {grid = expand.grid(
        kappa = exp(seq(-20, 10, len = 150)),
        rho = exp(seq(-15, 20, len = 150)))
    ddply(grid, .(kappa, rho), .progress = "text", function(row)
       {pd = preddata(builder.test.ts, model.diff, num(row))
        data.frame(expected.lls = mean(pd$pll), sdilogit = sd(ilogit(pd$pll)))})}))

First, a graph of expected proportion of LLs chosen.

ggplot(llsgrid.diff) +
    geom_tile(aes(log(kappa), log(rho), fill = expected.lls)) +
    scale_fill_gradient2(midpoint = .5) +
    scale_x_continuous(breaks = seq(-20, 10, 5), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(-15, 20, 5), expand = c(0, 0))

Now a figure showing the standard deviation of the inverse logits of the choice probabilities across the test set. The higher this number, the more heterogeneity there is across the test set; when it's 0, the choice probability is constant.

ggplot(llsgrid.diff) +
    geom_tile(aes(log(kappa), log(rho), fill = sdilogit)) +
    scale_fill_continuous(low = "black", high = "white", limits = c(0, .12), na.value = "red") +
    scale_x_continuous(breaks = seq(-20, 10, 5), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(-15, 20, 5), expand = c(0, 0))

So we have roughly three regions where the model is interesting enough that checking for mimicry makes sense, corresponding to the white bar and the two gray bars. Or we could just examine the intersection point and the space around it. Yeah, that's what I'll do.

msegrid = function(fitting.model, n.replicates = 5)
   {grid = expand.grid(
        kappa = exp(seq(-7, 7, len = 150)),
        rho = exp(seq(-8, 5, len = 150)))
    cbind(grid, mse = aaply(1 : nrow(grid), 1, .progress = "text", function(i)
        median(replicate(n.replicates,
              # I take the median of MSEs from several replicates,
              # instead of just using one, in order to reduce
              # noise in the plots.
            fit.model(model = fitting.model, ts =
                gendata(builder.test.ts, model.diff, num(grid[i,])))$mse))))}

plot.msegrid = function(g) ggplot(g) +
    geom_tile(aes(log(kappa), log(rho), fill = mse)) +
    scale_x_continuous(breaks = -7:7, expand = c(0, 0)) +
    scale_y_continuous(breaks = -8:5, expand = c(0, 0)) +
    inset.legend(c(1, 1))

msegrid.diff = cached(msegrid(model.diff))
msegrid.ghlp = cached(msegrid(model.ghlp))

plot.msegrid(msegrid.diff) +
     scale_fill_continuous(
         high = "black", low = "white", na.value = "red",
         limits = c(0, .6))

This seems to generally make sense. The upper part of the plot is white because rho is high and therefore there's little noise in the data; the lower part is dark gray because rho is low and therefore the data is mostly noise. I'm not sure what to make of the very pale gray spike.

plot.msegrid(msegrid.ghlp) +
     scale_fill_continuous(
         high = "black", low = "white", na.value = "red",
         limits = c(0, .6))

This looks weirder, which makes sense considering that the generating and fitting models no longer match. In general, ghlp has trouble with that central column.

Now for the difference grid.

plot.msegrid(transform(msegrid.ghlp, mse = mse - msegrid.diff$mse)) +
     scale_fill_gradient2() +
     geom_point(data = data.frame(fits.diff),
        aes(log(kappa), log(rho)), shape = "x", size = 5) +
     coord_cartesian(xlim = c(-7, 7))

model.diff generally does better, particularly in the central column. It's nice that the parameter vectors estimated from the subjects are in regions of reasonable favor towards diff. (Note that one subject, whose estimated kappa was very small, is omitted.)

Again, let's start with point estimates for the real data.

fits.ghlp = ordf(mse, `_data` = ddply(obsts, .(s), function(slice)
   {f = fit.model(slice, model.ghlp)
    data.frame(
        log.curve = log(f$par["curve"]),
        log.k = log(f$par["k"]),
        log.rho = log(f$par["rho"]),
        mse = f$mse)}))
row.names(fits.ghlp) = char(fits.ghlp$s)
fits.ghlp = ss(fits.ghlp, select = -s)
rd(d = 2, rbind(fits.ghlp,
    mapcols(fits.ghlp, function(v) quantile(v, c(.025, .975))),
    mapcols(fits.ghlp, range)))

Now, I can't very well use heatmaps and so on with a three-dimensional parameter space, so we have to more or less drop a parameter. Christian suggested dropping rho, since it is theoretically the least important parameter; makes sense to me. Let's freeze it at the median (which is the same thing as the antilog of the median of the logs).

rho.med = exp(median(fits.ghlp$log.rho))
c(
   "median of logs" = median(fits.ghlp$log.rho),
   "median" = rho.med)

Now let's get point estimates again with rho fixed at rho.med:

model.ghlp.mr = model.ghlp.frho(rho.med)
fits.ghlp.mr = ordf(mse, `_data` = ddply(obsts, .(s), function(slice)
   {f = fit.model(slice, model.ghlp.mr)
    data.frame(
        log.curve = log(f$par["curve"]),
        log.k = log(f$par["k"]),
        mse = f$mse)}))
row.names(fits.ghlp.mr) = char(fits.ghlp.mr$s)
fits.ghlp.mr = ss(fits.ghlp.mr, select = -s)
rd(d = 2, rbind(fits.ghlp.mr,
    mapcols(fits.ghlp.mr, function(v) quantile(v, c(.025, .975))),
    mapcols(fits.ghlp.mr, range)))

Removing outliers:

d = ss(fits.ghlp.mr, (log.curve < 100 | log.k < 100) & is.finite(log.k))
rd(d = 3, rbind(
    mapcols(d, function(v) quantile(v, c(.025, .975))),
    mapcols(d, range)))

llsgrid.ghlp.mr = cached(local(
   {grid = expand.grid(
        curve = exp(seq(-25, 10, len = 150)),
        k = exp(seq(-8, 4, len = 150)))
    ddply(grid, .(curve, k), .progress = "text", function(row)
       {pd = preddata(builder.test.ts, model.ghlp.mr, num(row))
        data.frame(expected.lls = mean(pd$pll), sdilogit = sd(ilogit(pd$pll)))})}))

ggplot(llsgrid.ghlp.mr) +
    geom_tile(aes(log(curve), log(k), fill = sdilogit)) +
    scale_fill_continuous(low = "black", high = "white", limits = c(0, .12), na.value = "red") +
    scale_x_continuous(breaks = seq(-25, 10, 5), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(-8, 4, 2), expand = c(0, 0))

msegrid.mr = function(fitting.model, n.replicates = 5)
   {grid = expand.grid(
        curve = exp(seq(-25, 10, len = 150)),
        k = exp(seq(-8, 4, len = 150)))
    cbind(grid, mse = aaply(1 : nrow(grid), 1, .progress = "text", function(i)
        median(replicate(n.replicates,
              # I take the median of MSEs from several replicates,
              # instead of just using one, in order to reduce
              # noise in the plots.
            fit.model(model = fitting.model, ts =
                gendata(builder.test.ts, model.ghlp.mr, num(grid[i,])))$mse))))}

plot.msegrid.mr = function(g) ggplot(g) +
    geom_tile(aes(log(curve), log(k), fill = mse)) +
    scale_x_continuous(breaks = seq(-25, 10, 5), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(-8, 4, 2), expand = c(0, 0)) +
    inset.legend(c(1, 0))

msegrid.mr.diff = cached(msegrid.mr(model.diff))
msegrid.mr.ghlp = cached(msegrid.mr(model.ghlp))

plot.msegrid.mr(msegrid.mr.ghlp) +
     scale_fill_continuous(
         high = "black", low = "white", na.value = "red",
         limits = c(0, .3), breaks = seq(0, .3, .1))

plot.msegrid.mr(msegrid.mr.diff) +
     scale_fill_continuous(
         high = "black", low = "white", na.value = "red",
         limits = c(0, .6))

The difference grid:

plot.msegrid.mr(transform(msegrid.mr.diff, mse = mse - msegrid.mr.ghlp$mse)) +
     scale_fill_gradient2(limits = c(-.1, .4)) +
     geom_point(data = data.frame(fits.ghlp.mr),
        aes(log.curve, log.k), shape = "x", size = 5) +
     coord_cartesian(xlim = c(-25, 5), ylim = c(-8, 2))

Each of the previous studies compared two discounting models to each other. In this study, we compared the most general of the discounting models, GH, with an attribute-based model, Diff. One might expect that even if discounting models are highly capable of mimicking each other, a model that uses an entirely different approach to intertemporal choice may not be so capable of mimicking a discounting model, nor vice versa.

We parametrized GH with the log-odds of choosing LL set to

$$\rho [r_L(1+\alpha d_L{)}^{-\beta /\alpha }-r_S(1+\alpha d_S{)}^{-\beta /\alpha }]$$

with α ≥ 0, β > 0, and ρ > 0. We parametrized Diff with the log-odds of choosing LL set to

$$\rho [r_L-r_S-\kappa (d_L-d_S)]$$

with κ > 0 and ρ > 0. We estimated model parameters by minimizing MSE with GNU R's function optim using its default Nelder-Mead method, along with some explicit logic to ensure that parameters were chosen to obtain MSE 0 when the input data was all-SS or all-LL. As our set of quartets for parameter-space exploration, we used Arfer and Luhmann's (2015) standard set of 100.

fig--g/pb-densities shows density plots of parametrically bootstrapped model fit (MSE) for each subject and each generating model. For most subjects and models, the mean fit favors the generating model. (The exceptions are subjects s53 and s61, for whom Diff achieves better mean fit when GH is the generating model, and subjects s50, s62, s67, s71, s72, and s75, for whom the reverse holds.) Still, one can see substantial visible overlap between bootstrap distributions. Collapsing across subjects, we find that when GH is the generating model, it achieves better fit than Diff in 80% of replicates, whereas when Diff is the generating model, it achieves better fit than GH in 68% of replicates. Thus, discrimination on the basis of model fit usually decides correctly, but has a non-negligible error rate, and is biased towards choosing GH. This bias is consistent with the models' different numbers of parameters: GH has 3 parameters whereas Diff has 2.

To explore mimicry with Diff as the generating model across its parameter space, we first needed to choose the rectangular section of the parameter space to inspect. Per-subject point estimates for original data had $\log \rho$ ranging from −8 to 1 and $\log \kappa$ ranging from −6 to 7, excepting one subject (s50) with an estimated $\log \kappa$ of −16. fig--g/diff-sd is a heatmap of the SD of inverse logits of choice probabilities across the Arfer and Luhmann (2015) quartets given Diff with the given parameter pair as the generating model. The darker each pixel of the heatmap, the closer the data is to being degenerate, with all choice probabilities being equal. Thus, the nontrivial portions of the parameter space (relative to the Arfer and Luhmann quartets) are given by the lighter regions. On this basis, we decided to inspect $\log \rho$ from −8 to 5 and $\log \kappa$ from −7 to 7. fig--g/dmsegrid is a heatmap of the difference of model fits of GH and Diff for data generated by Diff with each parameter pair, where positive MSE indicates greater fit by Diff. The original estimate for each subject is marked by a cross (all subjects are visible except s50). Observe that in nontrivial portions of the parameter space, Diff indeed achieves better fit, and subjects tend to cluster in areas where Diff achieves better fit.

To explore mimicry with GH as the generating model across its own parameter space, we again needed to choose a rectangular section of the parameter space to inspect. Since GH has 3 parameters, this required dropping a dimension. We opted to drop ρ because it is the least theoretically important parameter. We did this by making per-subject point estimates with all parameters free to vary, and calculating the median of the estimated ρ values, which was 0.274. We then redid these estimates with ρ fixed at the median, yielding $\log \alpha$ ranging from −22 to 2 and $\log \beta$ ranging from −6 to 0, excepting one subject (s50) with estimates of 137 and 129 (respectively) and another subject (s61) with a $\log \beta$ estimate of −∞ (i.e., with β = 0). fig--g/ghlp-mr-sd is a heatmap similar to fig--g/diff-sd, and fig--g/dmsegrid-mr is similar to fig--g/dmsegrid but with positive values meaning greater fit for GH. (Note that although all data was generated with a single value of ρ, ρ was free to vary when GH was fitted for fig--g/dmsegrid-mr.) The overall impression is similar to that given in the parameter-space exploration of Diff, but stronger. The only large advantages in terms of fit are in GH's favor, and subjects' estimates tend to be in regions of such advantages.

In summary, the parametric-bootstrapping analyses as well as the two parameter-space explorations all consistently suggest that Diff and GH are not always easy to distinguish, but on average, model fit strongly favors the generating model.