"The PCA analysis resulted in 4 factors above 1, but I stuck to the conceptual model of 3."

i had a look at your attached file, this 3-factor solution is not useful. you should check your 4-factor solution, maybe it makes more sense.

So, first things first.

1.Why PCA and not Principal Axis Factoring?

2. Varimax usually gives us the false impression that it makes the results of the factor loadings easier to read, but it is important that it isn't the sole purpose of varimax. Varimax presumes that you have an orthogonal structure, where presumably, factors are not correlated between each other. So it is expected that 1. Your items are supposed to load in a single factor and have loadings close to zero in other factors after rotating; 2. There's a tendency where usually while some of the items are going to have high factor loadings, other items are going to display a very small factor loading at that same factor as it "pulls" towards the extremes.

3. Your results are showing that you have several cross-loadings there and you commented that your variance shows that there are 4 factors with eigenvalues above 1. I understand that according to guidelines, we have to interpret factors according to some theory (known or temporary) and that residual non-explainable factors may show up. Still, beware of the amount of variance that you are leaving behind.

4. I strongly recommend that you verify the results with an oblique rotation and a more robust estimation method other than PCA. The reason why is that hardly social variables are going to behave like an orthogonal structure. It happens, yes, but your factor loadings are suggesting that this might not be the case of your data.

Another thing is that apparently your 4th dimension might be relevant, because squishing the variance in 3 factors might be provoking the cross-loadings in your model. It might not solve your problem magically, of course, but would be fun if come of the cross-loadings either get smaller or get totally reduced... It is just intuition.

Sorry to say, but it's not a good sign.  It is likely that either your questionnaire or the data is problematic.

I suggest you to check the wordings in your questionnaire to ensure that there is no misunderstanding.

Do the participants in pilot test / pre-test think the items/questions in the questionnaire are measuring the right thing (i.e. different aspects of the underlying variable you want to measure)?

Thanks for all the input.

Based on the feedback and my sparse knowledge, I did the following. Please, see if that's a correct way.

I ran parallel analysis and I got this result:

> fa.parallel(myData)

"Parallel analysis suggests that the number of factors = 7 and the number of components = 4"

So I tested both 4 and 7 factors, using 'oblimin' rotation now, and Principal Axis as factoring method:

> solution print(solution$loadings, cutoff=.3)

The result was much cleaner (see attached files).

Still, 4 items have similar (and low) cross-loadings in 4-factor solution and 3 items in 7-factor.

Ex.:

Item ID1: - 0.310 0.360 -

Item RU3: - 0.386 -0.337 -

What should I do with these items?

PS.: there was a lot of 'NA' in the data, I don´t know if this is problematic, since I imputed these missing values with the respective medians.

So that's actually good news. Some of the options you took simplified how variance is organized in the latent structure. So, still, one thing I'd recommend is for you to look at the proportion of N/As you have. If it is too high, this might be generating some noise.

It is good to see that the 4 factors structure has now 3 cross-loadings and the 7 factors has only 1 cross-loading. So I would try to fine tune the oblimin rotation, which has some options and I would look curiously also at a Maximum Likelihood estimation option. With ML estimation, the loadings are going to reduce even further I think.

With the current small cross-loadings, I would also think if there is a theoretical explanation for that to happen. Comparing the 4 and 7 factors solution, I prefer simple solutions. I don't use R since some time, but I know there are two different informations given by parallel if I'm not wrong which is the number of factors and how many of them should be actually retained.

Good progress there!

i think there are some problems in the data . whatever our educational research has such problems  .

in my opinion you can depend on high value   for example 

Item Pl2 -> will belong to  Factor 1: 0.542

Item AP2 -> will belong to Factor 1: 0.446;

Item TE2 -> Factor 1: will belong to  factore 2