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:

  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 +

> summary(daffit)

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

   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.

  • $\begingroup$ 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? $\endgroup$
    – user13907
    Oct 2 '14 at 12:44
  • $\begingroup$ 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. $\endgroup$
    – Glen_b
    Oct 2 '14 at 16:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.