r - How to get the standard error of the restricted least square regression in R

Question

I am estimating a restricted linear regression model

lm(TC~Q+PL+PK+PF)

under the linear restriction the coefficients of PL+PK+PF sum to one. I want both the regression coefficients and standard errors. How can I achieve this in R?

Answer 1

If the only constraint you have is that the estimates need to sum to 1 then you can construct this fairly easily by rewriting your model. Let's say we have 3 predictors X1, X2, X3 and we fit a model

y = B0 + B1*X1 + B2*X2 + B3*X3 + error

then notice that if we have the restriction that B1 + B2 + B3 = 1 that we could rewrite B3 = 1 - B1 - B2. In which case our model becomes

y = B0 + B1*X1 + B2*X2 + (1 - B1 - B2)*X3 + error
  = B0 + B1*(X1 - X3) + B2*(X2 - X3) + X3 + error
  = B0 + B1*P1 + B2*P2 + X3 + error

where we defined two new variables P1 = X1 - X3 and P2 = X2 - X3.

so if we can fit that model then we're in business. Notice that the estimated parameter for the variable P1 will be our estimate of B1, and our estimated parameter for the variable P2 will be our estimate of B2. We don't directly get an estimate of B3 but since B3 is just a function of B1 and B2 we'll be fine.

So let's generate some data and fit a model...

# Generate fake data and fix some parameters
set.seed(123412)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
b0 <- 2
b1 <- .2
b2 <- .5
b3 <- 1 - b1 - b2 # so that b1 + b2 + b3 = 1
e <- rnorm(n)
y <- b0 + b1*x1 + b2*x2 + b3*x3 + e


p1 <- x1 - x3
p2 <- x2 - x3

o <- lm(y ~ offset(x3) + p1 + p2)

Now we can easily get estimates for our parameters

b1hat <- coef(o)[2]
b2hat <- coef(o)[3]
b3hat <- 1 - b1hat - b2hat

If we want standard errors we can get SEs for the first two parameters from the output of summary (or by taking the square root of the diagonal of the output of vcov(o))

# Get standard errors for B1 and B2
summary(o)
sqrt(diag(vcov(o)))

# Get a standard error for B3 - which is just
# a linear combination of B1 and B2
D <- c(0, -1, -1)
b3hat.se <- sqrt( t(D) %*% vcov(o) %*% D)