Understanding Non-linear Regression and its Applications

Non-linear regression Xuhua Xia University of Ottawa http://dambe.bio.uottawa.ca

Non-linear regression All regression analyses are for finding the relationship between a dependent variable (y) and one or more independent variables (x), by estimating the parameters that define the relationship. Functional form known Non-linear relationships whose parameters can be estimated by linear regression: e.g, y = axb, y = abx, y = aebx Non-linear relationships whose parameters can be estimated by non- linear regression, e.g, bx y y ax e + + = = , - ( - ) x 1 Functional form unknown: lowess/loess. While lowess and loess are often treated as synonyms, some people do insist that they are different as prescribed below: lowess: a locally weighted linear least squares regression, generally involving a single IV loess: a locally weighted linear or quadratic least squares regression, involving one or more IVs Slide 2

Rationale of nonlinear regression Both linear and non-linear regression aim to find the parameter values that minimize the residual sum of squared deviation, RSS = [y E(y)]2, or maximize the likelihood function. For linear regression, a solution exists for intercept (a) and slope (b); for non-linear regression, such a solution often does not exist and we need to try various combination of parameter values that either minimize RSS, or maximize lnL Xuhua Xia Slide 4

By using nls Time 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 N 20 42 75 149 278 515 1018 2372 4416 6533 13068 19624 32663 57079 66230 87369 95274 109380 99875 129872 N K K 100000 0 = N t rt + ( ) N N e 0 0 60000 N 20000 0 2 4 6 8 10 Time 10 200000 10 (200000 10) 1.3864 10 200000 278 10 (200000 10) 1.3306 r = Initial values of the parameters to estimate: K: carrying capacity: 200000? N0: 10? r: 1.35? = 20 0.5 + r e = r = 2.5 + r e Xuhua Xia Slide 5

Use EXCEL solver to do estimates Time 0.5 N Pred SS K 200000 These (K, N0, and r) are our guestimates. Now refine them by using EXCEL solver (or by hand if you so wish) with both the LS and ML criteria. 20 19.63938311 0.130044542 42 38.56874494 11.77351128 75 75.73620696 149 148.6941255 0.093559197 278 291.8310018 191.2966097 515 572.3605864 3290.236869 1018 1121.042253 10617.70596 2372 2189.930426 33149.32971 4416 4256.168202 25546.20355 6533 8191.208515 13068 15476.73605 5802009.336 19624 28286.6258 75041085.77 32663 48889.91211 263312676.5 57079 77708.75379 425586741.5 66230 111033.0251 2007311061 87369 142048.5146 2989849314 95274 165601.2467 4945921632 9 109380 180870.7213 5110923225 9.5 99875 189780.4221 8082984924 10 129872 194662.7729 4197844257 N0 10 1 r 1.35 1.5 0.54200069 2 2.5 3 # assume Time and N in R workspace objNLS<-function(y,x,LS1LAD2ML3,expr){ nls(N~N0*K/(N0+(K-N0)*exp(-r*Time)) topfn <- function(df, ex) { eval(substitute(fn(df, ex))) } 3.5 4 4.5 5 2749655.48 5.5 6 6.5 7 7.5 8 fn <- function(dfr, expr) { eval(substitute(expr), dfr) } fn(df, a) fn(df, 2 * a + b) topfn(df, a) topfn(df, 2 * a + b) 8.5 RSS = 28107399387 Xuhua Xia Slide 6

nls (nonlinear LS) md<-read.table("http://dambe.bio.uottawa.ca/teach/bio8102_download/nlinLogistic.txt",header=T) attach(md) fit<-nls(N~N0*K/(N0+(K-N0)*exp(-r*Time)),start=c(K=150000,N0=10,r=1.35)) plot(N~Time) lines(fitted(fit)~Time) Parameters: Estimate Std. Error t value Pr(>|t|) K 1.232e+05 5.412e+03 22.759 3.59e-14 N0 2.708e+01 2.186e+01 1.239 0.232 r 1.151e+00 1.181e-01 9.753 2.23e-08 N K K 0 = N t rt + ( ) N N e 0 0 100000 60000 N N0 is statistically not significantly different from 0, but biologically must be greater than 0. 20000 0 2 4 6 8 10 Xuhua Xia Slide 7 Time

Fitting another equation In rapidly replicating unicellular eukaryotes such as the yeast, highly expressed intron-containing genes requires more efficient splicing sites than lowly expressed genes. GE: gene expression Natural selection will operate on the mutations at the slicing sites to optimize splicing efficiency (SE). Observation: SE increases with GE non-linearly, then levels off and appears to have reached a maximum. GE 1 2 3 4 5 6 7 8 9 10 11 12 13 13 15 16 SE 4.6 4.7 5.7 6.1 6.2 6.8 6.9 7.8 7 7.4 7.7 7.8 7.4 8 8 7.8 8.0 7.5 7.0 6.5 SE + + GE GE 6.0 = SE 1 5.5 5.0 4.5 5 10 15 GE Xuhua Xia Slide 8

Guesstimate initial values The minimum of E(SE) is when GE = 0. 4 + + GE GE = ( ) E SE 1 The maximum of E(SE) is / when GE is large, e.g., 15), / 8, i.e., 8 8 When GE = 6, SE 6.5. 7.5 + 4 (8 )6 1 + 7 = 6.5 0.278 6 8 8*0.278 2.22 6.5 SE Note that the relationship would be linear if = 0. If individual t-tests show that all parameters are all significantly different form 0, then we adopt the non-linear model. However, if a t-test shows that is not significantly different from 0, then one should compare two models, one with = 0 (linear relationship) and one with 0, and use likelihood ratio test or AIC or BIC for model selection (This is because failure to reject a null hypothesis is often not a good indicator that the null hypothesis is true. LRT and information-theoretic indices are better alternatives. 6 5.5 5 4.5 4 0 2 4 6 8 10 12 14 16 18 GE

R functions and output md<-read.table("nlinGESE.txt",header=T) attach(md) fit<-nls(SE~(a+b*GE)/(1+g*GE),start=c(a=4,b=2.22,g=0.278)) summary(fit) plot(GE,SE) lines(fitted(fit)~GE) 8.0 Parameters: Estimate SE t p a 3.33126 0.57853 5.758 0.0000662 *** b 1.92123 0.67520 2.845 0.0138 * g 0.20296 0.08311 2.442 0.0297 * 7.5 7.0 6.5 SE # testing against the simpler model # with g=0 (i.e., linear) fit2<lm(SE~GE) library(lmtest) lrtest(fit, fit2) 6.0 + 3.33126 1.92123 1 0.20296 + GE = SE 5.5 GE 5.0 4.5 5 10 15 GE

A general asymptotic relationship Sometimes we do not know the functional form but we do have some mechanical explanations for the relationship. So here is a general approach. Same problem as before, but we are not sure of the exact relationship between SE and GE GE 1 2 3 4 5 6 7 8 9 10 11 12 13 13 15 16 SE 4.6 4.7 5.7 6.1 6.2 6.8 6.9 7.8 7 7.4 7.7 7.8 7.4 8 8 7.8 8.0 7.5 7.0 6.5 SE 6.0 5.5 5.0 4.5 5 10 15 GE Xuhua Xia Slide 12

A general approach 1. y increases with x at decreasing rate, but stops increasing when x x0. When x < x0, use a polynomial to approximate, e.g., y = a + bx + cx2 2. When x x0), y reaches its maximum and does not increase any more, y = ymax for x x0 2 + + GE if GE GE GE GE GE 0 = ( ) E SE if 8.0 c 0 7.5 7.0 6.5 SE 6.0 5.5 4.5 5.0 5 10 15 GE

Guesstimate initial values + + GE 2 if GE GE GE GE GE GE0 is the GE value beyond which SE has reached maximum (c) and will no longer increase with GE. = 0 ( ) E SE if c 0 8 When GE=0 then SE = , so 4 7.5 For a short segment of GE, the relationship between SE and GE is approximately linear, i.e., SE a + bGE. When GE increases from 2 to 8, SE increases from 4.7 to 7.5, so (7.5-4.7)/(8-2) 0.47 7 6.5 SE 6 Given the linear approximation, with 4 and 0.47, then SE for GE = 12 should be 4+0.47 12 = 9.64, but the actual SE is only about 7.7. This must be due to the quadratic term GE2, i.e., 5.5 5 4.5 (7.7 9.6) = 122, so 4 - 0.02 0 2 4 6 8 10 12 14 16 18 GE Xuhua Xia Slide 14

A few more twists 2 + + GE We will find , , and that minimise RSS = [SE-E(SE)]2 if GE GE GE GE GE 0 = ( ) E SE if c 0 The continuity condition requires that We tell R to substitute various values for , , and , and find the set of values that minimizes RSS 2 0 = = + + GE ( ) when E SE c GE GE GE 0 0 The smoothness condition requires that ( ) E SE GE = + GE = 2 0 0 0 Note that GE0 and c are not parameters because they are functions of , , . = GE 0 2 2 2 0 = + + GE = c GE 0 4 Xuhua Xia Slide 15

EXCEL Solver LS: Alpha Beta Gamma x0 c (Max SE) RSS = ML: Alpha Beta Gamma x0 c (Max SE) lnL 3.921172558 0.604627194 -0.023715164 12.74769174 7.774973098 1.006599054 4.054513849 0.558630253 -0.020633006 13.53729679 7.835685612 17.98555704 Xuhua Xia Slide 16

R statements for LS md<-read.table("nlinGESE.txt",header=T) # Function for estimating the parameters by maximizing lnL # a: alpha, b: beta, g: gamma, x0: GE0 myF <- function(x,dfr) { a<- x[1] b<- x[2] g<- x[3] x0<- -b/2/g c<- a-b^2/4/g seg1Data<-subset(dfr,subset=(dfr$GE < x0)) EY<- a+b*seg1Data$GE+g*seg1Data$GE*seg1Data$GE sumD2<-sum((seg1Data$SE-EY)^2) seg2Data<-subset(dfr,subset=(dfr$GE >= x0)) sumD2<-sumD2 + sum((seg2Data$SE-c)^2) } # obtain solution by supplying the initial values for a, b, g, and the function sol<-optim(c(4,0.47,-0.02),myF,dfr=md) a<-sol$par[1] b<-sol$par[2] g<-sol$par[3] x0<- -b/2/g c<- a-b^2/4/g seg1Data<-subset(md,subset=(md$GE < x0)) EY1<- a+b*seg1Data$GE+g*seg1Data$GE*seg1Data$GE PredY<- c(EY1,rep(c,length(md$GE)-length(seg1Data$GE))) plot(md$GE,md$SE) lines(md$GE,PredY, col="red") abline(v=x0); abline(h=c)

LS Output $par [1] 3.92100405 0.60468725 -0.02371848 , , and 8.0 $value [1] 1.006599 RSS 7.5 $counts function gradient 102 NA 7.0 6.5 md$SE $convergence [1] 0 0 means success 6.0 c [1] 7.775032 5.5 5.0 x0 [1] 12.74718 4.5 5 10 15 md$GE

R statements for ML md<-read.table("nlinGESE.txt",header=T) # Function for estimating the parameters by maximizing lnL # a: alpha, b: beta, g: gamma, x0: GE0 myF <- function(x,dfr) { a<- x[1]; b<- x[2]; g<- x[3]; x0<- -b/2/g; c<- a-b^2/4/g seg1Data<-subset(dfr,subset=(dfr$GE < x0)) EY<- a+b*seg1Data$GE+g*seg1Data$GE*seg1Data$GE Res<-seg1Data$SE-EY lnL<-sum(dnorm(Res,0,sd(Res),log=TRUE)) seg2Data<-subset(dfr,subset=(dfr$GE >= x0)) Res2<-seg2Data$SE-c lnL2<-sum(dnorm(Res2,0,sd(Res2),log=TRUE)) lnL<- -(lnL +lnL2) } # obtain solution by supplying the initial values for a, b, g, and the function sol<-optim(c(4,0.47,-0.02),myF,dfr=md) a<-sol$par[1] b<-sol$par[2] g<-sol$par[3] x0<- -b/2/g c<- a-b^2/4/g seg1Data<-subset(md,subset=(md$GE < x0)) EY1<- a+b*seg1Data$GE+g*seg1Data$GE*seg1Data$GE PredY<- c(EY1,rep(c,length(md$GE)-length(seg1Data$GE))) plot(md$GE,md$SE) lines(md$GE,PredY, col="red") abline(v=x0) abline(h=c)

ML Output $par [1] 3.96764459 0.58297394 -0.02191621 $value [1] -0.0003214071 8.0 $counts function gradient 216 NA 7.5 7.0 $convergence [1] 0 6.5 SE 6.0 c [1] 7.84444 5.5 5.0 x0 [1] 13.3000 4.5 5 10 15 Xuhua Xia GE Slide 20

Robust regression LOWESS: robust local regression between Y and X, with linear fitting LOESS: robust local regression between Y and one or more Xs, with linear or quadratic fitting Used with relations that cannot be expressed in functional forms SAS: proc loess Data: Data set: monthly averaged atmospheric pressure differences between Easter Island, Chile and Darwin, Australia for a period of 168 months (NIST, 1998), suspected to exhibit 12-month (annual), 42-month (El Nino), and 25-month (Southern Oscillation) cycles (From Robert Cohen of SAS Institute) Xuhua Xia Slide 21

lowess in R md<-read.table("nlinGESE.txt",header=T) attach(md) smooth parameter (proportion of data points used): larger = more smooth, default=0.75 fit<-loess(SE~GE,span=0.75,degree=1|2) linear or quadratic, default is 1 summary(fit) pred<-predict(fit,GE,se=TRUE) # OR pred<-predict(fit,c(3,6),se=TRUE) plot(GE,SE) tricubic weighting (proportional to (1 - (dist/maxdist)3)3) lines(GE,pred$fit,col="red") par(mfrow=c(2,3)) for(span in seq(0.4,0.9,0.1)) { How would I know which span value to use? fit<-loess(SE~GE,span=span) pred<-predict(fit,GE) sTitle<-paste0("span = ",span) plot(GE,SE,main=sTitle) lines(GE,pred,col="red") } Xuhua Xia Slide 22

span = 0.4 span = 0.5 span = 0.6 8.0 8.0 8.0 7.5 7.5 7.5 7.0 7.0 7.0 6.5 6.5 6.5 SE SE SE 6.0 6.0 6.0 5.5 5.5 5.5 5.0 5.0 5.0 4.5 4.5 4.5 5 10 15 5 10 15 5 10 15 GE GE GE span = 0.7 span = 0.8 span = 0.9 8.0 8.0 8.0 7.5 7.5 7.5 7.0 7.0 7.0 6.5 6.5 6.5 SE SE SE 6.0 6.0 6.0 5.5 5.5 5.5 5.0 5.0 5.0 4.5 4.5 4.5 5 10 15 5 10 15 5 10 15 GE GE GE

Plotting the fitted values 8.0 fit<-loess(SE~GE,span=0.8) pred<-predict(fit,GE,se=T) pred $fit [1] 4.445761 ... 7.5 7.0 6.5 SE 6.0 $se.fit [1] 0.2785894 ... 5.5 5.0 $residual.scale [1] 0.3273702 4.5 5 10 15 GE $df [1] 10.77648 8 t<-qt(0.975,pred$df) ub<-pred$fit+t*pred$se.fit lb<-pred$fit-t*pred$se.fit plot(GE,SE) lines(GE,pred$fit) lines(GE,lb,col="red") lines(GE,ub,col="red") 7 SE 6 5 plot(GE,SE,ylim=c(min(lb),max(ub))) ... 4 5 10 15 GE