# R lm, Could anyone give me an example of the misleading case on “prediction from a rank-deficient”?

Could anyone give an example of "prediction from a rank-deficient fit may be misleading"?

From following example, it seems fine.

fit1=lm(mpg~wt,data=mtcars)
fit2=lm(mpg~wt+I(2*wt)+I(3*wt),data=mtcars)
all(predict(fit1,mtcars)==predict(fit2,mtcars))


Why and when "misleading" will happen?

• "May" is an important qualifier. It all depends on why the fit is rank-deficient. Sometimes that occurs due to mathematical redundancies, such as using proportions of a whole for regressors: those proportions must sum to unity. There is nothing special about this situation vis a vis prediction. Othertimes rank deficiency might occur simply because you have not yet made an observation that fills out the missing dimensions. When that's the case, you basically have no information (apart from some measurement error) about the missing dimensions. It's easy to imagine what the problems might be. – whuber Oct 26 '16 at 14:52
• @whuber could you give me one numerical example that the misleading will happen? – hxd1011 Oct 29 '16 at 21:41

Your second fit contains collinear factors. I know because you defined them that way! No two variables in a linear model can be "recreated" value-by-value using any linear combination of other variables in the model. Otherwise you arrive at rank deficiency. R is smart enough to just drop redundant factors, so it simply sets your second and third coefficients to missing in the second model.

When R does this, it assumes that you did not do so in such an obvious way. For instance, if you adjusted for 100 factors, you may arrive at deficiency and not know it. For that reason, when you predict, R reminds you that the prediction model is so overadjusted, there is almost certainly a problem of overfitting that will not give you reliable predictions. There are methods to accommodate high dimensional predictions.

Take this example of using n=50 observations to fit a 50 feature prediction model versus a 20 feature prediction model.

set.seed(1)

p <- 50
n <- 200

b <- matrix(rnorm(p))
x <- matrix(rbinom(p*n, 1, .3), n, p)
y <- rnorm(n, sweep(x, b, FUN=*, MARGIN = 2))
x <- as.data.frame(x)
X <- data.frame('y'=y, 'x'=x)
train <- rep.int(1:0, c(50, 150))==1

## all factors are "important", but fit2 with the first 20 preds more
## generalizable
fit1 <- lm(y ~ ., subset=train, data=X)
fit2 <- lm(y ~ ., subset=train, data=X[, 1:21])

pred1 <- predict(fit1, newdata = X[!train, ])
pred2 <- predict(fit2, newdata = X[!train, ])

var(y[!train] - pred1)
var(y[!train] - pred2)


results in

> var(y[!train] - pred1)
[1] 15.07587
> var(y[!train] - pred2)
[1] 1.610317

• I know, and I intentionally did that to create co-linear problem. My question was, even this exist, it seems no "misleading" predictions. – hxd1011 Oct 27 '16 at 20:42
• @hxd1011 see my edit, I am trying to explain that if this was done unconscionably (for instance, high dimensional prediction), the warning is one of overfitting. This means your predictions in NEW data will suck. I will add an example. – AdamO Oct 27 '16 at 20:44
• is this line y <- rnorm(b %*% x) wrong? can't execute it. Also, what does it mean? – hxd1011 Oct 27 '16 at 21:33
• Nice answer, +1 but I have to R-comment this: do.call and mapply for an R example showing something so basic is a really over-engineered way of doing things. – usεr11852 Oct 27 '16 at 21:34
• @hxd1011 Did some out of order fudging that confused R. I never quite understand R's object casting. Anyway, the coding example is reproducible although I hope it was easy to understand what was suggested in the coding. – AdamO Oct 27 '16 at 22:20

Ok, say you want to predict someones weight based on their height and chest size. For simplicity assume we have a training set of 3 people: $$p_1$$, $$p_2$$ and $$p_3$$. Now our data is:
Chest size: $$p_1=80cm, p_2=90cm, p_3=100cm$$
Height: $$p_1=160cm, p_2=180cm, p_3=200cm$$
weight: $$p_1=59kg, p_2=68kg, p_3=76kg$$

Note: height equals exactly 2 times the chest size in each training case.
Now we run a linear regression on this. A result you might get is:
$$0.75 \cdot \text{Chest size} + 0 \cdot \text{height}$$
Another result you could get is:
$$0 \cdot \text{Chest size} + 0.375 \cdot \text{height}$$
More concerning, you also could get:
$$3000.75 \cdot \text{Chest size} - 1500 \cdot \text{height}$$

All these give exactly the same result, so the cost-function doesn't care which one of these it finds. However now suppose you apply this to a new person, $$p_4$$, this person has as measurements:
$$Chest size = 81cm$$
$$Height=160cm$$
This is basically $$p_1$$, but 1cm more chest size. However if we use that third set of parameters we get:
$$81\cdot 3000.75 - 160 \cdot 1500 = 3060.75kg$$
This is obviously ridiculous, for some with the same height as $$p_1$$ but 1 cm more chest size, we shouldn't predict that person to be 3000kg heavier. So we can indeed conclude that these coefficients are obviously wrong. But how would our regression algorithm know that these coefficients are wrong? They give the same output, have the same cost as very reasonable coefficients.

The underlying problem is that in our training-set, height and weight are linearly dependent. This allows for an infinite number of coefficient pairs, all with the same outpout in each training case. So the cost of each of that infinite number of coefficient pairs is the same. However as soon as an example comes along where those 2 parameters have a different relation, all those infinite pairs of coefficients suddenly have very different outputs on this new individual.
As we saw here, as soon as someone comes along that doesn't have a height:chest size relation of 2, we are lost. Then what appeared as very good parameter suddenly can predict crazy weights. However when we don't know our parameters are linearly dependent (it might be coincidence), we might not know that the possibility is there for such crazy scenario's. We might get weights like $$(3000.75, -1500)$$. Those weights are very misleading. That is what the function is warning for. Misleading weights like in this case $$(3000.75, -1500)$$.

• +1, but could you give me a R example for the incorrect result? – hxd1011 Nov 4 '16 at 13:44

Rank-deficiency means that we are unable to uniquely determine the relationship between response and covariates. The QR factorization with pivoting strategy used by lm ensures a unique solution in numerical point of view, but estimated coefficients do not match the truth, typically because it constrains some coefficients to 0. When using this set of coefficients for prediction, we can be away from truth.

In the following I will present a long simulation study with R.

In practice there could only be four scenarios:

1. model matrix is full-rank, predictor matrix is full-rank;
2. model matrix is full-rank, predictor matrix is rank-deficient;
3. model matrix is rank-deficient, predictor matrix is full-rank;
4. model matrix is rank-deficient, predictor matrix is rank-deficient.

In machine language terminology, the model matrix is the one associated with training dataset, and the predictor matrix is the one associated with test dataset. The following toy function generates two datasets (each with n data) for a regression model overall three numeric covariates and an intercept. model.rank1.defect and predictor.rank1.defect specify whether we want those matrices to be rank-1 deficient. The reason for using relatively few covariates are commented. The way to achieve deficiency is also commented. Note that the true model is generated on the complete dataset (of 2 * n data) which has full rank. A noise-to-signal ratio 0.1 is used.

sim <- function (n = 1000, model.rank1.defect = FALSE, predictor.rank1.defect = FALSE) {

## since we only impose rank-1 deficiency, we'd better try small number of parameters
## so that the degree of deficiency is relatively high
## We guess that higher degree of deficiency is,
## the easier we are going to spot "misleading" result in scatter plot
p <- 3

#############################
## similate a model matrix ##
#############################

## generate a full rank model matrix
Xm <- matrix(runif(n * p), n, p)
## make it rank-1 deficient if required
if (model.rank1.defect) {
# take a random column, replace it by the sum of remaining two columns
ind <- sample(p, 1)
Xm[, ind] <- rowSums(Xm[, -ind])
}

#################################
## similate a predictor matrix ##
#################################

## generate a full rank predictor matrix
Xp <- matrix(runif(n * p), n, p)
## make it rank-1 deficient if required
if (predictor.rank1.defect) {
# take a random column, replace it by the sum of remaining two columns
ind <- sample(p, 1)
Xp[, ind] <- rowSums(Xp[, -ind])
}

#########################
## assume a true model ##
#########################

beta <- rnorm(p)  ## true coefficients for covariates
y.true <- 0.5 + rbind(Xm, Xp) %*% beta  ## intercept = 0.5
y <- y.true + rnorm(2 * n, 0, sqrt(0.1 * var(y.true)))  ## noise to signal ratio = 0.1

###########################################
## separate training and testing dataset ##
###########################################

train <- data.frame(y = y[1:n], X = I(Xm), y.true = y.true[1:n])
test <- data.frame(y = y[-(1:n)], X = I(Xp), y.true = y.true[-(1:n)])

###########################################
## return datasets and true coefficients ##
###########################################

list(train = train, test = test, coef = c(0.5, beta))
}


Basically we want to see how good / bad the prediction could be (compared with the truth) in four cases. The following function produces scattered plot for four cases, and note that good prediction means that the dots are scattered along y = x line.

inspect <- function (seed = 0, n = 1000) {

set.seed(seed)

case1 <- sim(n, model.rank1.defect = FALSE, predictor.rank1.defect = FALSE)
case2 <- sim(n, model.rank1.defect = FALSE, predictor.rank1.defect = TRUE)
case3 <- sim(n, model.rank1.defect = TRUE, predictor.rank1.defect = FALSE)
case4 <- sim(n, model.rank1.defect = TRUE, predictor.rank1.defect = TRUE)

par(mfrow = c(2, 2))

fit1 <- lm(y ~ X, data = case1$train) plot(case1$test$y.true, predict(fit1, case1$test), main = "case1", xlab = "true", ylab = "predicted")

fit2 <- lm(y ~ X, data = case2$train) plot(case2$test$y.true, predict(fit2, case2$test), main = "case2", xlab = "true", ylab = "predicted")

fit3 <- lm(y ~ X, data = case3$train) plot(case3$test$y.true, predict(fit3, case3$test), main = "case3", xlab = "true", ylab = "predicted")

fit4 <- lm(y ~ X, data = case4$train) plot(case4$test$y.true, predict(fit4, case4$test), main = "case4", xlab = "true", ylab = "predicted")

list(case1 = rbind(truth = case1$coef, estimated = coef(fit1)), case2 = rbind(truth = case2$coef, estimated = coef(fit2)),
case3 = rbind(truth = case3$coef, estimated = coef(fit3)), case4 = rbind(truth = case4$coef, estimated = coef(fit4)))
}


It is generally a good idea to set a reasonably big n, say 1000, for clear visualization. You are recommended to try several random seeds, but I will just show the result for 0, 1, 2.

inspect(seed = 0, n = 1000)
#$case1 # (Intercept) X1 X2 X3 #truth 0.5000000 -0.5380638 0.4965946 -1.441966 #estimated 0.5086424 -0.5182216 0.4801427 -1.462437 # #$case2
#          (Intercept)       X1        X2         X3
#truth        0.500000 1.683990 -1.544324 -0.1908871
#estimated    0.485759 1.690604 -1.529844 -0.1870299
#
#$case3 # (Intercept) X1 X2 X3 #truth 0.5000000 -0.8292477 0.04317512 0.4441213 #estimated 0.5144345 -1.2640980 0.47286008 NA # #$case4
#          (Intercept)         X1         X2        X3
#truth       0.5000000 -1.3502605 -0.2520004 0.9186201
#estimated   0.4857474 -0.4213176 -1.1679039        NA


inspect(seed = 1, n = 1000)
#$case1 # (Intercept) X1 X2 X3 #truth 0.500000 0.7391149 0.3866087 1.296397 #estimated 0.494829 0.7564111 0.3780248 1.302040 # #$case2
#          (Intercept)        X1       X2         X3
#truth       0.5000000 0.1490125 1.572552 -0.4313906
#estimated   0.4926926 0.1454275 1.575065 -0.4292703
#
#$case3 # (Intercept) X1 X2 X3 #truth 0.5000000 0.2994154 -0.4803364 2.456836 #estimated 0.5033469 2.7558496 -2.9404693 NA # #$case4
#          (Intercept)        X1         X2        X3
#truth       0.5000000 0.5066593  0.7627762 -2.231709
#estimated   0.4883394 2.7323985 -1.4522909        NA


inspect(seed = 2, n = 1000)
#$case1 # (Intercept) X1 X2 X3 #truth 0.5000000 -2.138476 -0.9007918 0.1499151 #estimated 0.4996853 -2.147597 -0.8861852 0.1555451 # #$case2
#          (Intercept)       X1         X2        X3
#truth       0.5000000 1.415801 -0.3443290 0.1640968
#estimated   0.4943133 1.424356 -0.3759823 0.2018281
#
#$case3 # (Intercept) X1 X2 X3 #truth 0.5000000 0.5502895 -0.05330089 -1.024089 #estimated 0.4945481 -0.4709953 0.96391344 NA # #$case4
#          (Intercept)        X1         X2        X3
#truth       0.5000000 -1.001211 -1.8184362 0.9240334
#estimated   0.5221117 -1.903361 -0.9189928        NA


Note that for case 1 and case 2, the quality of the prediction appears steadily good. But for case 3 and case 4, it can be good or bad.

• for random seed 0, prediction is good, but may be less accurate;
• for random seed 2 and seed 3, the prediction is terrible, and can even be wrong as the dots are scattered parallel to y = -x.

The mismatch between true and estimated coefficients in case 3 and case 4 is the root cause (NA coefficient is just 0 coefficient with 0 standard error).

A rank-deficient fit happens because your $x$-variables (features) are linearly dependent. This essentially means that you have "more" $x$-variables than observations. This can happen for two reasons. Either you have too few observations, or too many $x$-variables.

As in your example, having too many $x$-variables is harmless when it comes to prediction, although it may cause you to get the wrong regression coefficients. But having too few observations, while it gives you the same error message, will be a much bigger problem, because it means that R will not be able to find the regression coefficients, even when the data are exactly fitted by a linear model.

Here is an example:

# the model is y = V1 + V2 + V3 + V4

X <- matrix(rnorm(8*4), ncol=4)
y <- rowSums(X)
dat <- data.frame(cbind(X, y))

# case 1: not enough observations; everything goes wrong

model1 <- lm(y ~., data=dat[6:8,])
pred1 <- predict(model1, dat[1:5,])
mse1 <- sum((y[1:5] - pred1)^2)/5

# case 2: just right!

model2 <- lm(y ~., data=dat[1:5,])
pred2 <- predict(model2, dat[6:8,])
mse2 <- sum((y[6:8] - pred2)^2)/3

# case 3: too many features (like in your example) is still fine for
# prediction, even though you get the same error message as case 1.

model3 <- lm(y ~. + I(2*V1), data=dat[1:5,])
pred3 <- predict(model3, dat[6:8,])
mse3 <- sum((y[6:8] - pred2)^2)/3

• Having variables that are a perfect linear combination of other variables in the fitted data will not alter prediction as long as those same variables are still the same perfect linear combination of the other variables in the new data to predict on. If there are not, then any number is a valid prediction from the OLS/MLE estimates. – Cliff AB Oct 27 '16 at 23:32
• I think the statement " This essentially means that you have "more" x-variables than observations." is not correct. I can have rank deficiency on 1000 observations, and 2 identical features. – hxd1011 Oct 28 '16 at 13:46