#####Code for Principal Components Analysis


### Import Data


datumFull=read.csv(file.choose())


###start with just Mass and length


datum=data.frame(Mass=datumFull$Mass,Body=datumFull$BodyLength)
# need to scale and center data for later use
CentDatum=data.frame(Mass=scale(datum$Mass,scale=FALSE),Body=scale(datum$BodyLength,scale=FALSE))
ScaleDatum=data.frame(Mass=scale(datum$Mass),Body=scale(datum$BodyLength))
#plot(ScaleDatum)


### Analyze Data


# First look at relationship between body length and body mass
# Need independent variables to look at effect of morphometrics on some response
plot(datum$Body,datum$Mass) ### look at relationship between body length and body mass
help(princomp) ### use princomp, which does R-factor PCA (rarely use q-factor PCA in ecology)
results=princomp(datum) ### runs the PCA
summary(results) # gives proportion of variance explained by principle components
results$loadings # helps to determine what each principle component means 
# describes how each variable contributes to each principle component
#Note negatives


###plot results


plot(CentDatum) ### plot centered data
abline(0,) #insert loadings for PC1 y-axis/x-axis ### plot first principle component
abline(0,) #insert loadings for PC2 y-axis/x-axis ### plot second principle component
#Why not perpendicular? If orthagonal, would be perpendicular
# due to differences in scale


### plot results again, but using same aspect for x and y


plot(CentDatum,asp=1) ### plot centered data with x and y on same scale
abline(0,) #insert loadings for PC1 y-axis/x-axis ### plot first principle componenet
abline(0,) #insert loadings for PC2 y-axis/x-axis ### plot second principle component


### Principal component scores for each sample


results$scores ### use these for further analyses on orthagonal variables


### Analyze again but now include ear length


datum=data.frame(Mass=datumFull$Mass, Body=datumFull$BodyLength,Ear=datumFull$Ear)
plot(datum$Body,datum$Ear)
CentDatum=data.frame(Mass=scale(datum$Mass,scale=FALSE),Body=scale(datum$BodyLength,scale=FALSE),Ear=scale(datum$Ear,scale=FALSE))
ScaleDatum=data.frame(Mass=scale(datum$Mass),Body=scale(datum$BodyLength),Ear=scale(datum$Ear))
#plot(ScaleDatum)
results=princomp(datum)
summary(results) 
results$loadings
# Note that 3rd PCA (mostly ear) explains very small amount of variance
# Not because 'EAR' isn't important, but because 'Ear' has low variance relative to other data
# Not something that needs to be 'fixed', but if reducing number of variables, unfair to 'ear'
# Best to scale data
ScaleResults=princomp(ScaleDatum) ###PCA of scaled results. 
summary(ScaleResults)
ScaleResults$loadings ### However, can make interpretation of loadings more difficult