If in this problem I regress $x$ on $y$ instead than $y$ on $x$, do I need to use an error-in-variables model?

Question

I was trying to write an answer for this question:

Selection of data range changes coefficients too much in lmer (inverse regression)

Basically the OP has lots of data of Amplification vs Voltage (see the image below) for different tested devices and would like to estimate, for each device, the Voltage given an Amplification = 150.

Apparently, in literature there's a physical model which has been used for fitting data from these devices:

$$\text{Amplification}=\frac{1}{1-\left(\frac{\text{Voltage}}{p_0}\right)^{p_1}}+p_2+\epsilon$$

However, the model has Amplification as the dependent variable and Voltage as the predictor. This means that if we fit a nonlinear mixed model to the data with nlmer (for example), then we will get confidence intervals on the Amplification given Voltage, not the other way around.

I thought of simply regressing Voltage on Amplification, instead, i.e., invert the above relationship to get

$$\text{Voltage}=f(\text{Amplification};\theta_1,\theta_2,\theta_3)+\epsilon = \theta_1\left(1-\frac{1}{\text{Amplification}+\theta_2}\right)^{\theta_3} +\epsilon$$

and fit a nonlinear mixed model with nlmer doing something like

model <- function(x, theta1, theta2, theta3) {
  y <- theta1 * (1 - (1 / (x + theta2))) ^ theta3
}

starting_values <- c(theta1 = 1, theta2 = -180, theta3 = 0.8)
nlmer(APD_data, Voltage ~ model(Amplification, theta1, theta2, theta3) ~ 
        theta1 + theta2 + theta3 + theta1|Serial_number + theta2|Serial_number +
        theta3|Serial_number, start = starting_values, verbose = TRUE)

(Serial_number is the tag that identifies each APD device: think of it as a Subject variable).

However, the experiments were conducted (as far as I understood) by controlling Voltage and measuring Amplification. If I invert the two variables in the regression, does it mean I need to build an errors-in-variables model? In that case I give up - nonlinear mixed errors-in-variables models are way above my pay grade :)

Answer 1

The direction of regression may be important to prevent attenuation bias.

Your question about the regression $x \sim y$ versus $y \sim x$ has many angles. A problem which you might encounter is regression attenuation or regression dilution.

• This does not depend on which variable the experiment was controlled in the experiment or in which direction the causal relation has.
• It does depend on the error made in the variables.

It happens if the 'independent variable in the regression' has a large error.

The underlying mathematics does not care about the direction of the causal relationship, or which variable was controlled for, it cares about the errors in the variables.

$$ (y+\epsilon_y) = a + b(x + \epsilon_x) \qquad vs \qquad \frac{(y+\epsilon_y)-a}{b} = (x + \epsilon_x) $$

What we do when changing the "direction" of a regression, is not really changing the direction, but much more like ignoring either $\epsilon_x$ (situation left hand side) or $\epsilon_y$ (situation right hand side).

The variable that was controlled in the experiment, Voltage in this case, can also have a large error, even if it was "controlled". (you do not really set the voltage, you set some button or switch that controls the voltage and eventually you measure the voltage by reading it from a voltmeter or something). So this left hand side situation with the $\epsilon_x$ ignored $(y+\epsilon_y) = a + b(x)$ may be wrong and cause problems (that is attenuation).

'In this problem' you do not have to worry about direction of regression and attenuation.

In this problem there is not a large error, or at least not much noise. The curves are smooth with little jitter. If there is an error then it is a systematic error, but such errors are not linked to regression attenuation (which is about the random errors). Also such systematic errors have little to do with other issues in the mathematics. Except for possibly making some inverse regression ill-posed due to crossing some asymptote or creating negative values in roots logs etcetera.

These systematic errors are more like something that should be dealt with on the practical side (testing equipment, performing good calibration, etcetera).

Much more important is to use the proper model.

1. In the referenced question I have shown how the polynomial model is not working well in the large range https://stats.stackexchange.com/a/315546/164061 .

2. The nonlinear model is not doing much better $$\text{Amplification}=\frac{1}{1-\left(\frac{\text{Voltage}}{p_0}\right)^{p_1}}+p_2+\epsilon$$

   Or at least, I cant seem to make it converge. And I believe it is ill-posed.

3. All this work of trying this perfect fitting is a bit an overkill for the simple task of getting an estimate for the Voltage value at Amplification = 150. This is an interpolation problem, not a fitting problem! This can be done by getting the values of the nearest data points above and below the 150 and use the line between these points to make the estimate. There is no noise which makes this method work badly. If there would be noise then one could use a line trough several points. This interpolation working well, shows that the direction of the relationship is not really the issue here.

If one does wish to fit a reasonable curve then I believe it would be better to dig deeper into the mechanics of the system and use knowledge of the device to create a good curve fit, rather than some polynomial or simplistic literature curve, which are "just" experimental relationships which have little value for generalization and provide little information on the mechanics and inner workings of the devices (which may not be the prime goal of the experiments, but would be a bonus, and at least would increase the robustness of the fitting method).

It is hard to do this work by just gazing at the raw data without much information of the system, however, I was able to get a reasonable fit for a differential equation with two power terms.

$$\frac{\partial A}{\partial V} = a(A-1)^b + c(A-1)^d + \epsilon$$

This is a separable equation and can be solved both numerically (easy, e.g. by deSolve in R) or analytically (although this seems to involve the hypergeometric function $_2F_1$ which is not easy to fit directly).

---

Example of fitting with a differentiated function

Upon request the code to fit with the differentiated function.

I want to stress two points:

• It is indeed an interesting method to fit, although it a whole different question and not necessary (overkill) to make the predictions at Amplification 150. If there is more noise then it may even be be worse than a local linear or polynomial fit in the region around Amplification =150. The fit is only empirically determined and we can not guarantee that it helps to cancel out noise without introducing more bias in return. Some more information about the system (ab initio approach) would certainly help.

• I started out with this method of differentiation, in the first place to investigate and explore potential relationships based on differential equations. To use this method as a final fitting procedure is not always successful and recommendable, because the differentiation will amplify the noise. So it is not a general method to solve the problem and each problem has it's own quirks with a different approach (this also makes the answer to the question about 'regress $x$ on $y$ instead $y$ on $x$' not general and difficult to formulate).

code:

#### demonstrating fit based on differentiated data or differential equation
####
####    note that this is no production code
####    so no checks all sorts and several hard coded limits 


library(deSolve)
library(optimr)

# we make two plots side by side
layout(matrix(c(1,2), 1, 2, byrow = TRUE))

# getting data (for simplicity of the example we just take the 91200913 serial number)

A <- c(1.00252,1.00452,1.00537,1.0056,1.00683,1.0069,1.00847,1.00935,1.01157,1.01418,1.01914,1.0247,1.02919,1.03511,1.04545,1.07362,1.11549,1.17123,1.25019,1.36276,1.5104,1.69862,1.9518,2.26756,2.66278,3.14247,3.73163,4.46152,5.36262,6.49514,7.9227,9.73803,12.0663,15.0943,19.1004,20.0563,21.0672,22.142,23.2867,24.5037,25.8102,27.2024,28.6916,30.2968,32.0181,33.8775,35.8937,38.0569,40.4069,42.9713,45.7766,48.8312,52.2068,55.916,60.0356,64.6109,69.7152,75.4698,82.0003,89.4222,97.9493,107.807,119.441,133.2,149.796,170.113,195.89,229.058,273.481,335.96,431.682,593.091,918.112,1903.74)
V <- c(24.9681,29.9591,34.9494,44.9372,49.9329,54.9625,59.9639,64.9641,69.965,74.9663,79.969,84.9719,89.974,94.9752,99.9759,109.974,119.969,129.96,139.96,149.963,159.958,169.959,179.963,189.957,199.97,209.971,219.971,229.973,239.966,249.962,259.971,269.971,279.97,289.968,299.959,301.968,303.967,305.966,307.965,309.955,311.963,313.963,315.961,317.962,319.956,321.951,323.961,325.962,327.963,329.965,331.966,333.959,335.97,337.972,339.973,341.974,343.967,345.97,347.978,349.98,351.98,353.971,355.983,357.983,359.983,361.975,363.989,365.989,367.984,369.979,371.992,373.994,375.985,377.999)

# numerical differentiation (the noise in this example is not high, alternatively one could higher order methods e.g. a Savitky Golay filter)

dA <- A[-1]-A[-74]
dV <- V[-1]-V[-74]
mA <- 0.5*(A[-1]+A[-74])
mV <- 0.5*(V[-1]+V[-74])

# plotting the derivartive dA/dv with a function of (A-1)
#     substracting A by 1, resulting in A-1 was done 
#     to get an asymptote to zero instead of one.
#     this give probably more interresting relations on a log-log scale

col <- hsv(0,0,0.5)
plot(mA-1, dA/dV, 
     ylim=c(0.000001,2000), xlim=c(0.001,4000),
     pch=21,cex=0.5, col = col, bg = col,
     log="xy",xlab="amplification",ylab="d amplification / d Volt")
title("dA/dV as function of A \n data points and fit")


# this above plot looks so much like two seperate straight lines
#       lets try to fit a double power law
#
#       in the nls fit we use a scaling by the amplification
#       to give more weight on the lower values 
#       (this arbitrary scaling does introduce some subjectivity,
#       but it is just a practical way to avoid the alternatively
#       plotting of logarithms which require the use of the 'port' 
#       algorithm and limits to prevent logs of negative values.
#       so We can do it better, but we ar just lazy for this demonstration)

fit <- nls(dA/dV ~ a*(mA-1)^b+c*(mA-1)^d,
           start=c(a=0.026, b=0.9, c=0.0005, d=1.88),
           weights = (dA/dV)^-1,
           control = nls.control(maxiter = 2000 , minFactor = 10^-9))
# ploting

coefs <- coef(fit)
x <- 10^(c(-10:16)/3)
lines(x,coefs[1]*x^coefs[2]+coefs[3]*x^coefs[4],col=col)

text_string <- paste0("fitted line: dA/dV =",round(coefs[1],5),"A^",round(coefs[2],2)," + ",round(coefs[3],5),"A^",round(coefs[4],2))
text(10^-3,10^-6,text_string,pos = 4)


# going back to the case A vs V
#     - we will need to integrate the fit
#
#     - note that the fit of dA/dV vs A was done to explore the relationship
#       the differentiation is not a very robust operation (noisy) and usually
#       we would whish to obatain a symbolic integration of the differential equation
#       and use the result to relate A as a function of V rather than dA/dV as a function of A
#       
#       but this involves a hypergemoetric function 2F1 which happens to be difficult to fit
#       also the function is not that noisy at all so the fit in the space (dA/dV, A) works well
#
#     - in the end we do add some computation with the optim library to optimize the fit in the 
#       space A vs V, thus without performing differentiation 
#       (or at least, under the hood the optim algorithm will do some differtiation of the Loss function
#       to determine the convergence, but we do not do the noisy differntiation dA/dV)

# plotting data points
plot(V, A,
     xlim=c(0,400), ylim=c(1,4000), log="y",
     pch=21,cex=0.7, col = col, bg = col,
     xlab="Voltage",ylab="Amplification")
title("A as function of V \n data points and two fits")

legend(0,4000,
       c("data points","fit dA/dV vs A","fit A vs V"),
       col=c(8,4,2),
       pt.bg = c(col,0,0),
       lty=c(NA,2,2),
       pch=c(21,NA,NA))

# differential equation for use with deSolve
fitje <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dV <-  -1
    dA <-  -(a*(A-1)^b+c*(A-1)^d)
    list(c(dV, dA))
  })
}

# integration using deSolve (we put this in a function for easier use with optim)
integraal <- function(coefs) {
  parameters <- coefs
  state      <- c(V = max(V), A = max(A))
  times      <- seq(0, 1900, by = 0.1)

  out <- ode(y = state, times = times, func = fitje, parms = parameters)
  out
}

out <- integraal(coefs)

# plot deSolve result
lines(out[,"V"],out[,"A"],
      col=4,lty=2)


# solving the fit in the space V,A by putting the deSolve function into optim  
#         (this will take potentially take some time and is only recommended if we can not
#          find an analytical function to put into optim/nls/nlmer/etc)


f <- function(coefs) {
  A_o <- A
  model <-  integraal(coefs)
  A_m <- approx(x = model[,"V"], y = model[,"A"], xout = V)$y
  sum((A_o-A_m)^2)
}
f(coefs)  

   # fingers crossed that it won't diverge to NAs due to getting out of the range of the approx-function
   #
   # not that using the coefficient from the previous fit will help the convergence
   # performing this method sec without a good initial guess may end badly
model2 <- optim(coefs, f)

#plotting
coefs2 <- model2$par
out2 <- integraal(coefs2)

lines(out2[,"V"],out2[,"A"],
      col=2,lty=2)

Image generated by this code: