I'm trying to figure out how to reproduce in Python some work that I've done in SAS. Using this dataset, where multicollinearity is a problem, I would like to perform principal component analysis in Python. I've looked at scikit-learn and statsmodels, but I'm uncertain how to take their output and convert it to the same results structure as SAS. For one thing, SAS appears to perform PCA on the correlation matrix when you use PROC PRINCOMP
, but most (all?) of the Python libraries appear to use SVD.
In the dataset, the first column is the response variable and the next 5 are predictive variables, called pred1-pred5.
In SAS, the general workflow is:
/* Get the PCs */
proc princomp data=indata out=pcdata;
var pred1 pred2 pred3 pred4 pred5;
run;
/* Standardize the response variable */
proc standard data=pcdata mean=0 std=1 out=pcdata2;
var response;
run;
/* Compare some models */
proc reg data=pcdata2;
Reg: model response = pred1 pred2 pred3 pred4 pred5 / vif;
PCa: model response = prin1-prin5 / vif;
PCfinal: model response = prin1 prin2 / vif;
run;
quit;
/* Use Proc PLS to to PCR Replacement - dropping pred5 */
/* This gets me my parameter estimates for the original data */
proc pls data=indata method=pcr nfac=2;
model response = pred1 pred2 pred3 pred4 / solution;
run;
quit;
I know that the last step only works because I'm only choosing PC1 and PC2, in order.
So, in Python, this is about as far as I've gotten:
import pandas as pd
import numpy as np
from sklearn.decomposition.pca import PCA
source = pd.read_csv('C:/sourcedata.csv')
# Create a pandas DataFrame object
frame = pd.DataFrame(source)
# Make sure we are working with the proper data -- drop the response variable
cols = [col for col in frame.columns if col not in ['response']]
frame2 = frame[cols]
pca = PCA(n_components=5)
pca.fit(frame2)
The amount of variance that each PC explains?
print pca.explained_variance_ratio_
Out[190]:
array([ 9.99997603e-01, 2.01265023e-06, 2.70712663e-07,
1.11512302e-07, 2.40310191e-09])
What are these? Eigenvectors?
print pca.components_
Out[179]:
array([[ -4.32840645e-04, -7.18123771e-04, -9.99989955e-01,
-4.40303223e-03, -2.46115129e-05],
[ 1.00991662e-01, 8.75383248e-02, -4.46418880e-03,
9.89353169e-01, 5.74291257e-02],
[ -1.04223303e-02, 9.96159390e-01, -3.28435046e-04,
-8.68305757e-02, -4.26467920e-03],
[ -7.04377522e-03, 7.60168675e-04, -2.30933755e-04,
5.85966587e-02, -9.98256573e-01],
[ -9.94807648e-01, -1.55477793e-03, -1.30274879e-05,
1.00934650e-01, 1.29430210e-02]])
Are these the eigenvalues?
print pca.explained_variance_
Out[180]:
array([ 8.07640319e+09, 1.62550137e+04, 2.18638986e+03,
9.00620474e+02, 1.94084664e+01])
I'm at a bit of a loss on how to get from the Python results to actually performing Principal Component Regression (in Python). Do any of the Python libraries fill in the blanks to similarly to SAS?
Any tips are appreciated. I'm a little spoiled by the use of labels in SAS output and I'm not very familiar with pandas, numpy, scipy, or scikit-learn.
Edit:
So, it looks like sklearn won't operate directly on a pandas dataframe. Let's say that I convert it to a numpy array:
npa = frame2.values
npa
Here's what I get:
Out[52]:
array([[ 8.45300000e+01, 4.20730000e+02, 1.99443000e+05,
7.94000000e+02, 1.21100000e+02],
[ 2.12500000e+01, 2.73810000e+02, 4.31180000e+04,
1.69000000e+02, 6.28500000e+01],
[ 3.38200000e+01, 3.73870000e+02, 7.07290000e+04,
2.79000000e+02, 3.53600000e+01],
...,
[ 4.71400000e+01, 3.55890000e+02, 1.02597000e+05,
4.07000000e+02, 3.25200000e+01],
[ 1.40100000e+01, 3.04970000e+02, 2.56270000e+04,
9.90000000e+01, 7.32200000e+01],
[ 3.85300000e+01, 3.73230000e+02, 8.02200000e+04,
3.17000000e+02, 4.32300000e+01]])
If I then change the copy
parameter of sklearn's PCA to False,
it operates directly on the array, per the comment below.
pca = PCA(n_components=5,copy=False)
pca.fit(npa)
npa
Per the output, it looks like it replaced all of the values in npa
instead of appending anything to the array. What are the values in npa
now? The principal component scores for the original array?
Out[64]:
array([[ 3.91846649e+01, 5.32456568e+01, 1.03614689e+05,
4.06726542e+02, 6.59830027e+01],
[ -2.40953351e+01, -9.36743432e+01, -5.27103110e+04,
-2.18273458e+02, 7.73300268e+00],
[ -1.15253351e+01, 6.38565684e+00, -2.50993110e+04,
-1.08273458e+02, -1.97569973e+01],
...,
[ 1.79466488e+00, -1.15943432e+01, 6.76868901e+03,
1.97265416e+01, -2.25969973e+01],
[ -3.13353351e+01, -6.25143432e+01, -7.02013110e+04,
-2.88273458e+02, 1.81030027e+01],
[ -6.81533512e+00, 5.74565684e+00, -1.56083110e+04,
-7.02734584e+01, -1.18869973e+01]])
copy=False
, I get new values. Are those the principal component scores? $\endgroup$