Using principal component analysis (PCA) for feature selection in regression I have a dataset $D$ made of $m$ samples and $n$ features with $n \gg m$.
For each sample I have a score $s$ which I would like to be able to predict.
As the number of features is very high (compared to the number of samples),
I have performed PCA (Principal Component Analysis) and built a regression
model for $s$ based on the first 3 principal components (out of 45). I cross-validated the results and am happy with the performances of my model.
Now I would like to understand what is happening in terms of my original features. In particular I would like to understand which features are relevant in predicting $s$, and select them to build a model which does not involve PCA.
How can I perform this last step (if you have R based solutions it is better, but not essential)?
 A: PCA by its own doesn't do a feature selection, but it performs a transformation that will usually include all original variates: it is rare that PCA loadings (the coefficients for the transformation/projection) are 0. 
Thus, I'd say that PCA is more adapted to data where you expect the relevant information to be spread out over many/most of your variates. 
There are other regularization techniques which favor coefficients becoming 0, such as LASSO - which corresponds to selecting variates.
The Elements of Statistical Learning have a nice comparison of different regularization techniques.
There is also sparse PCA which combines the idea of maximum variance projection with favoring sparse solutions.
I may add that also sparse PLS exists. 
As you ask about R, several implementations of sparse PCA, sparse PLS, etc. exist. I cannot recommend any particular package as my data is usually of the type where you cannot expect a sparse solution - I therefore use normal PCA and PLS.

If you prefer to stay with the PCR you have, but would like to interpret the coefficients along the lines @SamLivingstone suggests:


*

*If your data is not variance scaled, looking at the coefficients is not enough, you need to take into account the size/range of the corresponding variate ("contribution").

*I'd resample the models (e.g. during a validation run) and compare the size of the coefficient/contribution of each variate to its spread. If you do that for the PCA, you need to take into account that PCA axes can flip (change sign). Also, you may want to consider models equal that are just rotated in the 3 dimensions you keep. 

*for a PCR, I'd look at the final PCR coefficients (those that transform the original data directly into the prediction of the dependent variable): that saves you all the hassle of aligning equivalent PCA projections. (Same for PLS)

