# Creating new function in R equivalent to the lm() for finding coefficients b0 and b1

I have a task to estimate the model to obtain b0 and b1 on the generated data (yt, zt and wt) using my created function and I'm not sure, how to do it. It's written that I should not use the existing functions such as lm(). Could someone explain and help me to understand that?

I found that b1 can be found like that:

In my calculation b1=0.1646209 which is the same what lm() function shows:

Call: lm(formula = y ~ z) Coefficients: (Intercept) z
1.0064 0.1646

In that case b0 can be expressed like that

b0 = Yt - b1*zt - ut

But when I try to find b0 I get a lot of values instead of one, so probably this equation is not correct. Could you help me with finding b0 (intercept)?

Model: yt = b0 + b1*zt + ut

Variables:

n <- 100

v <- rnorm(n ,mean=0, sd=0.27)

eps <- rnorm(n ,mean=0, sd=0.05)

u <- -0.5*v + eps

w <- rnorm(n ,mean=0, sd=1)

z <- 0.2*w + v

y <- 1 + 0.5*z + u

• I think you need to give us more details about what you did and what you found so people have a chance of putting you in the right direction. Dec 4 '20 at 16:31
• This looks like self study. There a good PowerPoint presentations and online lectures talking about this. This question is been around for a very long time and it’s an important question. Dec 10 '20 at 12:51

You are mixing several things up when you solve the equation $$y_t = \beta_0 + \beta_1 z_t + u_t$$

by writing $$\beta_0 = y_t - \beta_1 z_t - u_t$$

and inserting an estimate for $$\beta_1$$.

• The equation is only correct with the true $$\beta_0$$, $$\beta_1$$, and $$u_t$$. But you fill in an estimate for $$\beta_1$$.

If you would compute it with the true $$\beta_1$$, that is if you compute Yt - 0.5*zt - ut instead of Yt - b1*zt - ut. Then you should get an array with all 1's (the value of b0).

• In regression you do not know the true values u_t that generated the model. Instead you are finding $$\beta_0$$ and $$\beta_1$$ such that

$$u_t = y_t - \beta_0 - \beta_1 z_t$$

are as small as possible. You find the $$\beta_0$$ and $$\beta_1$$ that minimize the sum of squared $$u_t$$ (or similarly maximize the likelihood of Gaussian distributed error). And now these $$u_t$$ are not true errors (which you normally do not know and you only have them in your computer simulation) but instead, they are residuals.

The code below demonstrates three methods to perform the linear regression:

• Method 1: Direct computation of the coefficients in simple linear regression.
• Method 2: Solving the matrix equation $$\beta = (X^TX)^{-1} X^T y$$ This is what the function lm does. The method 1 can only solve simple linear regression. This method 2 can also deal with larger model-matrices $$X$$.
• Method 3: Minimizing the objective function with an itterative algorithm. This is more common for non-linear models and used by the function nls.

R-code:

### data
set.seed(1)
n <- 100
v <- rnorm(n ,mean=0, sd=0.27)
eps <- rnorm(n ,mean=0, sd=0.05)
u <- -0.5*v + eps
w <- rnorm(n ,mean=0, sd=1)
z <- 0.2*w + v
y <- 1 + 0.5*z + u

plot(z,y)

####################################
######## solution with lm ##########
####################################

mod <- lm(y ~ z)

### result: Intercept = 0.9937834  Slope = 0.2064167
mod$coefficients #################################### ######## solution method 1 ######### ######## direct formulas ########### #################################### slope = cov(z,y)/var(z) intercept = mean(y) - slope * mean(z) ### result: Intercept = 0.9937834 Slope = 0.2064167 c(intercept, slope) #################################### ######## solution method 2 ######### ######## matrix formula ############ #################################### X <- cbind(rep(1,length(z)),z) beta <- solve(t(X) %*% X) %*% t(X) %*% y ### result: Intercept = 0.9937834 Slope = 0.2064167 beta #################################### ######## solution method 3 ######### ######## minimization of cost ###### #################################### f_to_minimize <- function(b) { ### estimated y y_est <- b[1] + b[2] * z ### sum of squared residuals SSR <- sum((y_est-y)^2) return(SSR) } ### optim minimizes the function ### you can use several algorithms begin_par <- c(0,0) mod_opt <- optim( par = begin_par , fn = f_to_minimize, method = "BFGS") ### result: Intercept = 0.9937814 Slope = 0.2064134 ### (these results are slighly different because this is an ### approximation of the estimate/optimum. You can tune the optim function to get closer) mod_opt$par
• Wow! It's such an amazing help, thank you so much for all explanations and given examples!!
– Gina
Dec 4 '20 at 15:03

For your model $$y_t = \beta_0 + \beta_1 z_t + u_t$$, the formulas for the OLS coefficients estimates $$\hat{\beta_0}$$ and $$\hat{\beta_1}$$ are:

• $$\hat{\beta_1} = \frac{\widehat{Cov}(y,z)}{\widehat{Var}(z)}$$, as you correctly pointed out, where $$\widehat{Cov}$$ and $$\widehat{Var}$$ are the sample covariances and variances respectively; in R: cov(y,z), var(y) where y and z are vectors;

• $$\hat{\beta_0}=\bar{y}-\hat\beta_1\cdot \bar{z}$$, where $$\bar y$$ and $$\bar z$$ are the sample averages; in R: mean(y) and mean(z) where y and z are vectors;

These values are found by minimizing $$\sum_{i=1}^{n}{(y_i-(\hat\beta_0+\hat \beta_1 z_i))^2}$$ with respect to $$\hat\beta_0$$ and $$\hat\beta_1$$ (i.e. set partial derivatives to 0 and solve).

.

Extra content: I would suggest using the matrix form though, since it generalizes to any number of coefficient estimates: $$\hat\beta_{OLS} = (X^TX)^{-1}X^Ty$$, where $$X$$ is the regressor matrix, i.e. first column is just 1s, second column is regressor 1 (here: vector z), third column regressor 2 etc. $$y$$ is the response vector. The rows are the observations. This operation would yield a vector $$\hat\beta_{OLS}$$ with the first element corresponding to $$\hat \beta_0$$, second element to $$\hat \beta_1$$ etc. In R: beta_OLS = solve(t(X)%*%X)%*%t(X)%*%y

• Thank you so much, very useful!
– Gina
Dec 4 '20 at 15:03
• @Gina I made a mistake (edited now): the columns of the regressor matrix $X$ are the regressors, not the rows. The rows are simply the observations from 1 to n. Dec 4 '20 at 15:42