0

I am dealing with the relationship:

y = h * R + x * v * h

where:

x = (N - M) * exp(-Q * u) + M

which gives the principal equation:

y = h * R + v * h * (N - M) * exp(-Q * u) + v * h * M

All uppercase letters are constants, and all lowercase letters are variables.

I have real data for all the variables, but I either do not know the values of the constants (R and Q), or want to check the ability of the data to fit the values of the constants (N and M). I want to use nls() to fit the equation using the data for the variables, to estimate these constant parameters.

How do I write code using the nls() function to depict the principal equation, to allow estimation of the parameters R, N, Q, and M from the mock measurement data? (Mock measurement data = lower cases letters with _j suffix, see below.)

To create mock data:

library(dplyr)
library(ggplot2)

### Generate mock data

# Equations:
# y = h*R + x*v*h
# x = (N-M)*exp(-Q*u) + M
# y = h*R + ((N-M)*exp(-Q*u) + M)*v*h
# y = h*R + v*h*(N-M)*exp(-Q*u) + v*h*M

### Variables have varying periodicity,
# and so can be approximated via different functions,
# with unique noise added to each to simulate variability:

# Variability for each variable
n <- 1000 # number of data points
t <- seq(0,4*pi,length.out = 1000)
a <- 3
b <- 2

y.norm <- rnorm(n)
u.norm <- rnorm(n)
u.unif <- runif(n)
v.norm <- rnorm(n)
v.unif <- runif(n)
amp <- 1

# Create reasonable values of mock variable data for all variables except h;
# I will calculate from known fixed values for R, N, Q, and M.

y <- 1.5*a*sin(b*t)+y.norm*amp-10 # Gaussian/normal error
u <- ((1*a*sin(11*b*t)+u.norm*amp)+(0.5*a*sin(13*b*t)+u.unif*amp)+7)/2
v <- 1/((2*a*sin(11*b*t)+v.norm*amp)+(1*a*sin(13*b*t)+v.unif*amp)+20)*800-25

# Put vectors in dataframe
dat <- data.frame("t" = t, "y" = y, "u" = u, "v" = v)

### Create reasonable values for constants:

R=0.5
N=1.12
Q=0.8
M=1

### Define final variable based on these constants and the previous
# mock variable data:

dat$h = y/(R + v*(N-M)*exp(-Q*dat$u))

### Gather data to plot relationships:

dat_gathered <- dat %>%
  gather(-t, value = "value", key = "key")

### Plot data to check all mock variables:

ggplot(dat_gathered, aes(x = t, y = value, color = key)) + geom_line()

# Add small error (to simulate measurement error):

dat <- dat %>%
  mutate(h_j = h + rnorm(h, sd=0.05)/(1/h)) %>%
  mutate(u_j = u + rnorm(u, sd=0.05)/(1/u)) %>%
  mutate(v_j = v + rnorm(v, sd=0.05)/(1/v)) %>%
  mutate(y_j = y + rnorm(y, sd=0.05)/(1/y))
photosynthesis
  • 95
  • 1
  • 11

1 Answers1

1

nls appears to work OK, but it looks like the solution (in terms of parameters) is non-unique ... or I made a mistake somewhere.

## parameter values chosen haphazardly
n1 <- nls(y ~ h_j*(R + v_j*((N-M)*exp(-Q*u_j)+M)),
    start=list(R=1,N=2,M=1,Q=1),
    data=dat)

## starting from known true values
true_vals <- c(R=0.5,N=1.12,Q=0.8,M=1)
n2 <- update(n1, start=as.list(true_vals))

round(cbind(coef(n1),coef(n2),true_vals),3)
              true_vals
R 0.495 0.495      0.50
N 0.120 0.120      1.12
M 0.001 0.818      0.80
Q 0.818 0.001      1.00

Using AIC() on the two fits shows they have essentially equivalent goodness of fits (and the predictions are almost identical), which suggests that there's some symmetry in your model that allows M and Q to be interchanged. I haven't thought about/looked at the equation hard enough to know why this would be the case.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453