install.packages('stringr')
install.packages('lsr')
install.packages('stats')
install.packages(c("ggplot2", "ggpubr", "tidyverse", "broom", "AICcmodavg"))
install.packages("dplyr")
install.packages("ggpubr")
install.packages("chron")

# *** FOLLOWING CODE IS RUN ONLY ONCE. GO DOWN TO WHERE THE NEW FILES ARE CREATED ***

library(stringr)
library(lsr)
library(ggplot2)
library(ggpubr)
library(tidyverse)
library(broom)
library(AICcmodavg)
library(chron)

setwd('c:/util/articles')
FileNameExpG = 'IntroverExtraverFromDrFink.csv'
dataExpGrp <- read.csv(FileNameExpG)
head(dataExpGrp)
dim(dataExpGrp)
datExpGrp <- dataExpGrp[-c(21:145)] # first 3 columns are ID, chart name, Julian day and zodiac longitudes. remove speed, rt asc etc.
head(datExpGrp)
dim(datExpGrp)

FileNameScores = 'nScoRNAR.dat'
dataScores <- read.table(FileNameScores)
head(dataScores)
dim(dataScores)

dataScores$name = paste(dataScores$V2,dataScores$V3) # merge 2 parts of the name into the name
datScores <- dataScores[-c(1:3)] # remove "N" in 1st column and the 2 parts of the name
colnames(datScores)[1] = "Scores" # change column name from "V4" to "Scores"
head(datScores)

write.csv(datScores,"dataScores.txt")
write.csv(datExpGrp,'datPlanPos.txt')

# Now run WinMerge and manually remove items in dataScores.txt that are not in datPlanPos.txt and vice versa
# After manual editing, the files were renamed to dataScoresFinal.txt and datPlanPosFinal.txt

# ******************************************************************************************************

# *** THE NEW FILES ARE CREATED IN THE CODE ABOVE. RUN ONLY FOLLOWING CODE TO REPLICATE THE ANALYSIS ***

# ******************************************************************************************************

# THESE LINES ARE ABOVE TOO AND ARE HERE TOO FOR CONVENENCE WHEN REPLICATING THE ANALYSIS:
setwd('c:/util/articles')
library(stringr)
library(lsr)
library(ggplot2)
library(ggpubr)
library(tidyverse)
library(broom)
library(AICcmodavg)
library(dplyr)
library(chron)

datSco = read.csv("dataScoresFinal.txt")
datExpG = read.csv("datPlanPosFinal.txt")
dim(datSco)
dim(datExpG)  # number of rows in datSco and ExpG must match! 
head(datSco)
head(datExpG)
# Assumed also that the manual editng with WinMerge have the names of the charts in the same order. Perhaps this should be verified

datExpG$Name = str_trim(datExpG$Name)  # remove leading space. Note: datSco$name does not have a leading space.
datSco$name[1]
datExpG$Name[1]

if (nrow(datSco) == nrow(datExpG)) {
  AllOK = TRUE
} else {
  AllOK = FALSE
}

if (AllOK) {
  NamesMatch = TRUE  # assume true
  
  for (row in 1:nrow(datSco)) {
    if (datSco$name[row] != datExpG$Name[row]) {
      NamesMatch = FALSE
    }
  }
}

if ((AllOK) && (NamesMatch)) {
  ZodSignName = c("01","02","03","04","05","06","07","08","09","10","11","12")
  ZodSignsArry <- matrix(ncol=14,nrow=nrow(datExpG))
  datExpG$Scores = datSco$Scores  # Add Scores field to datExpG. Now datExpG has planet positions and scores
  
  # create array of zodiac sign of planets Sun - Pluto
  for (row in 1:nrow(datExpG)) {
    for (j in 5:18) {
      st = ZodSignName[as.integer(datExpG[row,j] / 30) + 1]  # "01" = Aries "02"=Taurus . . . "12"=Pisces
      ZodSignsArry[row,j-4] = st
    }
  }
  
  datExpG$SunSign = ZodSignsArry[,1]  # add Sun sign column to datExpG
  datExpG$MooSign = ZodSignsArry[,2]  # add Moo sign column to datExpG
  datExpG$MerSign = ZodSignsArry[,3]  # add Mer sign column to datExpG
  datExpG$VenSign = ZodSignsArry[,4]  # add Ven sign column to datExpG
  datExpG$MarSign = ZodSignsArry[,5]  # add Mar sign column to datExpG
  datExpG$JupSign = ZodSignsArry[,6]  # add Jup sign column to datExpG
  datExpG$SatSign = ZodSignsArry[,7]  # add Sat sign column to datExpG
  datExpG$UraSign = ZodSignsArry[,8]  # add Ura sign column to datExpG
  datExpG$NepSign = ZodSignsArry[,9]  # add Nep sign column to datExpG
  datExpG$PluSign = ZodSignsArry[,10] # add Plu sign column to datExpG
  datExpG$TNoSign = ZodSignsArry[,11] # add TNo sign column to datExpG
  datExpG$AscSign = ZodSignsArry[,12] # add Asc sign column to datExpG
  datExpG$MCSign  = ZodSignsArry[,13] # add  MC sign column to datExpG
  datExpG$VtxSign = ZodSignsArry[,14] # add Vtx sign column to datExpG
}

# create columns for elements (fire/earth/air/water)
# first add the columns to datExpG and then the element of the sign will be put in the columns:

nr = nrow(datExpG)
nrvector = rep("X",nr)
datExpG = cbind(datExpG,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector,nrvector)
# add column names to the columns added in the code above:
colnames(datExpG)[37] = "SunElem"
colnames(datExpG)[38] = "MooElem"
colnames(datExpG)[39] = "MerElem"
colnames(datExpG)[40] = "VenElem"
colnames(datExpG)[41] = "MarElem"
colnames(datExpG)[42] = "JupElem"
colnames(datExpG)[43] = "SatElem"
colnames(datExpG)[44] = "UraElem"
colnames(datExpG)[45] = "NepElem"
colnames(datExpG)[46] = "PluElem"
colnames(datExpG)[47] = "TNoElem"
colnames(datExpG)[48] = "AscElem"
colnames(datExpG)[49] = "MCElem"
colnames(datExpG)[50] = "VtxElem"

for (row in 1:nrow(datExpG)) {
  for (j in 23:36) {              # 14 column of zodiac sign positions of Sun-Pluto, TNo< Asc, MC and Vtx
    if ( (datExpG[row,j]=="01") |
         (datExpG[row,j]=="05") | 
         (datExpG[row,j]=="09") 
    ) { datExpG[row,j+14] = "F" }  # new row of Element. In this case Fire
    else if 
    ( (datExpG[row,j]=="02") |
      (datExpG[row,j]=="06") | 
      (datExpG[row,j]=="10") 
    ) { datExpG[row,j+14] = "E" }  # new row of Element. In this case Earth
    else if 
    ( (datExpG[row,j]=="03") |
      (datExpG[row,j]=="07") | 
      (datExpG[row,j]=="11") 
    ) { datExpG[row,j+14] = "A" }  # new row of Element. In this case Air
    else if 
    ( (datExpG[row,j]=="04") |
      (datExpG[row,j]=="08") | 
      (datExpG[row,j]=="12") 
    ) { datExpG[row,j+14] = "W" }  # new row of Element. In this case Water
  }
}

nrnumbervector = rep(0,nr)
datExpG = cbind(datExpG,nrnumbervector)
colnames(datExpG)[51]="AgeOfPerson"

for (row in 1:nrow(datExpG)) {
  iColonPos = unlist(gregexpr(':',datExpG$Name[row]))
  YearStr = substr(datExpG$Name[row],iColonPos-2,iColonPos-1)     # Year that extraversion-introversion test was taken
  SemesterStr = substr(datExpG$Name[row],iColonPos-3,iColonPos-3) # S=Spring F=Fall M=Summer semester that extraversion-introversion test was taken

  if (SemesterStr == "F")  {
    imonth = 11      # fall semester about Sept - Dec so Set to Nov 1, about the middle of the semester
    iday = 1
  } else if (SemesterStr == "S") {
    imonth = 3       # spring semester about Jan - Apr so Set to Mar 1, about the middle of the semester
    iday = 1
  } else if (SemesterStr == "M") {
    imonth = 8       # summer semester about May or Jun to Aug or Sep so Set to Aug 1, about the middle of the semester
    iday = 1
  } else {
    imonth = 7       # this line should never be reached because SemesterStr should be "F", "S", or "M" but just in case to prevent the R script from crashing or producing bizarre results
    iday = 1         # July 1 is about the middle of the year so date can be only 6 months off at the most.
  }

  iyear = strtoi(YearStr) + 1900
  datExpG$AgeOfPerson[row] = (julian(imonth,iday,iyear) - (datExpG$Jul_Day[row] - 2440588.0)) / 365.25 # age of person in years when the person took the extraversion-introversion test. Note: 2440588.0 = julian day on Jan 1 1970, which is the date that the R routine for julian day begins. The person's age is accurate within about 2 months because the time of taking the test is known within about 2 months, so this is a very accurate age of the person for the purposes of this reearch
}
  
dim(datExpG)
head(datExpG)

for (j in 23:36) {
  cat("\n\n\n")  # skip 3 lines
  print(paste(colnames(datExpG)[j],":",sep=""))    # print column name so we know what the variable is for the summary command that is in the next line of the script
  print(summary(aov(Scores ~ datExpG[,j], datExpG)))
}

# combine Sun and Moon and Interaction:
cat("\n\n\n")  # skip 3 lines
print(paste(colnames(datExpG)[23],"-",colnames(datExpG)[24],":",sep=""))    # print column name so we know what the variable is
print(anova(lm(Scores ~ datExpG[,23]*datExpG[,24], datExpG)))

# combine Sun, Moon and Ascendant and Interaction:
cat("\n\n\n")  # skip 3 lines
print(paste(colnames(datExpG)[23],"-",colnames(datExpG)[24],"-",colnames(datExpG)[34],":",sep=""))    # print column name so we know what the variable is
print(anova(lm(Scores ~ datExpG[,23]*datExpG[,24]*datExpG[,34], datExpG)))

# histogram of introversion-extraversion scores:
dev.new() # make new resizable window for the graph
hist(datExpG$Scores,breaks=23,xlim=c(0,25),xlab='Introversion-Extraversion Scores')
# save file:
png('IntExtScores_Histogram.png')
hist(datExpG$Scores,breaks=23,xlim=c(0,25),xlab='Introversion-Extraversion Scores')
dev.off()

# shapiro.test(datExpG$Scores)

ProducePlot <- function(FilNam){
  dev.new() # make new resizable window for the graph
  plot(model, las = 1, cex.axis = 0.6)
  png(FilNam, width=480, height=960)
  plot(model, las = 1, cex.axis = 0.8)
  dev.off()
}

ProducePlotShortHeight <- function(FilNam){
  dev.new() # make new resizable window for the graph
  plot(model, las = 1, cex.axis = 0.6)
  png(FilNam, width=480, height=480)
  plot(model, las = 1, cex.axis = 0.8)
  dev.off()
}

model = TukeyHSD( aov(Scores ~ SunSign, datExpG), "SunSign", conf.level=.95)
ProducePlot('IntExtScores_SunSign.png')

model = TukeyHSD( aov(Scores ~ MooSign, datExpG), "MooSign", conf.level=.95)
ProducePlot('IntExtScores_MooSign.png')

model = TukeyHSD( aov(Scores ~ MerSign, datExpG), "MerSign", conf.level=.95)
ProducePlot('IntExtScores_MerSign.png')

model = TukeyHSD( aov(Scores ~ VenSign, datExpG), "VenSign", conf.level=.95)
ProducePlot('IntExtScores_VenSign.png')

model = TukeyHSD( aov(Scores ~ MarSign, datExpG), "MarSign", conf.level=.95)
ProducePlot('IntExtScores_MarSign.png')

model = TukeyHSD( aov(Scores ~ JupSign, datExpG), "JupSign", conf.level=.95)
ProducePlot('IntExtScores_JupSign.png')

model = TukeyHSD( aov(Scores ~ AscSign, datExpG), "AscSign", conf.level=.95)
ProducePlot('IntExtScores_AscSign.png')

model = TukeyHSD( aov(Scores ~ SatSign, datExpG), "SatSign", conf.level=.95)
ProducePlot('IntExtScores_SatSign.png')
dev.new() # make new resizable window for the graph
counts <- table(datExpG$SatSign)
barplot(counts, main="Saturn in Signs", xlab="Zodiac Signs", col="darkblue")
# save file:
png('IntExtScores_SatSignCount.png')
barplot(counts, main="Saturn in Signs", xlab="Zodiac Signs", col="darkblue")
dev.off()

bplotcolors = c("cadetblue1","cornflowerblue","darkblue","aquamarine3",
                "chartreuse","bisque","darkgoldenrod1","chocolate",
                "brown1","brown4","aliceblue","darkgray")
Nums1to12 = c(1,2,3,4,5,6,7,8,9,10,11,12)
AgeBySign = table(datExpG$SatSign, as.integer(datExpG$AgeOfPerson))
barplot(AgeBySign,col=bplotcolors)
legend(35,160,Nums1to12,fill=bplotcolors)
# save file:
png('IntExtScores_SatSignByAgeOfPerson.png')
barplot(AgeBySign,col=bplotcolors)
legend(35,160,Nums1to12,fill=bplotcolors)
dev.off()

DrawLineOfMeansOnPlot <- function(){
  AvgAgeOfSigns = rep(0,dim(AgeBySign)[1])  # create vector (array) of 12 numbers
  for (i in 1:dim(AgeBySign)[1]) {
    NumPeople = sum(AgeBySign[i,]) #
    AvgAgeOfSigns[i] = sum(as.integer(colnames(AgeBySign)) * AgeBySign[i,]) / NumPeople
  }
  Nums1to12Vector = seq(from = as.integer(rownames(AgeBySign)[1]), to = as.integer(rownames(AgeBySign)[1])+dim(AgeBySign)[1]-1) # for example, if planet is Pluto and it is only in signs Cancer to Libra then Nums1to12Vector = 4 5 6 7 (for 4th to 7th signs of the zodiac)
  lines(Nums1to12Vector,AvgAgeOfSigns,lty=1,lwd=3,col="red")
  dev.off()
}

# Saturn:
dev.new()
png('IntExtScores_SatSignByAgeOfPerson_ScatterPlotWithMeans.png',width=960,height=640)
plot(jitter(as.integer(datExpG$SatSign),1.5), datExpG$AgeOfPerson, xlab="Saturn Signs", type="p", ylab="Age of Person", pch=20, col="black") # + geom_uitter(width=2.0,height=0.0) # scatterpolot of x=Saturn sign and y=age in years
DrawLineOfMeansOnPlot()

# Uranus:
model = TukeyHSD( aov(Scores ~ UraSign, datExpG), "UraSign", conf.level=.95)
ProducePlot('IntExtScores_UraSign.png')
dev.new()
png('IntExtScores_UraSignByAgeOfPerson_ScatterPlotWithMeans.png',width=960,height=640)
AgeBySign = table(datExpG$UraSign, as.integer(datExpG$AgeOfPerson))
plot(jitter(as.integer(datExpG$UraSign),1.5), datExpG$AgeOfPerson, xlab="Uranus Signs", type="p", ylab="Age of Person", pch=20, col="black") # + geom_jitter(width=2.0,height=0.0) # scatterpolot of x=Saturn sign and y=age in years
DrawLineOfMeansOnPlot()

# Neptune:
model = TukeyHSD( aov(Scores ~ NepSign, datExpG), "NepSign", conf.level=.95)
ProducePlot('IntExtScores_NepSign.png')
dev.new()
png('IntExtScores_NepSignByAgeOfPerson_ScatterPlotWithMeans.png',width=960,height=640)
AgeBySign = table(datExpG$NepSign, as.integer(datExpG$AgeOfPerson))
plot(jitter(as.integer(datExpG$NepSign),1.5), datExpG$AgeOfPerson, xlab="Neptune Signs", type="p", ylab="Age of Person", pch=20, col="black") # + geom_jitter(width=2.0,height=0.0) # scatterpolot of x=Saturn sign and y=age in years
DrawLineOfMeansOnPlot()

# Pluto: 
model = TukeyHSD( aov(Scores ~ PluSign, datExpG), "PluSign", conf.level=.95)
ProducePlot('IntExtScores_PluSign.png')
dev.new()
png('IntExtScores_PluSignByAgeOfPerson_ScatterPlotWithMeans.png',width=960,height=640)
AgeBySign = table(datExpG$PluSign, as.integer(datExpG$AgeOfPerson))
plot(jitter(as.integer(datExpG$PluSign),1.5), datExpG$AgeOfPerson, xlab="Pluo Signs", type="p", ylab="Age of Person", pch=20, col="black") # + geom_jitter(width=2.0,height=0.0) # scatterpolot of x=Saturn sign and y=age in years
DrawLineOfMeansOnPlot()

# Anova for Elements (Fire - Earth - Air - Water)
for (j in 37:50) {
  cat("\n\n\n")  # skip 3 lines
  print(paste(colnames(datExpG)[j],":",sep=""))    # print column name so we know what the variable is for the summary command that is in the next line of the script
  print(summary(aov(Scores ~ datExpG[,j], datExpG)))
}

model = TukeyHSD( aov(Scores ~ TNoSign, datExpG), "TNoSign", conf.level=.95)
ProducePlot('IntExtScores_TNoSign.png')

model = TukeyHSD( aov(Scores ~ AscSign, datExpG), "AscSign", conf.level=.95)
ProducePlot('IntExtScores_AscSign.png')

model = TukeyHSD( aov(Scores ~ MCSign, datExpG), "MCSign", conf.level=.95)
ProducePlot('IntExtScores_MCSign.png')

model = TukeyHSD( aov(Scores ~ VtxSign, datExpG), "VtxSign", conf.level=.95)
ProducePlot('IntExtScores_VtxSign.png')

model = TukeyHSD( aov(Scores ~ SunElem, datExpG), "SunElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_SunElem.png')

model = TukeyHSD( aov(Scores ~ MooElem, datExpG), "MooElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_MooElem.png')

model = TukeyHSD( aov(Scores ~ MerElem, datExpG), "MerElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_MerElem.png')

model = TukeyHSD( aov(Scores ~ VenElem, datExpG), "VenElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_VenElem.png')

model = TukeyHSD( aov(Scores ~ MarElem, datExpG), "MarElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_MarElem.png')

model = TukeyHSD( aov(Scores ~ JupElem, datExpG), "JupElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_JupElem.png')

model = TukeyHSD( aov(Scores ~ SatElem, datExpG), "SatElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_SatElem.png')

model = TukeyHSD( aov(Scores ~ UraElem, datExpG), "UraElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_UraElem.png')

model = TukeyHSD( aov(Scores ~ NepElem, datExpG), "NepElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_NepElem.png')

model = TukeyHSD( aov(Scores ~ PluElem, datExpG), "PluElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_PluElem.png')

model = TukeyHSD( aov(Scores ~ TNoElem, datExpG), "TNoElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_TNoElem.png')

model = TukeyHSD( aov(Scores ~ AscElem, datExpG), "AscElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_AscElem.png')

model = TukeyHSD( aov(Scores ~ MCElem, datExpG), "MCElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_MCElem.png')

model = TukeyHSD( aov(Scores ~ VtxElem, datExpG), "VtxElem", conf.level=.95)
ProducePlotShortHeight('IntExtScores_VtxElem.png')

# *** START ANCOVA *****
# BELOW IS ANCOVA ANALYSIS WHERE THE AGE OF THE PERSON IS THE COVARIATE CONTROLLED TO SEE IF THE ELEMENT (Fir/Ear/Air/Wat) affect Extraversion-Introversion with Age as the controlled covariate

print("ANCOVA ANALYSIS STARTS HERE (AGE IS THE COVARIATE\n\n")

print("A requirement of Ancova is that the covariate AgeOfPerson and the signs and elements of planets are not correlated (p-value below > .05). Note: If, for example, signs correlate with AgeOfPerson but elements do not, then Ancova can be used for elements with AgeOfPerson as the covariate but zodiac signs cannot be used in Ancova with AgeOfPerson as the covariate.")

for (j in 23:50) {
  cat("\n\n\n")  # skip 3 lines
  print(paste(colnames(datExpG)[j],":",sep=""))    # print column name so we know what the variable is
  # print(summary(aov(datExpG$AgeOfPerson ~ datExpG[,j], datExpG)))
  print(summary(aov(AgeOfPerson ~ datExpG[,j], datExpG)))
}

# Ancova: see if signs of planets predict extraveresion-introversion scores with age of the person controlled

library(car)

for (j in 23:50) {
  cat("\n\n\n")  # skip 3 lines
  print(paste(colnames(datExpG)[j],":",sep=""))    # print column name so we know what the variable is
  ancova_model <- aov(Scores ~ datExpG[,j] + AgeOfPerson, data = datExpG)
  AncovaModel = Anova(ancova_model, type="III")
  print(AncovaModel) 
}

# Age by extraversion score plot:
dev.new()
#png('IntExtScores_SatSignByAgeOfPerson_ScatterPlotWithMeans.png',width=960,height=640)
plot(datExpG$AgeOfPerson, jitter(datExpG$Scores,1.8), xlab="Age of Person", type="p", ylab="Introv-Extrav Score", pch=20, col="black") # + geom_uitter(width=2.0,height=0.0) # scatterpolot of x=Saturn sign and y=age in years

# Draw line of means
ScoreByAge = table(datExpG$Scores, as.integer(datExpG$AgeOfPerson))
AvgScoreByAge = rep(0,nrow(ScoreByAge))

for (i in 1:ncol(ScoreByAge)) {  # loop through ages of people, which ranges from 17 to 57 years old
  NumPeopleOfThisAge = sum(ScoreByAge[,i]) # number of people of this age
  iSumOfScores = 0
  for (j in 1:nrow(ScoreByAge)) {  # loop through the scores for people of this age
    iSumOfScores = iSumOfScores + as.integer(rownames(ScoreByAge)[j])*ScoreByAge[j,i] # iSumOfScores = total of scores for people of this age
  }
  AvgScoreByAge[i]=iSumOfScores/NumPeopleOfThisAge
}

lines(as.integer(colnames(ScoreByAge)),AvgScoreByAge,lty=1,lwd=3,col="red")