1Physiology 472/572 - 2009 - Quantitative modeling of biological systems Lecture 26: Fitting models to data Fitting a straight line: regression of y on x • given data points (xi,yi), where y is dependent on x • simplest approach is to fit a straight line y = a + bx • fitted values of y: y'i = a + bxi • residuals: ei = yi - y'i = yi - a - bxi • choose a and b to make the residuals as small as possible • they can't all be zero: make sum of squares as small as possible • i.e., minimize 2in1iin1i2i]bxa[yeE −−==∑∑==with respect to a and b • the minimum must be a stationary point • 0xbay]bxa[yn1so0]bxa[y2aEin1iiin1ii=−−=−−=−−−=∂∂∑∑== • 2in1ii])xb(x)y[(yE −−−=∑=, 0)x(x])xb(x)y[(y2bEiin1ii=−−−−−=∂∂∑= • ∑∑==−−−=n1i2iin1ii)x(x)x(x)y(yb and xby a−= x0123456y012345(x i ,yi )ei2Comments • sample variance: 2n1ii)x(x1-n1var(x)∑=−= • sample covariance: )yy)(x(x1-n1y)covar(x,in1ii−−=∑= • then var(x)y)covar(x, b = • regression of x on y: slope is var(y)y)covar(x, b =′ • if the two regression lines were the same, we would expect bb' = 1, but in fact 22ry)var(x)var(y)covar(x, bb ==′where r is the correlation coefficient and r2 ≤ 1 • results depends on whether x or y is considered the independent variable • value of r can be used to test for existence of a correlation • but - always look at the data! - 'Anscombe's quartet' (1973) Property Value x 9.0 var(x) 10.0 y 7.5 var(y) 3.75 R 0.816 regression line y = 3 + 0.5x These properties are the same for all four data sets x0123456y012345regression of x on yregression of y on x3Fitting linear models • suppose x and y are clearly correlated but the relationship is not a straight line • try to fit using a more general relationship • example y = a + bx + cx2 or y = a1f1(x) + a2f2(x) + a3f3(x) + … + amfm(x) • fi may be chosen for convenience or from prior knowledge ('a model') of the system • called a linear model because the dependence on the unknown parameters is linear • the least squares approach can be used • define 2immn1ii11i))]x(fa...)x(fa([yE ++−=∑= and set 0aE...aEm1=∂∂==∂∂ • this gives a linear system of equations that can be solved for a1, …, am • built in to many software packages Fitting nonlinear models • often the dependence on unknown parameters is nonlinear • examples y = axb, P = P0ekt, CKCVVmmax+=, C = Ae-at + Be-bt • two main approaches: transform to a linear model, or solve the nonlinear problem Transform to a linear model • for y = axb, P = P0ekt, take logs: log y = log a + b log x, ln P = ln P0 + kt • for CKCVVmmax+=, maxmaxmV1C1VKV1+= (Lineweaver-Burk plot) • but for C = Ae-at + Be-bt, no such transformation is possible Solving a nonlinear problem • a general nonlinear problem can be written y = f(x, a1, a2, … , am)4• as before, define ∑=′−=n1i2ii]y[yEwhere )a,...,a,a,x(fym21ii=′ • suppose we have some initial estimates of the parameters a10, a20, … , am0 • try to find better estimates a1 = a10 + Δa1, a2 = a20 + Δa2, etc. • mm11m21im0m101iiafa...afa)a,...,a,a,x(f)aa,...,aa,x(fy∂∂Δ++∂∂Δ+=Δ+Δ+≈′ • equations 0aEj=∂∂ can expressed as system of linear equations for Δa1, …, Δam • solve to obtain updated approximations to parameters, repeat until converged • this method is known as Gauss-Newton • often fails because E increases with successive iterations Other methods for nonlinear optimization • modified Gauss-Newton: use aj = aj0 + λΔaj where λ < 1, useful near the solution • gradient descent method: calculate the gradient ⎟⎟⎠⎞⎜⎜⎝⎛∂∂∂∂=∇m1aaE,...,aEE, and update the estimates using j0jjaEaa∂∂+=γwhere γ is a constant that may need to be adjusted • Levenberg-Marquardt method uses a mixture of the above two methods Graph shows contours of E and directions of next step in iteration for each