I have tried to reproduce some research (using PCA) from SPSS in R. In my experience, principal() function from package psych was the only function that came close (or if my memory serves me right, dead on) to match the output. To match the same results as in SPSS, I had to use parameter principal(..., rotate = "varimax"). I have seen papers talk about how they did PCA, but based on the output of SPSS and use of rotation, it sounds more like Factor analysis.

Question: Is PCA, even after rotation (using varimax), still PCA? I was under the impression that this might be in fact Factor analysis... In case it's not, what details am I missing?

This question is largely about definitions of PCA/FA, so opinions might differ. My opinion is that PCA+varimax should not be called either PCA or FA, bur rather explicitly referred to e.g. as "varimax-rotated PCA".

I should add that this is quite a confusing topic. In this answer I want to explain what a rotation actually is; this will require some mathematics. A casual reader can skip directly to the illustration. Only then we can discuss whether PCA+rotation should or should not be called "PCA".

One reference is Jolliffe's book "Principal Component Analysis", section 11.1 "Rotation of Principal Components", but I find it could be clearer.



Let $\mathbf X$ be a $n \times p$ data matrix which we assume is centered. PCA amounts (see my answer here) to a singular-value decomposition: $\mathbf X=\mathbf{USV}^\top$. There are two equivalent but complimentary views on this decomposition: a more PCA-style "projection" view and a more FA-style "latent variables" view.

According to the PCA-style view, we found a bunch of orthogonal directions $\mathbf V$ (these are eigenvectors of the covariance matrix, also called "principal directions" or "axes"), and "principal components" $\mathbf{US}$ (also called principal component "scores") are projections of the data on these directions. Principal components are uncorrelated, the first one has maximally possible variance, etc. We can write: $$\mathbf X = \mathbf{US}\cdot \mathbf V^\top = \text{Scores} \cdot \text{Principal directions}.$$

According to the FA-style view, we found some uncorrelated unit-variance "latent factors" that give rise to the observed variables via "loadings". Indeed, $\widetilde{\mathbf U}=\sqrt{n-1}\mathbf{U}$ are standardized principal components (uncorrelated and with unit variance), and if we define loadings as $\mathbf L = \mathbf{VS}/\sqrt{n-1}$, then  $$\mathbf X= \sqrt{n-1}\mathbf{U}\cdot (\mathbf{VS}/\sqrt{n-1})^\top =\widetilde{\mathbf U}\cdot \mathbf L^\top = \text{Standardized scores} \cdot \text{Loadings}.$$ (Note that $\mathbf{S}^\top=\mathbf{S}$.) Both views are equivalent. Note that loadings are eigenvectors scaled by the respective eigenvalues ($\mathbf{S}/\sqrt{n-1}$ are eigenvalues of the covariance matrix).

(I should add in brackets that PCA$\ne$FA; FA explicitly aims at finding latent factors that are linearly mapped to the observed variables via loadings; it is more flexible than PCA and  yields different loadings. That is why I prefer to call the above "FA-style view on PCA" and not FA, even though some people take it to be one of FA methods.)  

Now, what does a rotation do? E.g. an orthogonal rotation, such as varimax. First, it considers only $k<p$ components, i.e.: $$\mathbf X \approx \mathbf U_k \mathbf S_k \mathbf V_k^\top = \widetilde{\mathbf U}_k \mathbf L^\top_k.$$ Then it takes a square orthogonal $k \times k$ matrix $\mathbf T$, and plugs $\mathbf T\mathbf T^\top=\mathbf I$ into this decomposition: $$\mathbf X \approx \mathbf U_k \mathbf S_k \mathbf V_k^\top = \mathbf U_k \mathbf T \mathbf T^\top \mathbf S_k \mathbf V_k^\top = \widetilde{\mathbf U}_\mathrm{rot} \mathbf L^\top_\mathrm{rot},$$ where rotated loadings are given by $\mathbf L_\mathrm{rot} = \mathbf L_k \mathbf T$, and rotated standardized scores are given by $\widetilde{\mathbf U}_\mathrm{rot} = \widetilde{\mathbf U}_k \mathbf T$. (The purpose of this is to find $\mathbf T$ such that $\mathbf L_\mathrm{rot}$ became as close to being sparse as possible, to facilitate its interpretation.)

Note that what is rotated are: (1) standardized scores, (2) loadings. But not the raw scores and not the principal directions! So the rotation happens in the latent space, not in the original space. This is absolutely crucial.

From the FA-style point of view, nothing much happened. (A) The latent factors are still uncorrelated and standardized. (B) They are still mapped to the observed variables via (rotated) loadings. (C) The amount of variance captured by each component/factor is given by the sum of squared values of the corresponding loadings column in $\mathbf L_\mathrm{rot}$. (D) Geometrically, loadings still span the same $k$-dimensional subspace in $\mathbb R^p$ (the subspace spanned by the first $k$ PCA eigenvectors). (E) The approximation to $\mathbf X$ and the reconstruction error did not change at all. (F) The covariance matrix is still approximated equally well:$$\boldsymbol \Sigma \approx \mathbf L_k\mathbf L_k^\top = \mathbf L_\mathrm{rot}\mathbf L_\mathrm{rot}^\top.$$

But the PCA-style point of view has practically collapsed. Rotated loadings do not correspond to orthogonal directions/axes in $\mathbb R^p$ anymore, i.e. columns of $\mathbf L_\mathrm{rot}$ are not orthogonal! Worse, if you [orthogonally] project the data onto the directions given by the rotated loadings, you will get correlated (!) projections and will not be able to recover the scores. [Instead, to compute the standardized scores after rotation, one needs to multiply the data matrix with the pseudo-inverse of loadings $\widetilde{\mathbf U}_\mathrm{rot} = \mathbf X (\mathbf L_\mathrm{rot}^+)^\top$. Alternatively, one can simply rotate the original standardized scores with the rotation matrix: $\widetilde{\mathbf U}_\mathrm{rot} = \widetilde{\mathbf U} \mathbf T$.] Also, the rotated components do not successively capture the maximal amount of variance: the variance gets redistributed among the components (even though all $k$ rotated components capture exactly as much variance as all $k$ original principal components).

Here is an illustration. The data is a 2D ellipse stretched along the main diagonal. First principal direction is the main diagonal, the second one is orthogonal to it. PCA loading vectors (eigenvectors scaled by the eigenvalues) are shown in red -- pointing in both directions and also stretched by a constant factor for visibility. Then I applied an orthogonal rotation by $30^\circ$ to the loadings. Resulting loading vectors are shown in magenta. Note how they are not orthogonal (!).



An FA-style intuition here is as follows: imagine a "latent space" where points fill a small circle (come from a 2D Gaussian with unit variances). These distribution of points is then stretched along the PCA loadings (red) to become the data ellipse that we see on this figure. However, the same distribution of points can be rotated and then stretched along the rotated PCA loadings (magenta) to become the same data ellipse.

[To actually see that an orthogonal rotation of loadings is a rotation, one needs to look at a PCA biplot; there the vectors/rays corresponding to original variables will simply rotate.] 



Let us summarize. After an orthogonal rotation (such as varimax), the "rotated-principal" axes are not orthogonal, and orthogonal projections on them do not make sense. So one should rather drop this whole axes/projections point of view. It would be weird to still call it PCA (which is all about projections with maximal variance etc.).

From FA-style point of view, we simply rotated our (standardized and uncorrelated) latent factors, which is a valid operation. There are no "projections" in FA; instead, latent factors generate the observed variables via loadings. This logic is still preserved. However, we started with principal components, which are not actually factors (as PCA is not the same as FA). So it would be weird to call it FA as well.

Instead of debating whether one "should" rather call it PCA or FA, I would suggest to be meticulous in specifying the exact used procedure: "PCA followed by a varimax rotation".



Postscriptum. It is possible to consider an alternative rotation procedure, where $\mathbf{TT}^\top$ is inserted between $\mathbf{US}$ and $\mathbf V^\top$. This would rotate raw scores and eigenvectors (instead of standardized scores and loadings). The biggest problem with this approach is that after such a "rotation", scores will not be uncorrelated anymore, which is pretty fatal for PCA. One can do it, but it is not how rotations are usually being understood and applied.
