How to systematically remove collinear variables (pandas columns) in Python? Thus far, I have removed collinear variables as part of the data preparation process by looking at correlation tables and eliminating variables that are above a certain threshold. Is there a more accepted way of doing this? Additionally, I am aware that only looking at correlation amongst 2 variables at a time is not ideal, measurements like VIF take into account potential correlation across several variables. How would one go about systematically choosing variable combinations that do not exhibit multicollinearity?
I have my data within a pandas data frame and am using sklearn's models. 
 A: You can try use below code:
from statsmodels.stats.outliers_influence import variance_inflation_factor

def calculate_vif_(X):

    '''X - pandas dataframe'''
    thresh = 5.0
    variables = range(X.shape[1])

    for i in np.arange(0, len(variables)):
        vif = [variance_inflation_factor(X[variables].values, ix) for ix in range(X[variables].shape[1])]
        print(vif)
        maxloc = vif.index(max(vif))
        if max(vif) > thresh:
            print('dropping \'' + X[variables].columns[maxloc] + '\' at index: ' + str(maxloc))
            del variables[maxloc]

    print('Remaining variables:')
    print(X.columns[variables])
    return X

It works, but I don't like the performance of that approach
A: I tried SpanishBoy's answer and found serval errors when running it for a data-frame. Here is a debugged solution. 
from statsmodels.stats.outliers_influence import variance_inflation_factor    

def calculate_vif_(X, thresh=100):
cols = X.columns
variables = np.arange(X.shape[1])
dropped=True
while dropped:
    dropped=False
    c = X[cols[variables]].values
    vif = [variance_inflation_factor(c, ix) for ix in np.arange(c.shape[1])]

    maxloc = vif.index(max(vif))
    if max(vif) > thresh:
        print('dropping \'' + X[cols[variables]].columns[maxloc] + '\' at index: ' + str(maxloc))
        variables = np.delete(variables, maxloc)
        dropped=True

print('Remaining variables:')
print(X.columns[variables])
return X[cols[variables]]

I also had no issues with performance, but have not tested it extensively. 
A: Thanks SpanishBoy - It is a good piece of code.
@ilanman: This checks VIF values and then drops variables whose VIF is more than 5. By "performance", I think he means run time. The above code took me about 3 hours to run on about 300 variables, 5000 rows.
By the way, I have modified it to remove some extra loops. Also, i've made it a bit cleaner and return the dataframe with reduced variables. This version reduced my run time by half!
My code is below- Hope it helps.
from statsmodels.stats.outliers_influence import variance_inflation_factor    

def calculate_vif_(X, thresh=5.0):
    X = X.assign(const=1)  # faster than add_constant from statsmodels
    variables = list(range(X.shape[1]))
    dropped = True
    while dropped:
        dropped = False
        vif = [variance_inflation_factor(X.iloc[:, variables].values, ix)
               for ix in range(X.iloc[:, variables].shape[1])]
        vif = vif[:-1]  # don't let the constant be removed in the loop.
        maxloc = vif.index(max(vif))
        if max(vif) > thresh:
            print('dropping \'' + X.iloc[:, variables].columns[maxloc] +
                  '\' at index: ' + str(maxloc))
            del variables[maxloc]
            dropped = True

    print('Remaining variables:')
    print(X.columns[variables[:-1]])
    return X.iloc[:, variables[:-1]]

