######################################################
#Lecture 7 Script
#Factor analysis and interaction contrasts
#Anastasia Ejova 8/12/2014
#R101/PSY532
######################################################

#FACTOR ANALYSIS 
#Example 1
#Analysis steps 1 and 3 described in slides; SSFA dataset created based on these steps.

#Example of what was done for Step 2
hist(SS$PostListPractice)
hist(SS$PostListDeserve)
  #Scores are definitely not normally distributed.

#Example of what was done for Step 3
SSPracNoZeros <- subset(SS, subset = PostListPractice != 0)
library(lsr)
correlate(SSPracNoZeros, corr.method = "spearman", test = TRUE, p.adjust.method = "none")
#Look at CORRELATIONS column for PostListPractice: Still some significant correlations with 
#other PostList variables. Thus, kept this variable in the analysis.

#Did this for all variables, including PostListBet (not in the SS database we have
#been using here), and PostListCHANCE
SSChanceNoZeros <- subset(SS, subset = PostListCHANCE != 10)
correlate(SSChanceNoZeros, corr.method = "spearman", test = TRUE, p.adjust.method = "none")
#Look at CORRELATIONS column for PostListCHANCE: No significant correlations with 
#other PostList variables. Thus, excluded from factor analysis.

#SSFA data set thus contains the scores for 13 of the original 15 items in PostList
#intended to measure illusion of control.

#Step 4
library(psych)
KMO(SSFA) #overall MSA close to 1: good
cortest.bartlett(SSFA) #significant: good

#Step 5 - Scree plots and parallel analysis - psych package
#still using psych package, but need nFactors package for parallel analysis
#install.packages("nFactors")
library(nFactors)
fa.parallel(SSFA, 
            n.obs = NULL, #we are using a data frame and not a correlation matrix
            fm="pa", #use principal axis factoring (rather than maximum likelihood) for the factor analysis
            fa="both", #show both, the principal components and PAF scree plots
            main = "Parallel Analysis Scree Plots", #title of plot
            n.iter=20,error.bars=FALSE,SMC=FALSE,ylabel=NULL,show.legend=TRUE, #defaults
            sim=FALSE) #do not show an additional analysis involving simulations

#A parallel analysis in the nFactors package
parAnC <- nScree(SSFA, model = "components")
summary(parAn) #this version suggests 2 components, as does the Pc scree plot (xs in the 
#plot above)

parAnF <- nScree(SSFA, model = "factors")
summary(parAnF) #this version for factor analysis suggests no multi-factor structure

#So, will examine the two-factor and three-factor solutions

#Step 6 - Principal components 
#With orthogonal rotation
#install.packages("GPArotation")
library(GPArotation) #necessary for getting rotated solutions

#2 factors
PCorth2 <- principal(SSFA, nfactors = 2, rotate = "varimax", scores = FALSE)
print.psych(PCorth2, #show loadings from analysis
            cut = 0.3, #do not show loadings less than 0.3
            sort = TRUE) 
#interpretable in line with the literature on
#the illusions of natural and supernatural control
#Factor 1 represents natural IoC, while factor 2 represents supernatural IoC

#3 factors
PCorth3 <- principal(SSFA, nfactors = 3, rotate = "varimax", scores = FALSE)
print.psych(PCorth3, cut = 0.3, sort = TRUE) #components 1 and 3 show some cross-loading

#With oblique rotation
#2 factors
PCobl2 <- principal(SSFA, nfactors = 2, rotate = "oblimin", scores = FALSE)
print.psych(PCobl2, cut = 0.3, sort = TRUE) 
#Notice the "component correlation"
#matrix at the end

#3 factors
PCobl3 <- principal(SSFA, nfactors = 3, rotate = "oblimin", scores = FALSE)
print.psych(PCobl3, cut = 0.3, sort = TRUE)

#Step 7 - Principal axis factoring
#With orthogonal rotation
#You can infer the code for yourself based on the code below.
#Not run here because the factors are clearly correlated given the previous analysis,
#and known theory

#With oblique rotation
#2 factors
PAFobl2 <- fa(SSFA, fm = "pa", #specifying PAF as the fitting procedure
              nfactors = 2, rotate = "oblimin")
print.psych(PAFobl2, cut = .3, sort = TRUE)
#By removing cut = .3, you get the numbers for the equations on Slide 5
#Again, a very interpretable solution where
#Factor 1 represents natural IoC, while factor 2 represents supernatural IoC

#3 factors
PAFobl3 <- fa(SSFA, fm = "pa", #specifying PAF as the fitting procedure
              nfactors = 3, rotate = "oblimin")
print.psych(PAFobl3, cut = 0.3, sort = TRUE)
#The RMSEA index for this model is lower, suggesting a better fit. However, the 
#two-factor model is more interpretable, and the RMSEA of that is still within an
#acceptable range. RMSR same across models.
#http://www.psych.umass.edu/uploads/people/79/Fit_Indices.pdf (RMSEA info)
#http://www.ats.ucla.edu/stat/seminars/muthen_08/part2.pdf (RMSR info)

#Step 8 - Deleting items with multiple loadings
#Missed because (a) our solution is interpretable, and (b) we have very few items,
#and wish to keep as much information in the analysis as possible

#Step 9 - Interpretation
#Done throughout

#Step 10 - Factor scores (for final solution: 2-factor PAF with oblique rotation)
PAFobl2 <- fa(SSFA, fm = "pa", #specifying PAF as the fitting procedure
              nfactors = 2, rotate = "oblimin", 
              scores = "regression") #default method, suitable for correlated factors
              #according to Reading 3
              #Other options available (e.g., "Bartlett")
              #see the factor.scores function in the psych package as an extension

PAFobl2$scores
#Can add this vector to our SS database for use as the outcome variable (instead of 
#scale averages)

SS$NatIoCScore <- PAFobl2$scores[,1]
SS$SupIoCScore <- PAFobl2$scores[,2]

#Step 11 - Cronbach's alpha
#Create data frame of variables loading on first factor (natural IoC)
SSNatIoC <- data.frame(SS$PostListPractice, SS$PostListGoalieMoves, SS$PostListSkill, SS$PostListComputerGames,
                       SS$PostListKickDirn, SS$PostListSoccerKnowl, SS$PostListStrategy, SS$PostListPlayers)

library(psych) #load package again, since a number of other packages have a package "alpha"
alpha(SSNatIoC) # Calculate alpha: std.alpha = .88, quite good

#Create data frame of variables loading on first factor (natural IoC)
SSSupIoC <- data.frame(SS$PostListDeserve, SS$PostListLuckyPerson, SS$PostListLuckyPlay, 
                       SS$PostListLuckyMoments, SS$PostListLuck)
alpha(SSSupIoC) #std.alpha = .85, which is again quite good

#--
#Example 2: Based on Reading 6 ("Exploratory factor analysis" section and "Appendix B")
#We are aiming to determine the factor structure of 62 items about erroneous gambling-related beliefs.
#100 statements described in Reading 6 were presented to participants, who expressed their agreement/
#disagreement on a 7-point scale, in which -3 represented "strongly disagree", 0 represented "neither
#agree nor disagree", and 3 represented "strongly agree". An example of a statement is
#"My gambling wins are evidence that I have skill and knowledge related to gambling". We hypothesised
#that all the statements express some form of illusion of control, and should therefore exhibit
#a two-factor structure, with factors corresponding to the "illusion of natural control" and the "illusion
#of supernatural control". Only 62 of the original 100 items are being analysed, partly because of the
#removal of some items in Step 3 of the factor analysis (see below and "Removal of weak items" section in 
#the paper). Also, statements relating to "cycles of luck" are removed because a confirmatory factor
#analysis (see paper) showed these items to be related to both factors. 

#This factor analysis is not covered in the lecture slides because it was conducted as a follow-up to 
#a confirmatory factor analysis, which we have not covered. However, the supplementary code for removing
#weak items in Step 3 could prove useful to you in your own analyses, and you will see that the factor
#analysis steps 4 to 10 are the same with this increased number of items (as compared to only 13 items
#in Examnple 1).

#Step 2
hist(GRBFull$No1GRB_ChoicesAffect)
hist(GRBFull$No2GRB_GamblingWins)
#Variables still not normally distributed enough to warran ML fitting procedure
#So, for factor analysis, will use PAF

#Step 3
#See supplementary code below. The code creates mat.3 and mat.8, which record, for each item, the number of
#other items with which it had a correlation of 0.3 or less (mat.3) or a correlation of 0.8 or more (mat.8). 
#In creating the GRB database for the factor analysis, items that correlated with 60 or more items at a 
#magnitude of less than 0.3 were removed. Items correlating 0.8 or more with each other would have been
#inspected, and one of the pair would have been removed. However, there were no such items (mat.8 consists
#entirely of 0s). See Reading 6 for which other items were excluded to create the 62 variables in the 
#GRB data set used in this analysis.

#Step 4
#From this point forward, using the dataset GRB with 62 retained items.
library(psych)
KMO(GRB) #close to 1: good
cortest.bartlett(GRB) #significant: good

#Step 5
fa.parallel(GRB, 
            n.obs = NULL, #we are using a data frame and not a correlation matrix
            fm="pa", #use principal axis factoring (rather than maximum likelihood) for the factor analysis
            fa="both", #show both, the principal components and PAF scree plots
            main = "Parallel Analysis Scree Plots", #title of plot
            n.iter=20,error.bars=FALSE,SMC=FALSE,ylabel=NULL,show.legend=TRUE, #defaults
            sim=FALSE) #do not show an additional analysis involving simulations
#3 factors recommended for Principal Components and 5 for Principal Axis Factoring

#A parallel analysis in the nFactors package
parAn2C <- nScree(GRB, model = "components")
summary(parAn2C)
#9 components recommended

parAn2F <- nScree(GRB, model = "factors")
summary(parAn2F)
#5 factors recommended

#Step 6 - Try 2-factor (hypothesised), 3-factor, and 9-factor solutions
#Will do all with oblique rotation
Ex2PCobl2 <- principal(GRB, nfactors = 2, rotate = "oblimin", scores = FALSE)
print.psych(Ex2PCobl2, cut = 0.3, sort = TRUE) 
#Component 2 does contain most references to luck, lucky charms and so on (see ppaer)

Ex2PCobl3 <- principal(GRB, nfactors = 3, rotate = "oblimin", scores = FALSE)
print.psych(Ex2PCobl3, cut = 0.3, sort = TRUE)
#Components 2 and 3 are clearly representing the illuisons of supernatural and natural control, 
#respectively. Component 1 seems to relate to the gambler's fallacy (see paper for what this is).

Ex2PCobl9 <- principal(GRB, nfactors = 9, rotate = "oblimin", scores = FALSE)
#no convergence: so, not a good fit using PCA, even though recommended by parallel analysis

#Step 7 - Try 2-factor (hypothesised), 3-factor (comparison with PCA) and 5-factor solutions
#Will do all with oblique rotation
Ex2PAFobl2 <- fa(GRB, fm = "pa", nfactors = 2, rotate = "oblimin")
print.psych(Ex2PAFobl2, cut = .3, sort = TRUE)
#Even clearer split into items about the "natural" and "supernatural", compared to the PCA solution.

Ex2PAFobl3 <- fa(GRB, fm = "pa", nfactors = 3, rotate = "oblimin")
print.psych(Ex2PAFobl3, cut = .3, sort = TRUE)
#As for PCA: F1 - Gambler's fallacy, F2 - Illusion of supernatural control, F3 - Illusion of natural control
#Note that the RMSEA is identical to the two-factor model

Ex2PAFobl5 <- fa(GRB, fm = "pa", nfactors = 5, rotate = "oblimin")
print.psych(Ex2PAFobl5, cut = .3, sort = TRUE)
#Not very interpretable from a theoretical point of view

#Step 8
#Assuming that my chosen final soltuion is Ex2PAFobl2 (2-factor PAF), the cross-loading/weak items are
#No62GRB_BadVibes, No76GRB_LuckRubsOff, No36GRB_LearningExperience, No33GRB_Q12_ShowMe, 
#No4GRB_KnowlAndSkill_MostPeople, No59GRB_LuckyTechnique, and No75GRB_GenerallyLucky 
#We could develop a questionnaire where we do not include these items. However, running the analyses
#again without them is not necessary for the purposes of the paper.

#Step 9
#Done throughout, and commented on in the paper.

#Step 10
#Factor scores not necessary, since no further analyses were planned using natural and supernatural IoC
#as outcome measures. However, we do need them for the 

#Step 11
#Example for Natural IoC factor (the code is somewhat tedious)
GRBNatIoC <- data.frame(GRB$No40GRB_FreeSpins, GRB$No53GRB_WinComingSoon, 
                        GRB$No21GRB_MaintainStreak,
                        GRB$No12GRB_ChangingResponses, GRB$No46GRB_MaintainStrategy, GRB$No34GRB_Q3_Familiar,
                        GRB$No41GRB_Q29, GRB$No18GRB_MachineTime, GRB$No6GRB_SameBet, GRB$No47GRB_Q24_LongerLosing,
                        GRB$No99GRB_Q15_ContestOfWills, GRB$No19GRB_PersonTypeWins, GRB$No38GRB_IdentifyingPattern,
                        GRB$No74GRB_Q22_BornLucky, GRB$No9GRB_Q25_SystemForPlay, GRB$No2GRB_GamblingWins, 
                        GRB$No30GRB_PokerMachChallenge, GRB$No56GRB_RandomPatterns, GRB$No10GRB_WinOnce,
                        GRB$No37GRB_FollowStrategies, GRB$No22GRB_ChangeFlow, GRB$No29GRB_PreviousBets, 
                        GRB$No42GRB_NearWin, GRB$No93GRB_MakeMoney, GRB$No96GRB_WinBack, GRB$No35GRB_Q20_WinningLotteryNos,
                        GRB$No32GRB_Q16_SecretsToLearn, GRB$No1GRB_ChoicesAffect, GRB$No3GRB_KnowlAndSkill_Likelihood,
                        GRB$No100GRB_Q10_Competition, GRB$No45GRB_DontWantToMiss, GRB$No39GRB_TimeOfDay,
                        GRB$No49GRB_BoundToFollow, GRB$No31GRB_Q21_Sportsperson, GRB$No60GRB_Q11_Fate,
                        GRB$No62GRB_BadVibes, GRB$No97GRB_Dedication, GRB$No76GRB_LuckRubsOff, GRB$No36GRB_LearningExperience,
                        GRB$No13GRB_BeatSystem, GRB$No33GRB_Q12_ShowMe, GRB$No11GRB_ControlOverPredicting,
                        GRB$No4GRB_KnowlAndSkill_MostPeople)
library(psych)
alpha(GRBNatIoC) #.96 - again, good
#If you get error messages, check that you haven't included a variable in the data frame twice!

######################################################

#INTERACTION CONTRASTS
#Two-way interactions
#The final model from Assignment 4
Sp1 <- Spending[(Spending$Spend < 10000),] #removing outliers
library(MASS)
glmnba <- glm.nb(Spend ~ Online*YoungMale*IncomeCat, data = Sp1)
summary(glmnba)

library(car)
Anova(glmnba, type = 2)
#Significant Online:YoungMale interaction. We will analyse this interaction
#using the specialised package, phia.

#install.packages("phia")
library(phia)

intMeans <- interactionMeans(glmnba, factors = c("Online", "YoungMale"))
intMeans #values are obtained from the slope coefficients (Estimates in summary)
#We can see the very high expenditure among young males playing online (690 Czech crowns)

plot(intMeans, multiple = TRUE) #for a visual representation, we can plot the above mean values
#Shows the main effects as well
#multiple can be set to FALSE if there are too many panels

#Simple effects analysis
#e.g., Testing the effect of "Online" (shifting from No to Yes) across different levels of 
#YoungMale (No, Yes) in a post-hoc manner
testInteractions(glmnba, fixed="YoungMale", across="Online")
#So, the shift has a significant effect among young males but not everyone else.

#Interaction contrasts
#notice change of "fixed" argument to "pairwise"
#Better option than simple effects analysis (above) which tends to conflate main effects
#and interaction effects. See slide for further explanation.
testInteractions(glmnba, pairwise ="YoungMale", across="Online")

#--

#More on interacation contrasts (see slides for additional information)
#Test of Hypothesis 2 in SS database
#Create categories of prior beliefs
SS$CatPreDBC_IOC <- cut(x = SS$PreDBC_IOC, breaks = 2, labels = c("Low", "High"))
summary(SS$CatPreDBC_IOC)

glmnbb <- glm.nb(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType *PreDBC_Total, 
                 data = SS)
summary(glmnbb)
Anova(glmnbb, Type = 2)

intMeansb <- interactionMeans(glmnbb, 
                              factors = c("SeqCond", "PostHowManySingleCaptionType"))
intMeansb
#The Ascending condition is the only one in which memory for wins is unaffected by wording.

levels(SS$SeqCond) <- c("Desc", "U", "Asc", "Flat") #So that the hyphen in "U-shaped" isn't confused for
#a minus

custom.contr <- contrastCoefficients(
  SeqCond ~ Asc - (Desc + U + Flat)/3, #Contrasting Ascending with all the others across levels of CaptionType
  SeqCond ~ Asc - Desc, 
  data=SS, normalize=TRUE)

custom.contr
names(custom.contr$SeqCond) <- c("Asc vs. all", "Asc vs. Desc") #giving the columns of the 
#contrast matrix more understandable names
custom.contr

testInteractions(glmnbb,custom=custom.contr)
#Can see from the results that the pattern of wording (CaptionType) effects is different in the Ascending
#condition

######################################################
#Supplementary code
#Working with the original data set for Example 2: GRBFull (contains 100 items)

store.3 <- -1 #equal -1 if this is the first participant
store.8 <- -1

#Record no of problematic corrs for all participants (creating vectors store.3 and store.8)
#and record, in a separate text file, details on which PARTICULAR correlations were <.3 
#and >.8 for each participant

for (i in 1:100) {
  
  count.3 <- 0
  count.8 <- 0
  
  store.3.id <- -1
  store.8.id <- -1
  
  cormat <- cor(x = GRBFull[i], y = GRBFull[,-i], method = "spearman")
  
  for (j in 1:99) {
    if (cormat[j] < .3) {
      count.3 <- count.3 + 1
      if (store.3.id == -1) {
        store.3.id <- j
      } else {
        store.3.id <- c(store.3.id, j)
      }
    }
  }
  
  for (k in 1:99) {
    if (cormat[k] > .8) {
      count.8 <- count.8 + 1
      if (store.8.id == -1) {
        store.8.id <- k
      } else {
        store.8.id <- c(store.8.id, k)
      }
    }
  }
  
  if (store.3 == -1) {
    store.3 <- count.3
  } else {
    store.3 <- c(store.3, count.3)
  }
  
  if (store.8 == -1) {
    store.8 <- count.8
  } else {
    store.8 <- c(store.8, count.8)
  }
  
  #Add to a text file identifying the corrs >=.3 (i.e., valid) and >=.8 
  #for each participant
  
  Nos1to99 <- c(1:99)
  
  if (count.8 != 0 & count.3 != 0) {
    store.valid <- Nos1to99[-c(store.3.id, store.8.id)]
  }
  
  if (count.8 == 0 & count.3 != 0) {
    store.valid <- Nos1to99[-c(store.3.id)]
  }
  
  if (count.8 != 0 & count.3 == 0) {
    store.valid <- Nos1to99[-c(store.8.id)]
  }
  
  if (count.8 == 0 & count.3 == 0) {
    store.valid <- Nos1to99
  }
  
  P.id <- i
  
  write(P.id, file = "CorrsID.txt", append = TRUE)
  cat("\nID of corrs >= .3 but <= .8\n", file = "CorrsID.txt", append = TRUE)
  write(store.valid, file = "CorrsID.txt", append = TRUE)
  cat("\n", file = "CorrsID.txt", append = TRUE)
  
  cat("\nID of corrs > .8\n", file = "CorrsID.txt", append = TRUE)
  write(store.8.id, file = "CorrsID.txt", append = TRUE)
  cat("\n\n", file = "CorrsID.txt", append = TRUE)
  
}

#Make vectors store.3 and store.8 into matrices
ItemID <- 1
rownames <- "Q 1"
for (q in 1:99) {
  ItemID <- ItemID + 1
  row <- paste("Q", as.character(ItemID))
  rownames <- c(rownames, row)
}

mat.3 <- as.matrix(store.3, nrow = 100, ncol = 1, dimnames = list(rownames, "no < .3"))
mat.8 <- as.matrix(store.8, nrow = 100, ncol = 1, dimnames = list(rownames, "no > .8"))

#Remove unnecessary variables
rm(count.3, count.8, store.3.id, store.8.id, i, j, k, ItemID, rownames, row, Nos1to99, P.id, q, store.valid)
