# How can I optimize coefficients of an arbitrary model?

This might be terribly easy but I'm probably lacking the keywords to search for. Assume the following (dummy) data:

condition   session assumed_score  final_score
1         1           110           95
1         2           90            95
1         3           80            95
1         4           120           95
1         5           130           95
2         1           80            77
2         2           90            77
2         3           80            77
2         4           50            77
2         5           90            77
3         1           60            35
3         2           40            35
3         3           30            35
3         4           30            35
3         5           40            35


In each testing condition, there are people taking part in 5 (or more, or less) sessions. They perform a task in each session. After all sessions, they give a final score (final_score) to their tasks. This is the ground truth. But we also (based on some measure) have assumed scores for each of the individual test sessions (assumed_score), based on observation of the test participants.

Now, I would like to be able—in the future—to predict the final score from the assumed scores. So, one "dumb" model would be to just average them:

$$avg_i = \frac{1}{N} \sum assumed_{ni}$$

where $N$ is the number of sessions and $i$ is the condition.

So now I have the averages, and I can build a simple linear model:

$$final = \beta_0 + \beta_1 avg$$

I could do this in R like so, and get a model that should fit the assumed scores to the final scores with least squares:

d = ddply(d, .(condition), transform, avg = mean(assumed_score))
lm(final_score ~ avg, data = d)


This works, because there are no coefficients in the average.

But what if I wanted a more complex model? Say I'd like to give more weight on the sessions that have a lower score. Then my model would be:

$$final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}}$$

So here, I'd need to find an optimal parameter $\beta$ to fit my final scores against the assumed scores. But I have no idea how to do that with R.. I could try a fixed $\beta$, and then do:

custom <- function(d, beta = 0.1) {
N = length(d)
sum = 0
for (a in d) {
sum = sum + a^beta
}
avg = sum / N
avg = avg^(1/beta)
}
d = ddply(d, .(condition), transform, custom = custom(assumed_score))
lm(final_score ~ custom, data = d) # => nice, this has lower R^2!


But that's hardly a good option to find a good value of $\beta$ without trying it manually.

How should I proceed?

One problem is you haven't specified the error structure.

Your model isn't actually $final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}}$. The observed $y$-values won't equal the right hand side -- if that were the case you wouldn't need statistics.

You might mean something like:

$final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}}+\epsilon$

where $\text{Var}(\epsilon)=\sigma^2I$.

In which case you probably want nonlinear least squares.

You would probably want to put your data in this (wide) form:

dw
condition assumed_s1 assumed_s2 assumed_s3 assumed_s4 assumed_s5 final
1         1        110        120         90         90         30    95
2         2         90        130         80         60         30    77
3         3         80         80         50         40         40    35


since then

daffit <- nls(final ~ b0 +
((assumed_s1^b+assumed_s2^b+assumed_s3^b+assumed_s4^b+assumed_s5^b)/5)^(1/b),
data=dw,start=list(b0=0,b=1))

> summary(daffit)

Formula: final ~ b0 + ((assumed_s1^b + assumed_s2^b + assumed_s3^b
+ assumed_s4^b + assumed_s5^b)/5)^(1/b)

Parameters:
Estimate Std. Error t value Pr(>|t|)
b0  -17.232     31.944  -0.539    0.685
b     3.855     10.441   0.369    0.775

Residual standard error: 18.7 on 1 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 1.586e-06


If you know roughly what $\beta_0$ or $\beta$ should be (perhaps from some simpler or approximate model fit, for example, or from experience of typical values), use them as start values.

It's still possible to do it with a different number of sessions per condition ($N_i$ rather than $N$).

With your particular function, you could still do it with nls since, for example, (if you keep $\beta>0\$**) $0^\beta=0$, you can just add a column for $N_i$, pad out the number of terms to the maximum (with 0), then put $(\sum\{...\}/N_i)^{(1/\beta)}$ (where "$...$" has a sum of terms like now).

** That can be done fairly simply via a reparameterization.

Some more complicated functions might require using more general optimization functions like optim, but you can get pretty clever with nls (cleverer than I have shown here).

Of course, none of this is really specific to R - many other stats packages have general optimization and/or flexible nonlinear least squares.

• Ah! But that means I can't have a variable number of N? That was the idea of using a summation, so that it'd work in all cases, even those where my model would use N as input. Do you think this is possible to achieve, somehow? Oct 2 '14 at 12:44
• Yes. I answered here but it was pretty long so I moved it up to the answer. See the update at the end of my answer. Oct 2 '14 at 16:24