# Fitting a non-linear least squares regression to estimate coefficients of parameters in R

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.

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