<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/second-step-with-non-linear-regression-adding-predictors/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/second-step-with-non-linear-regression-adding-predictors/"> R Programming – DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In this post we will see how to include the effect of predictors in non-linear regressions. In other words, letting the parameters of non-linear regressions vary according to some explanatory variables (or predictors). Be sure to check the <a rel="nofollow" href="https://datascienceplus.com/first-steps-with-non-linear-regression-in-r/" target="_blank">first post on this</a> if you are new to non-linear regressions.</p>
<p>The example that I will use throughout this post is the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Logistic_function#In_ecology:_modeling_population_growth" target="_blank">logistic growth function</a>, it is often used in ecology to model population growth. For instance, say you count the number of bacteria cells in a petri dish, in the beginning the cell counts will increase exponentially but after some time due to limits in resources (be it space or food), the bacteria population will reach an equilibrium. This will produce the classical S-shaped, non-linear, logistic growth function. The logistic growth function has three parameters: the growth rate called “r”, the population size at equilibrium called “K” and the population size at the beginning called “n0”.</p>
<h2>First Example: Categorical Predictor</h2>
<p>Say that we want to test the effect of food quality on the population dynamic of our bacteria. To do so we use three different type of food and we count the number of bacteria over time. Here is how to code this in R using the function <code>nlsList</code> from the package <a rel="nofollow" href="https://cran.r-project.org/web/packages/nlme/index.html" target="_blank">nlme</a>:</p>
<pre>
#load libraries
library(nlme)
#first try effect of treatment on logistic growth
Ks &lt;- c(100,200,150) 
n0 &lt;- c(5,5,6)
r &lt;- c(0.15,0.2,0.15)
time &lt;- 1:50
#this function returns population dynamics following
#a logistic curves
logF &lt;- function(time,K,n0,r){
  d &lt;- K * n0 * exp(r*time) / (K + n0 * (exp(r*time) - 1))
  return(d)
}
#simulate some data
dat &lt;- data.frame(Treatment=character(),Time=numeric(),
                  Abundance=numeric())
for(i in 1:3){
  Ab &lt;- logF(time = time,K=Ks[i],n0=n0[i],r=r[i])
  tmp &lt;- data.frame(Treatment=paste0(&quot;T&quot;,i),Time=time,
                    Abundance=Ab+rnorm(time,0,5))
  #note that random deviates were added to the simulated 
  #population density values
  dat &lt;-rbind(dat,tmp)
}
#the formula for the models
lF&lt;-formula(Abundance~K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1)) | Treatment)

#fit the model
(m &lt;- nlsList(lF,data=dat,start=list(K=150,n0=10,r=0.5)))
<em>Call:
Model: Abundance ~ K * n0 * exp(r * Time)/(K + n0 * (exp(r * Time) - 1)) | Treatment Data: dat Coefficients:
K n0 r
T1 105.2532 4.996502 0.1470740
T2 199.5149 4.841359 0.2006150
T3 150.2882 5.065196 0.1595293
Degrees of freedom: 150 total; 141 residual
Residual standard error: 5.085229
</em>
</pre>
<p>The code simulated population values using three sets of parameters (the r, K and n0’s). Then we specified the non-linear regression formula, using the pipe “|” symbol to explicitly ask for fitting different parameters to each Treatment. The model output gives us the estimated parameters for each Treatment. We can now plot the model predictions against the real data:</p>
<pre>
#derive the predicted lines
dat$Pred &lt;-predict(m)
#plot
ggplot(dat,aes(x=Time,y=Abundance,color=Treatment))+
       geom_point()+geom_line(aes(y=Pred))
</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE0.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE0-490x368.png?w=450&ssl=1" alt="" class="wp-image-7105" data-recalc-dims="1"></a></p>
<h2>Second Example: Continuous Predictor</h2>
<p>We can also model the effect of continuous predictors on the parameters of the non-linear regression. For instance, we might assume that bacterial colonies grow faster in warmer temperature. In this case we could model the effect of temperature on growth rate (assuming that it is linear) and use the fitted growth rate to model the number of bacteria, all this within one model, pretty neat.<br>
To do so I will use another package: <a rel="nofollow" href="https://cran.r-project.org/web/packages/bbmle/index.html" target="_blank">bbmle</a> and its function <code>mle2</code>:</p>
<pre>
#load libraries
library(bbmle)
library(reshape2)
#parameter for simulation
K &lt;- 200
n0 &lt;- 5
#the gradient in temperature
Temp &lt;- ceiling(seq(0,20,length=10))
 

#simulate some data
mm &lt;- sapply(Temp,function(x){ 
             rpois(50,logF(time = 1:50,K = K,n0=n0,r = 0.05 + 0.01 * x))})
#sample data from a poisson distribution with lambda parameter equal
#to the expected value, note that we do not provide one value for the
#parameter r but a vector of values representing the effect of temperature
#on r
#some reshaping
datT &lt;- melt(mm)
names(datT) &lt;- c(&quot;Time&quot;,&quot;Temp&quot;,&quot;Abund&quot;)
datT$Temp &lt;- Temp[datT$Temp]

#fit the model
(mll &lt;- mle2(Abund~dpois(K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1))),
             data=datT,parameters = list(r~Temp),
             start=list(K=100,n0=10,r=0.05)))
<em>Call:
mle2(minuslogl = Abund ~ dpois(K * n0 * exp(r * Time)/(K + n0 * (exp(r * Time) - 1))), start = list(K = 100, n0 = 10, r = 0.05), data = datT, parameters = list(r ~ Temp))
Coefficients:
K n0 r.(Intercept) r.Temp 200.13176811 4.60546966 0.05195730 0.01016994 Log-likelihood: -1744.01 </em>
</pre>
<p>The first argument for <code>mle2</code> is either a function returning the negative log-likelihood or a formula of the form “response ~ some distribution(expected value)”. In our case we fit a poisson distribution with expected values coming form the logistic equation and its three parameters. The dependency between the growth rate and the temperature is given using the “parameter” argument. The basic model output gives us all the estimated parameters, of interest here are the “r.(Intercept)” which is the growth rate when the temperature is at 0 and “r.Temp” which is the increase in the growth rate for every increase of temperature by 1.</p>
<p>Again we can get the fitted values for plotting:</p>
<pre>
datT$pred &lt;- predict(mll)
ggplot(datT,aes(x=Time,y=Abund,color=factor(Temp)))+
       geom_point()+geom_line(aes(y=pred))
</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE1-490x368.png?w=450&ssl=1" alt="" class="wp-image-7104" data-recalc-dims="1"></a></p>
<h2>Third Example: More on Continuous Predictors</h2>
<p>The cool thing with <code>mle2</code> is that you can fit any models that you can imagine, as long as you are able to write down the log-likelihood functions. We will see this with an extension of the previous model. It is a common assumption in biology that species should have some optimum temperature, hence we can expect a bell-shape relation between population equilibrium and temperature. So we will add to the previous model a quadratic relation between the population equilibrium (“K”) and temperature.</p>
<pre>
#simulate some data
mm &lt;- sapply(Temp,function(x){
 rpois(50,logF(time = 1:50,K = 100+20*x-x**2,n0=n0,r = 0.05 + 0.01 * x))})
#note that this time both K and r are given as vectors to represent their
#relation with temperature

#some reshaping
datT &lt;- melt(mm)
names(datT) &lt;- c(&quot;Time&quot;,&quot;Temp&quot;,&quot;Abund&quot;)
datT$Temp &lt;- Temp[datT$Temp]

#the negative log-likelihood function
LL &lt;- function(k0,k1,k2,n0,r0,r1){
  K &lt;- k0 + k1*Temp + k2*Temp**2
  r &lt;- r0 + r1*Temp
  lbd &lt;- K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1))
  -sum(dpois(Abund,lbd,log=TRUE))
}

(mll2 &lt;- mle2(LL,data=datT,start=list(k0=100,k1=1,k2=-0.1,n0=10,r0=0.1,r1=0.1)))
<em>Call:
mle2(minuslogl = LL, start = list(k0 = 100, k1 = 1, k2 = -0.1, n0 = 10, r0 = 0.1, r1 = 0.1), data = datT)
Coefficients:
k0 k1 k2 n0 r0 r1 100.47610915 20.09677780 -1.00836974 4.84940342 0.04960558 0.01018838 Log-likelihood: -1700.91 </em>
</pre>
<p>This time I defined the negative log-likelihood function to fit the model. This function takes as argument the model parameter. In this case we have 3 parameters to describe the quadratic relation between K and temperature, 1 parameter for the constant n0 and 2 parameters to describe the linear relation between r and temperature. The important thing to understand is what the last line is doing. Basically the <code>dpois</code> call return the probabilities to get the observed abundance values (“Abund”) based on the expected mean value (“lbd”). The <code>LL</code> function will return a single negative log-likelihood value (the sum) and the job of <code>mle2</code> is to minimize this value. To do so <code>mle2</code> will try many different parameter combination each time calling the LL function and when it reaches the minimum it will return the parameter set used. This is Maximum Likelihood Estimation (MLE). As before the model output gives us the estimated parameters, and we can use these to plot the fitted regressions:</p>
<pre>
cc &lt;- coef(mll2)
#sorry for the awful coding, to get the model predicted values we need
#to code the relation by hand which is rather clumsy ...
datT$pred &lt;- melt(sapply(Temp,function(x) {logF(time=1:50,K=cc[1]+cc[2]*x+cc[3]*x**2,n0=cc[4],r=cc[5]+cc[6]*x))})[,3]
ggplot(datT,aes(x=Time,y=Abund,color=factor(Temp)))+
       geom_point()+geom_line(aes(y=pred))
</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE2.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2017/06/MLE2-490x368.png?w=450&ssl=1" alt="" class="wp-image-7103" data-recalc-dims="1"></a></p>
<h2>Conclusion</h2>
<p>I showed three way by which you can include the effect of predictors on the parameters of some non-linear regression: nlsList, mle2 with formula and mle2 with customized negative log-likelihood function. The first option is super easy and might be enough in many cases, while the last option though more complex gives you (almost) infinite freedom in the models that you can fit.</p>
<p>Happy non-linear modelling!</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/weather-forecast-with-regression-models-part-4/" title="Permanent Link toWeather forecast with regression models – part 4" target="_blank">Weather forecast with regression models – part 4</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/weather-forecast-with-regression-models-part-3/" title="Permanent Link toWeather forecast with regression models – part 3" target="_blank">Weather forecast with regression models – part 3</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/weather-forecast-with-regression-models-part-2/" title="Permanent Link toWeather forecast with regression models – part 2" target="_blank">Weather forecast with regression models – part 2</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/weather-forecast-with-regression-models-part-1/" title="Permanent Link toWeather forecast with regression models – part 1" target="_blank">Weather forecast with regression models – part 1</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/weighted-linear-support-vector-machine/" title="Permanent Link toWeighted Linear Support Vector Machine" target="_blank">Weighted Linear Support Vector Machine</a></li>
</ol></strong>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://datascienceplus.com/second-step-with-non-linear-regression-adding-predictors/"> R Programming – DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/H2iAwrG5tn4" height="1"></a></div>