CODEWITHSUNDEEP
rnon-linear-regression
nicholasflamel
nicholasflamel

Reputation: 33

Run nonlinear (nls) regression that contains summation in R

I'm trying to run a nonlinear regression of the form y_i = sum_t(x_{it}^b)

on the following (simplified) dataset:

require(dplyr)
set.seed(2019)    
df <- data.frame(t = rep(1:4, each = 4),
                 x = sample(1:16, replace=F))

df <- df %>%
group_by(t) %>%
mutate(y = sum(x^2))

Which gives:

       t     x     y
   <int> <int> <dbl>
 1     1    13   396
 2     1    11   396
 3     1     5   396
 4     1     9   396
 5     2     1   626
 6     2    12   626
 7     2    16   626
 8     2    15   626
 9     3    10   361
10     3    14   361
11     3     7   361
12     3     4   361
13     4     8   113
14     4     6   113
15     4     2   113
16     4     3   113

That is, there are 16 uniques values of x but only 4 unique observations of y, and each y is determined by the summation of x^2 that share the same t. I want to run an nls regression along the lines of:

fit <- nls(y ~ sum(x^b), data = df, start=list(b = 2))

With the hope that b will be 2, but I don't know how to write the equation for the fit such that x is summed by groups (of t), instead of being summed altogether.

Thank you.

Upvotes: 3

Views: 246

Answers (1)

LocoGris
LocoGris

Reputation: 4480

I see two things:

  1. For grouping I will use nlsList which it is much more convenient for your model (Fitting nls to grouped data R)

  2. Your model is too much perfect, your "forgot" the error! Because much of the math require to inverse matrices and all this stuff, with your "perfect" model there are issues. Just add a little error! (As we say in Spain: Perfection is fascist!)

    library(nlsList)
set.seed(2019)  
df <- data.frame(t = rep(1:4, each = 4),
                 x = sample(1:16, replace=F))

df <- df %>%
  group_by(t) %>%
  mutate(y = sum(x^2)+runif(4,0,1))

rec.hyp <- nlsList(y ~ sum(x^b) | t,
                   data=df,
                   start=c(b=2),
                   na.action=na.omit)
coef(rec.hyp)

Results

1 2.000428
2 2.000314
3 2.000486
4 2.002057

Upvotes: 2

Related Questions