Load the tumor growth data set from the url http://benzekry.perso.math.cnrs.fr/DONNEES/data_exam.csv into a dataframe
df = read.csv('http://benzekry.perso.math.cnrs.fr/DONNEES/data_exam.csv', sep=";")
Load the time vector in a variable time
time = df$Time
Load the volume data in a variable V
V = df$Volume
Linear least-squares¶
We will first assume a constant error model (i.e. $\sigma_j=\sigma,\, \forall j$) and an exponential structural model: $$ V\left(t; \left(V_0, \alpha \right)\right) = V_0 e^{\alpha t}. $$ We can transform the problem so that it reduces to a linear regression.
Define a variable y as the log of V
y=log(V)
Using the formula seen in class, build the least-squares matrix $M$ for fitting y
Solve the system corresponding to the linear regression
ones = rep(1, length(time))
M = cbind(ones, time)
theta = solve(t(M)%*%M, t(M)%*%y)
print(theta)
Plot the regression line together with the data
time_plot = seq(time[1], tail(time, 1), by=0.1)
y_plot = theta[1] + theta[2]*time_plot
plot(time, y)
lines(time_plot, y_plot)
Considering that the number of injected cells is $10^6$ cells, which corresponds to $V_0 = 1$ mm$^3$, and looking at the fit, what do you conclude about the validity of the exponential model?
The estimate of $\sigma^2$ is given by $$ s^2 = \frac{1}{n-2}\sum_{j=1}^n\left(y_j - M\hat{\theta}\right)^2 $$ with $\hat{\theta}$ the vector of optimal parameters just found and $n$ is the number of time points.
If $$residuals = y-M\hat{\theta}$$ is the vector of residuals, then $s^2$ can be computed as $$ s^2 = \frac{1}{n-2}residuals^T\cdot residuals $$ with $residuals^T$ the tranpose of the vector $residuals$. Using these considerations, compute $s^2$.
residuals = y - M%*%theta
s2 = 1/(length(time)-2)*t(residuals)%*%residuals # MLE estimator of error model variance
s2 = s2[1,1]
Deduce the estimation of the covariance matrix of the parameter estimates, given by $$ s^2 \left(M^T M\right)^{-1} $$
cov_est = s2*solve(t(M)%*%M)
cov_est
Compute the standard errors on the parameter estimates.
se = sqrt(diag(cov_est))
se
Use the built-in ordinary linear least-squares function lm() to verify the results
lm_res = lm(y ~ 1 + time)
summary(lm_res)