Mutual Information larger than smaller one of both entropies? Lot's of questions and good answers on mutual information e.g. here out there and I think I get the concept which is also nicely explained on wikipedia. But I'm nervous that R's infotheo package is outputting some funny results... or I'm not getting it after all ;-)
In all self contained examples with synthetic data sets I tried I could not provoke this behaviour so I will not be able to provide a fully self contained example but after running
MI.matrix <- mutinformation(df)

I am getting the following results:
> MI.matrix[1,1]
[1] 0.4276829

So the entropy of the 1st variable (lets call it X) is ~0.43
> MI.matrix[440,440]
[1] 0.24493

So the entropy of the 440th variable (lets call it Y) is ~0.24
> MI.matrix[1,440]
[1] 0.4276829

So I(X;Y)=H(X) with H(X)>H(Y)?! If you refer to the image on the wiki article, how can the shared portion be larger than the smaller one of the entropies?
Note: I have also seen examples where H(X)=0, H(Y)=y and I(X;Y)=y where to my understanding I(X;Y) should be 0 as well.
Edit: after some more debugging I found the following:

*

*variable pairs which behave funny are one with lots of NAs

*While infotheo throws no errors with NAs and I could not find a hint in the documentation that NAs are not allowed I get a very different and correct result when removing NAs.

Below you can see a table with three variables, one version with all NAs included (df.temp) which has 1282 cases and one with complete cases only (df.temp.completecases) which is down to 15 cases.

Computing the MI matrix for both of them yields the following:

So now the I(X;Y) values seem to make more sense but the H(X) values along the diagonal obviously suffer from kicking out all X rows where Y is NA. I think I'll end up having to to compute the whole thing twice: once the entropies on individual columns and a 2nd run on pairwise complete columns to get the correct mutual information.
Guess my question has changed to: anyone's got an idea how to efficiently limit the pairwise MI computations to complete cases for a given pair only? Something similar to use=pairwisecomplete.obs in the cor function.
 A: You are correct, mutual information between X and Y cannot be larger than either the entropy of X or the entropy of Y. Ever. To debug where you've gone astray: instead of passing in a 2D matrix, try passing in pairs of vectors as separate X and Y arguments, and see if you can make sense of the results you get that way. Also pay attention to the method argument; empirical entropy estimates can be very biased unless you have a huge dataset; but it looks like corrections for this are supported.
A: With a little help from a former colleague I think I have a solution incl. a self contained example:
v1 <- c(1,2,3,4,5,NA,NA,NA,NA,NA)
v2 <- c(1,NA,3,NA,5,NA,7,NA,9,NA)
v3 <- c(NA,2,3,NA,NA,6,7,NA,7,NA)
v4 <- c(NA,NA,NA,NA,NA,6,7,8,9,10)
df <- cbind.data.frame(v1,v2,v3,v4)


ColPairMap<-function(df){
t <- data.frame(matrix(ncol = ncol(df), nrow = ncol(df)))
colnames(t) <- colnames(df)
rownames(t) <- colnames(df)
for (j in 1:ncol(df)) {
               for (i in 1:ncol(df)) {
                                c(1:ncol(df))
                                if (nrow(df[complete.cases(df[,c(i,j)]),])>0) {
                                    t[j,i] <- natstobits(mutinformation(df[complete.cases(df[,c(i,j)]),j], df[complete.cases(df[,c(i,j)]),i]))
                                } else {
                                    t[j,i] <- 0
                                }
               }
}
return(t)
}


ColPairMap(df)

However, performance is not exactly great with large data sets so if anyone's got an idea how to improve on that...
A: I've playing a bit more with this and the solution described above can still produce I(x;y)>max(H(x),H(y)). The reason is:

*

*in that nested loop for i!=j we can get fewer complete.cases than for i==j

*as a result I(x;y) can sometimes be computed on fewer rows than H(x) or H(y) which can result in I(x;y)>max(H(x),H(y))

To illustrate a small example:
ColPairMap<-function(df){ #function for pairwise removal of incomplete cases and subsequent MI computation per pair of variables
    t <- matrix(0, ncol = ncol(df), nrow = ncol(df), dimnames = list(colnames(df), colnames(df)))
    df <- as.matrix(df)
    for (j in 1:ncol(df)) {
        for (i in j:ncol(df)) {
            compl_cases <- complete.cases(df[, c(i, j)])
            if (sum(compl_cases) > 0) {
                t[j,i] <- natstobits(mutinformation(df[compl_cases, j], df[compl_cases, i], method="emp")) #using method="mm" can cause >1 after normalization!
            }
        }
    }
    lt <- lower.tri(t)
    t[lt] <- t[lt] + t(t)[lt]
    return(t)
}

# Easy, no missing data
V1<-c(6,4,1,5,21,5,21,21,21,9)
V2<-c(1,10,6,9,21,9,21,21,21,20)
test1<-cbind.data.frame(V1,V2)
ColPairMap(test1)

# Missing data in both variables but the same rows!
# Consequently the entropy for a single variable is being computed on the exact same rows as the mutual information!
# As a result, mutual information can never be larger than any of the two individual entropies!
V1<-c(6,4,1,NA,21,5,21,21,NA,9)
V2<-c(1,10,6,NA,21,9,21,21,NA,20)
test3<-cbind.data.frame(V1,V2)
ColPairMap(test3)

# Missing data in both variables but in different rows!
# This leads to the entropy for a single variable being computed on a larger set of rows than the mutual information!
# As a result, mutual information can be larger than both individual entropies!
V1<-c(6,4,1,5,NA,5,21,21,21,9)
V2<-c(1,10,6,NA,21,9,21,21,NA,20)
test2<-cbind.data.frame(V1,V2)
ColPairMap(test2)

How would you go about this? Accept this typically small imperfection or compute individual H(x) for every different y and normalize directly inside the function or ...?
A: To complete this topic I would like to share my final version which combines the idea provided by BioLiason with parallelized loops using foreach. On small data sets this is in fact slower but on large data sets this is significantly faster, especially when running on some proper H/W.
library(foreach)
library(parallel)
library(doParallel)
library(infotheo)

n <- 250 #creates an nXn matrix, the larger the more compute time is required
df <- (discretize(matrix(rnorm(n*n,n,n/10),ncol=n)))

## pairwise complete mutual information via nested for loop ##
start_for <- Sys.time()
ColPairMap<-function(df){
t <- data.frame(matrix(ncol = ncol(df), nrow = ncol(df)))
colnames(t) <- colnames(df)
rownames(t) <- colnames(df)
for (j in 1:ncol(df)) {
               for (i in 1:ncol(df)) {
                                c(1:ncol(df))
                                if (nrow(df[complete.cases(df[,c(i,j)]),])>0) {
                                    t[j,i] <- natstobits(mutinformation(df[complete.cases(df[,c(i,j)]),j], df[complete.cases(df[,c(i,j)]),i]))
                                } else {
                                    t[j,i] <- 0
                                }
               }
}
return(t)
}
ColPairMap(df)
end_for <- Sys.time()
end_for-start_for


## pairwise complete mutual information via nested foreach loop ##
start_foreach <- Sys.time()
ncl <- max(2,floor(detectCores()*0.75)) #number of cores
clst <- makeCluster(n=ncl,type="TERR") #create cluster
e <- new.env() #new environment to export libraries to cores
e$libs <- .libPaths()
clusterExport(clst, "libs", envir=e) #export required packages to all cores
clusterEvalQ(clst, .libPaths(libs)) #export required packages to all cores
clusterEvalQ(clst, { #export required packages to all cores
    library(infotheo)
})
registerDoParallel(cl = clst) #register cluster

t <- foreach (j=1:ncol(df), .combine="c") %:% #parallellized nested loop for computing normalized pairwise complete MI between all columns
        foreach (i=j:ncol(df), .combine="c", .packages="infotheo") %dopar% {
            combine="c"
            compl_cases <- complete.cases(df[,c(i,j)])
            if (sum(compl_cases) > 0) {
                natstobits(mutinformation(df[compl_cases,][,j], df[compl_cases,][,i]))
            } else {
                0
            }
        }

RCA_MI_Matrix <- matrix(0, ncol = ncol(df), nrow = ncol(df), dimnames = list(colnames(df), colnames(df))) #set-up empty matrix for MI values
RCA_MI_Matrix[lower.tri(RCA_MI_Matrix, diag=TRUE)] <- t #fill lower triangle with MI values from nested loop
RCA_MI_Matrix[upper.tri(RCA_MI_Matrix)] <- t(RCA_MI_Matrix)[upper.tri(RCA_MI_Matrix)] #mirror lower triangle of matrix into upper one
end_foreach <- Sys.time()
end_foreach-start_foreach

