4
$\begingroup$

I am trying to fit hierarchical mixture model (using ML and MCMC, but this shouldn't matter) where the linear predictor part contains 17 independent variables. These are habitat variables: for each habitat type I have one variable saying the proportions of the area in 100 m circle which belongs to that particular habitat type.

The thing is that these 17 predictor variables sum up to 1 (i.e. simplex).

Could this be a problem with 1) fitting the model 2) model selection 3) predictions? This is not exactly collinearity (there is no correlation coefficient over 0.4 or under -0.4), but the variables are linearly dependent - the each one could be derived from all the others. If there is too much of a certain habitat, there cannot be a lot of other habitat types.

EDIT: The correlogram is here (the number is correlation coefficient multiplied by 100 and rounded. Only significant p < 0.05 coefficients are displayed).

enter image description here

EDIT 2: please do not assume that the variables are correlated. They are slightly in my case, but in general the variables can be linearly dependent but with no correlation! Look at this artificialy generated example:

set.seed(1063)
x <- rmultinom(17, rep(1000, 17), rep(1/17, 17))
envV <- x/1000

(If you have different RNG, please download the generated matrix: http://pastebin.com/sK55w3Y2)

Now the columns of envV are linearly dependent, as they sum up to 1 (see apply(envV, 2, sum)), but they are not correlated. See:

cor.mtest <- function(mat, conf.level = 0.95){
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat <- lowCI.mat <- uppCI.mat <- matrix(NA, n, n)
    diag(p.mat) <- 0
    diag(lowCI.mat) <- diag(uppCI.mat) <- 1
    for(i in 1:(n-1)){
        for(j in (i+1):n){
            tmp <- cor.test(mat[,i], mat[,j], conf.level = conf.level)
            p.mat[i,j] <- p.mat[j,i] <- tmp$p.value
			lowCI.mat[i,j] <- lowCI.mat[j,i] <- tmp$conf.int[1]
            uppCI.mat[i,j] <- uppCI.mat[j,i] <- tmp$conf.int[2]
        }
    }
    return(list(p.mat, lowCI.mat, uppCI.mat))
}

cor1 <- cor.mtest(envV)

number_of_correlated_variables <- sum(cor1[[1]] < 0.05 & !diag(nrow(cor1[[1]])))
number_of_correlated_variables  # reports 0

EDIT 3: It is interesting and strange that the matrix can be actually inverted: if I do solve(t(as.matrix(envV)) %*% as.matrix(envV)) both on my real predictors and the artificial generated ones in EDIT 2, the inverse matrix will get computed with no error or warning (both with solve and ginv()).

Also:

> is.singular.matrix(t(as.matrix(envV)) %*% as.matrix(envV))
[1] FALSE
$\endgroup$
15
  • 4
    $\begingroup$ Curious: If so, then the matrix is not of full rank. If your software tells you it is, it is wrong and you should find better software! $\endgroup$ Dec 30, 2014 at 17:12
  • 1
    $\begingroup$ Curious: It is also the possibility you are misusing R. I guess you are confused. $\endgroup$ Dec 30, 2014 at 17:50
  • 1
    $\begingroup$ @Curious, you get "different" results with lm() and solve(), because you were supplyingh the design matrix without the intercept into solve(), while lm() adds it by default, as Khashaa noted $\endgroup$
    – Aksakal
    Dec 30, 2014 at 19:05
  • 1
    $\begingroup$ The sum sum to 1 thing isn't necessarily an issue. If this forum was open I could post a proof for linear independence in $\mathbf{X}$ (because you aren't including a constant). We do this kind of stuff in regression all the time with fixed effects. where we don't include an intercept and all the effects in the model, every row sums to one but the results of the regressions are still valid (and well published in academic research). $\endgroup$ Dec 30, 2014 at 19:54
  • 3
    $\begingroup$ There are a couple of problems which I do find concerning in your simulation. 1) you have only 17 observations for 17 variables which essentially means you have zero degrees of freedom. 2) From the way you set up your simulation it looks like all the observations for each variable sum to 1, which is different from saying that the sum of all variables for each observation equals one. The latter is more pertinent to the discussion of this thread. Which one is it? $\endgroup$ Dec 30, 2014 at 19:55

2 Answers 2

4
$\begingroup$

The problem is that you have only 16 predictors, not 17. You can take any 16 variables and compute the remaining. This is the case of perfect milticollinearity. The solution is to simply regress on any subset of 16 variables.

The perfect multicollinearity (also rank deficiency) is the problem of identification in OLS. It's basically a technicality where you'd need to invert the design matrix $X'X$, and it's impossible when one of your variables is the linear combination of others.

Look up this wiki page for the definitions of the perfect multicollinearity condition: $\lambda_0 + \lambda_1 X_{1i} + \lambda_2 X_{2i} + \cdots + \lambda_k X_{ki} = 0$, which is what you have in your data with $k=17$, $\lambda_i=1, i\in [1,17]$ and $\lambda_0=-1$

UPDATE 3 Thanks to @Khashaa, he noticed that you did not use the intercept in your test, that's why it seems to pass. Here's the correct test code:

> x=runif(100*16)/17;
> x1=matrix(x,100,16)
> envV=cbind(x1,1-apply(x1,1,sum),rep(1,100))> apply(envV,1,sum)
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> hat=t(envV)%*%envV
> solve(hat)
Error in solve.default(hat) : 
  system is computationally singular: reciprocal condition number = 4.17876e-18

UPDATE 2 The reason why you don't see the correlations high is because you have a linear relationship with 17 variables, so pair-wise correlations don't have to be high. Consider this $x_{17}=1-\sum_{i=1}^{16}x_i, x_i\sim\mathcal{N}(0,1)$, if you look at the pair-wise correlations, you get $Corr[x_{17},x_i]=\frac{1}{\sqrt{16}}$, a relatively low number. So, in this case you're not going to notice multicollinearity by looking at pair-wise correlations.

$\endgroup$
17
  • $\begingroup$ I am afraid you are messing up linear dependence with correlation (collinearity). The fact that the variables sum up to 1 doesn't mean they are correlated - see the example I've added in the edit. $\endgroup$
    – Tomas
    Dec 30, 2014 at 16:16
  • $\begingroup$ @Curious, why don't you check the definition of multicollinearity and convince yourself that OP's data fits it? $\endgroup$
    – Aksakal
    Dec 30, 2014 at 16:19
  • $\begingroup$ 1) I am the OP :) 2) First sentence of the definition of Multicollinearity at wikipedia says "statistical phenomenon in which two or more predictor variables in a multiple regression model are highly correlated" - which they are not, see my EDIT 2 $\endgroup$
    – Tomas
    Dec 30, 2014 at 16:23
  • 1
    $\begingroup$ @Curious I think @Aksakal assumed that intercept is included in your predictors. But I see no trace of it in envV. $\endgroup$
    – Khashaa
    Dec 30, 2014 at 17:08
  • 2
    $\begingroup$ @Heisenberg lm automatically includes intercept thus creates multicollinearity. Sometimes people mistakenly include all the dummies in the regression, but fortunately R is smart enough to spot the error. $\endgroup$
    – Khashaa
    Dec 30, 2014 at 17:44
1
$\begingroup$

It could be an issue. Thinking of the design matrix, $X$ the last column could be expressed as $1 - $ sum of other columns. That means $X^TX$ would not be invertible.

Two options that might help, the first of which is probably going to be easier:

  • Drop a column. As its value is perfectly defined by the remaining data, you lose no information.
  • Use (Bayesian equivalent of) lasso regression, which will effectively choose one to drop for you.

Low entries in the correlogram are unlikely to signify we can relax. The proportions must be negatively related in some sense as for one to get larger, another must get smaller.

$\endgroup$
9
  • $\begingroup$ It is interesting and strange that the matrix can be actually inverted: if I do solve(t(as.matrix(envV)) %*% as.matrix(envV)) both on my real predictors and the artificial generated ones in EDIT 2, the inverse matrix will get computed with no error or warning (both with solve and ginv()) $\endgroup$
    – Tomas
    Dec 30, 2014 at 16:28
  • $\begingroup$ @Curious, show us the code for artificially generated matrix test. $\endgroup$
    – Aksakal
    Dec 30, 2014 at 16:30
  • $\begingroup$ @Aksakal see EDIT 2. $\endgroup$
    – Tomas
    Dec 30, 2014 at 16:33
  • $\begingroup$ +1 thanks for the lasso tip! As for dropping the column, I don't think this is a good solution. If the habitat I drop is important for abundance of the species (response variable), I will loose its prediction power and significance of this coefficient (since all the other ones will have no important meaning for the species). $\endgroup$
    – Tomas
    Dec 30, 2014 at 16:59
  • 1
    $\begingroup$ @Curious, by dropping the column you will NOT lose anything. All the information from that column is contained in the remaining 16 columns. $\endgroup$
    – Aksakal
    Dec 30, 2014 at 18:37

Not the answer you're looking for? Browse other questions tagged or ask your own question.