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?
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
y <- rnorm(b %*% x) wrong? can't execute it. Also, what does it mean?
$\endgroup$
Commented
Oct 27, 2016 at 21:33
do.call and mapply for an R example showing something so basic is a really over-engineered way of doing things.
$\endgroup$
Commented
Oct 27, 2016 at 21:34
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)$.
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:
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.
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
