0
$\begingroup$

I need to fit a non-linear least squares regression to estimate coefficients which are governed by an equation that have parameters obtained from generating data series.

Data:

q <- c(2524260, 2417498, 2490707, 2376776, 2465064, 2373954, 2410375, 2557140, 2471365, 2523763, 2512031, 2269907, 2504000, 2446307, 2488677, 2385386, 2512441, 2495089, 2414405, 2513172, 2455169, 2525967, 2580019, 2356978, 2623914, 2583873, 2633253, 2560033, 2629156, 2644937, 2625503, 2735719, 2661917, 2770080, 2790052, 2520509, 2827524, 2758730, 2790849, 2660158, 2756222, 2744748, 2721449, 2806078, 2726220, 2812110, 2823104, 2653842, 2824983, 2678610, 2767844, 2628164, 2721210, 2722631, 2626213, 2720187, 2678247, 2746543, 2723628, 2493267, 2797117, 2689967, 2778965, 2688422, 2753997, 2769670, 2763285, 2911373, 2897003, 3025419, 2993178, 2751215, 3069284, 2976108, 3098558, 2972738, 3138242, 3206177, 3153978, 3314676, 3259912, 3391473, 3384552, 3067183, 3396193, 3329471, 3432069, 3317023, 3412343, 3467241, 3398645, 3570623, 3496408, 3620708, 3591672, 3343735, 3563212, 3373797, 3299885, 3226449, 3387987)

p <- c(1.95, 2.43, 2.46, 2.95, 2.84, 2.85, 3.32, 3.54, 3.34, 3.33, 3.33, 3.81, 4.17, 4.04, 3.83, 3.62, 3.43, 3.62, 3.68, 3.64, 4.24, 4.71, 6, 4.9, 4.66, 4.58, 4.59, 4.05, 3.91, 3.92, 3.78, 4.12, 3.48, 2.99, 2.87, 2.83, 2.61, 2.85, 2.78, 2.84, 2.77, 2.66, 2.34, 2.09, 1.93, 2.28, 1.99, 1.73, 1.92, 1.92, 2.59, 2.82, 2.82, 2.99, 2.98, 2.55, 3.59, 3.3, 2.85, 2.88, 3.1, 3.15, 2.98, 2.98, 2.9, 2.98, 2.88, 3.01, 2.82, 3.87, 2.67, 2.69, 2.8, 2.8, 2.97, 2.83, 2.96, 3, 3.28, 4.09, 4.04, 3.11, 2.69, 2.95, 2.65, 2.64, 2.4, 2.37, 2.22, 2.56, 2.33, 2.65, 2.22, 2.02, 1.91, 1.79, 1.74, 1.75, 1.63, 1.77, 2.3)

sigma <- c(0.22, 0.48, 0.03, 0.49, 0.11, 0.01, 0.47, 0.22, 0.2, 0.01, 0.45, 0.48, 0.36, 0.13, 0.21, 0.21, 0.19, 0.19, 0.06, 0.04, 0.6, 0.47, 1.29, 1.1, 0.24, 0.08, 0.01, 0.54, 0.14, 0.01, 0.14, 0.34, 0.64, 0.49, 0.12, 0.04, 0.22, 0.24, 0.07, 0.06, 0.07, 0.11, 0.32, 0.25, 0.16, 0.35, 0.29, 0.26, 0.19, 0.45, 0.67, 0.23, 0.45, 0.17, 0.01, 0.43, 1.04, 0.29, 0.45, 0.03, 0.22, 0.05, 0.17, 0.45, 0.08, 0.08, 0.1, 0.13, 0.19, 1.05, 1.2, 0.02, 0.11, 0.45, 0.17, 0.14, 0.13, 0.04, 0.28, 0.81, 0.05, 0.93, 0.42, 0.26, 0.3, 0.01, 0.24, 0.03, 0.15, 0.34, 0.23, 0.32, 0.43, 0.2, 0.11, 0.12, 0.05, 0.01, 0.12, 0.14, 0.53)

b <- c(2.526, 2.296, 2.048, 2.493, 2.498, 2.963, 2.807, 2.918, 3.5, 3.687, 3.444, 3.35, 3.314, 3.773, 4.161, 4.068, 3.806, 3.641, 3.413, 3.618, 3.654, 3.64, 4.277, 4.542, 5.163, 4.486, 4.608, 4.536, 4.594, 4.025, 3.899, 3.921, 3.801, 4.235, 3.509, 2.929, 2.755, 2.747, 2.591, 2.856, 2.769, 2.805, 2.753, 2.639, 2.378, 2.281, 2.044, 2.233, 1.929, 1.812, 2.014, 2.083, 2.634, 2.761, 2.722, 2.903, 3.064, 2.875, 3.591, 3.291, 2.907, 2.992, 3.189, 3.236, 2.994, 2.955, 2.905, 3.006, 2.916, 3.064, 2.776, 3.153, 2.657, 2.701, 2.724, 2.834, 2.943, 2.793, 2.908, 2.897, 3.208, 4.11, 3.903, 3.11, 2.68, 2.805, 2.598, 2.593, 2.33, 2.303, 2.174, 2.508, 2.339, 2.625, 2.283, 2.03, 1.84, 1.73, 1.76, 1.81, 1.7)

rho <- c(0.05, 0.23, 0, 0.24, 0.01, 0, 0.22, 0.05, 0.04, 0, 0.2, 0.23, 0.13, 0.02, 0.04, 0.04, 0.04, 0.04, 0, 0, 0.36, 0.22, 1.66, 1.21, 0.06, 0.01, 0, 0.29, 0.02, 0, 0.02, 0.12, 0.41, 0.24, 0.01, 0, 0.05, 0.06, 0, 0, 0, 0.01, 0.1, 0.06, 0.03, 0.12, 0.08, 0.07, 0.04, 0.2, 0.45, 0.05, 0.2, 0.03, 0, 0.18, 1.08, 0.08, 0.2, 0, 0.05, 0, 0.03, 0.2, 0.01, 0.01, 0.01, 0.02, 0.04, 1.1, 1.44, 0, 0.01, 0.2, 0.03, 0.02, 0.02, 0, 0.08, 0.66, 0, 0.86, 0.18, 0.07, 0.09, 0, 0.06, 0, 0.02, 0.12, 0.05, 0.1, 0.18, 0.04, 0.01, 0.01, 0, 0, 0.01, 0.02, 0.28)

delta <- c(0.35, 0.38, 0.03, 0.45, 0.12, 0.04, 0.4, 0.04, 0.35, 0.11, 0.02, 0.5, 0.4, 0.12, 0.24, 0.19, 0.21, 0.21, 0.06, 0.01, 0.6, 0.43, 1.46, 0.26, 0.17, 0.03, 0.05, 0.54, 0.12, 0.02, 0.14, 0.32, 0.76, 0.52, 0.06, 0.08, 0.14, 0.26, 0.08, 0.07, 0.04, 0.09, 0.3, 0.29, 0.35, 0.24, 0.24, 0.2, 0.11, 0.09, 0.51, 0.19, 0.06, 0.27, 0.08, 0.51, 0.71, 0.29, 0.44, 0.03, 0.11, 0.04, 0.26, 0.01, 0.06, 0.08, 0.13, 0.09, 0.24, 1.09, 0.48, 0.03, 0.1, 0.08, 0.14, 0.11, 0.17, 0.09, 0.38, 0.88, 0.07, 0.79, 0.42, 0.27, 0.16, 0.04, 0.19, 0.04, 0.08, 0.39, 0.18, 0.31, 0.4, 0.26, 0.12, 0.05, 0.01, 0.01, 0.18, 0.07, 0.53)

I proceed to create initial values and attempt to run an nls model:

initial = list(
Va = ?
Vgamma_p = ?
Vgamma_sigma = ?
Vgamma_b = ?
Vgamma_rho = ?
Vgamma_delta = ?
Va_p = ?
Va_sigma = ?
Va_b = ?
Va_rho = ?
Va_delta = ?
)  

model = nls(q ~ (Va + Vgamma_p*p + Vgamma_sigma*sigma + Vgamma_b*b + Vgamma_rho*rho + Vgamma_delta*delta)/((1 + Va_p + Va_sigma*sigma + Va_b*b + Va_rho*rho + Va_delta*delta)*p),start=initial)

I've assigned question marks to the initial values for my coefficients here because I honestly have no idea what values I should use. I've tried a lot of different starting values, but all give me some error. One common error I am getting is: "step factor 0.000488281 reduced below 'minFactor' of 0.000976562"

I've tried to work around that error by playing with the nls.control function, but I am not entirely sure if this is appropriate:

model = nls(q ~ (Va + Vgamma_p*p + Vgamma_sigma*sigma + Vgamma_b*b + Vgamma_rho*rho + Vgamma_delta*delta)/((1 + Va_p + Va_sigma*sigma + Va_b*b + Va_rho*rho + Va_delta*delta)*p),start=initial, control = nls.control(maxiter = 10^6, minFactor = 10^-20))

I've received many more kinds of errors depending on which initial values I've used. These include: "Missing value or an infinity produced when evaluating the model", ""singular gradient matrix at initial parameter estimates" and "NA/NaN/Inf in foreign function call (arg 1)"

I know that the source of the problem must be the choice of initial values, but I've been stuck on this for a few days trying to figure out what to do.

It may be suggested that the equation does not provide a good fit for q, but it is essential that I must use this equation. I also tried differencing q (which was a suggestion made by my supervisor), but I still end up with these problems.

Any assistance would be really appreciated.

$\endgroup$
0
2
$\begingroup$

A few points:

  • we can avoid having to specify the parameters which enter linearly (i.e. all the ones in the numerator) by using the plinear algorithm. In that case specify the numerator as a matrix such that the ith column multiplies the ith linear parameter.

  • the Va_p parameter is causing the problem making the solution non-identifiable and over-parameterized. If we multiply every numerator parameter by k and multiply every denominator parameter also by k except transform Va_p to k*(Va_p+1)-1 then this multiplies both numerator and denominator by k giving the same predictions for every non-zero value of k thus we need to drop the Va_p term from the denominator.

  • An alternate way to implement the above is to use the same initial and nls statement as shown in the question but drop Va_p from initial while setting all other initial values to 1 and fixing Va_p to 0, i.e. Va_p <- 0, before running nls.)

The approach outlined in the first two points lead to significantly shorter code as follows:

initial <- list(Va_sigma = 1, Va_b = 1, Va_rho = 1, Va_delta = 1)  
nls(q ~ cbind(1, sigma, b, rho, delta) /
  ((1 + Va_sigma * sigma + Va_b * b + Va_rho * rho + Va_delta * delta) * p),
  start = initial, algorithm = "plinear")

giving:

Nonlinear regression model
  model: q ~ cbind(1, sigma, b, rho, delta)/((1 + Va_sigma * sigma + Va_b *     b + Va_rho * rho + Va_delta * delta) * p)
   data: parent.frame()
  Va_sigma       Va_b     Va_rho   Va_delta      .lin1 .lin.sigma     .lin.b 
   4.41090    0.21459    0.59320   -4.28024   -0.01483    7.80228    1.27647 
  .lin.rho .lin.delta 
   4.64760   -9.20623 
 residual sum-of-squares: 1.278

Number of iterations to convergence: 18 
Achieved convergence tolerance: 6.61e-06
$\endgroup$

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.