6
$\begingroup$

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.

enter image description here

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 :)

| cite | improve this question | | | | |
$\endgroup$
  • $\begingroup$ It depends on what kind of relation this 'voltage given an amplification = 150' is. Because of the errors you will get regression dilution (not only a problem for the direction of the relation), but this is not creating a bias if the task is to make a point prediction based on the estimated model (the measurement to which the point estimate applies has the same error and regression dilution). On the other hand, I wonder if the task is really to make those point estimates for the same situation, and the 'voltage at amp=150' may really be more like a sort of parameter estimate. $\endgroup$ – Sextus Empiricus Nov 15 '17 at 22:24
  • $\begingroup$ @MartijnWeterings I'm not sure why regression dilution is not a problem for predicting from the estimated model... Amplification=150 is not necessarily a measurement. I'm also very interested by your comment about this prediction problem being really an estimation problem in disguise - do you mean one of the parameters in the original model may correspond to 'voltage at amp = 150'? For sure it's not the voltage where Amplification becomes infinite - that's larger. Anyway, if you're right, what would you suggest to do? Build a model which explicitly contains that voltage as a parameter? $\endgroup$ – DeltaIV Nov 16 '17 at 7:59
  • $\begingroup$ Like done with growth curves that have halfway points or maximum points (or compare also with Michaelis-Menten kinetics), you could be able to parameterize those curves with some of it's points. For this type of curve it may not be easy to do it explicitly. Yet, it may be done implicitly, if the 150 amp voltage value is more used to characterize the device (for use elsewhere). Saying, 'this is a device with a 150 amp voltage of ...' is the same sort of like saying 'this is a device with a $p_1$ of ....'. $\endgroup$ – Sextus Empiricus Nov 16 '17 at 8:52
  • $\begingroup$ The contrast is with using the 150 amp voltage value, $V_{150}$ , as a prediction for the amplification. The $E(amp|V_{150})$ will be 150 (if you put it in the same setting, that is same error in the measurements of Voltage and amplification, so it is not valid for giving the value to, say, a customer that buys the device and makes different errors ). $\endgroup$ – Sextus Empiricus Nov 16 '17 at 8:53
  • 1
    $\begingroup$ The problem has remained around in my mind. I have to find time to cast it into a nicely formulated form, and also try some things out with the data of the specific problem. $\endgroup$ – Sextus Empiricus Nov 19 '17 at 13:06
6
$\begingroup$

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$$

example of fit with differential equation

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: demonstrating fitting

| cite | improve this answer | | | | |
$\endgroup$
  • 1
    $\begingroup$ Impressive answer (+1). I will study it in detail, but for starters, why do you think the nonlinear model is ill-posed? $\endgroup$ – DeltaIV Nov 25 '17 at 20:49
  • 2
    $\begingroup$ it doesn't fit (at least i cant get it to fit) $\endgroup$ – Sextus Empiricus Nov 25 '17 at 21:32
  • 2
    $\begingroup$ The data clearly indicate asymptotes at amplifiction=1 and Voltage ~ 380. This fixes the coefficient $p_0$ and $p_2$, leaving only the power term $p_1$ to determine the rest of the curve, which does not succeed. Then those $p_0$ and $p_2$, which should have the values of the asymptotes, are gonna be changed in order to make the function still somewhat fitting. The best you get is a reasonable fit, but with the underlying principle (some realistic function) being violated because the fit only makes sense as an empirical relation and not theoretic. More likely is no fit due to bad convergence. $\endgroup$ – Sextus Empiricus Nov 25 '17 at 23:25
  • 1
    $\begingroup$ Voltage ~ 380 is the industrial current voltage... interesting: I wonder if it might be related to the experimental setup. Amplification ~ 1 instead is intuitive - when there's no voltage applied, no amplification results. One last clarification - how do you fit the differential equation with the noise term? Could you please include the code you used to make the fit? Since $a,b,c,d$ are unknown, I think you need to estimate them while numerically integrating the ODE at the same time... $\endgroup$ – DeltaIV Nov 26 '17 at 14:24
  • 2
    $\begingroup$ I have added the code. I changed one thing in the fit which is the weights which were previously defined by a function of the domain A^-2.5 and are now defined by a function of the image (dA/dV)^-1. $\endgroup$ – Sextus Empiricus Nov 28 '17 at 13:43

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.