Weekly Log

Week 1

Reading Notes

Book: Using Geochemical Data (Rollinson and Pease)

1. Normalising values can be chosen from 2 sets:
   i)CI chondrites
   ii)Ordinary chrondrites (volatile free, tend to have higher values).
   Therefore, it is important to cite the right data set. See pg.130 for dataset.
   Shale normalisation is used for upper crust sediments while Chondrite normalisation will be used in our project since we want to compare igneous rocks to primitive solar system.

2. REE in minerals in Igneous Rocks:
   i)Feldspars: low concentration, but when in felsic rocks, Eu anomaly is positive
   ii)Olivine: very low concentration, unlikely to fractionate REE while fractionation/partial melting.
   iii)Pyroxenes: orthopyroxenes have lower REE than clinopyroxenes. In mafic rocks, REE are mostly present in clinopyroxene. Pyroxenes have small negative anomalies but maybe be higher in clinopyroxenes.
   iv)Hornblende: enriched in middle REE, have higher concentration of HREE that forms a U-shaped REE pattern in melt.
   v)Garnet: high abundance of HREE to the extent that they are compatible whilst LREE aren’t.
   vi)Muscovite: enriched in REE, especially in felsic rocks. negative Eu anomaly.
   vii)Biotite: less enriched in REE, negative Eu anomaly.
   viii)Zircon: deplete a melt in HREE.
   ix)Other accessory phases: strong influence on REE pattern even if they have low abundance. They have high partition coefficients which result in disproportionate influence and depletion in REE.

3. Multi-element Diagrams (MeD):
   i)Chondrite-Normalised MeD: advantageous since it uses directly measured values instead of calculated ones. It might not be very suitable for sediment and evolved igneous rocks ( do they mean metamorphic?). They have plotted REE and and some other trace elements in order of increasing mantle compatibility on log scale using normalised values with McDonough and and Sun(1995) and Palme and ONeil(2014) as reference values.
   ii)In evolved igenous rocks, anomalies may be controlled by specific mineral as mentioned in point 2 above where LREE are controlled by allanite and HREE by garnet.

Code for Week 1: Spider Plot Function

In a Spider Diagram, we will be representing the log normalised concentrations of Rare Earth Elements. The following Elements will be plotted: La-Lanthanum, Ce- Cerium, Pr- Praseodymium, Nd- Neodymium, Pm- Promethium, Sm- Samarium, Eu- Europium, Gd- Gadolinium,, Tb- Terbium, Dy- Dysprosium, Ho- Holmium, Er- Erbium, Tm- Thulium, Yb- Ytterbium, Lu- Lutetium

Data

Samples

Spider Plot Function

#this is a function for plotting spider plots for REE data
plot_spider <- function(samples,reference_values,logscale=TRUE,colors=NULL,main= "Spider Plot for Normalised REE Values",xlab="Rare Earth Elements",ylab ="Normalised values")
{
  # these are the REE we will be covering generally
  REE_symbols <-  c("La","Ce","Pr","Nd","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu")
  if (is.null(colors)) 
  {
  colors <- c("black","red","yellow","lightblue","green","orange","pink","cyan","purple","brown","gold","violet","navy")
  }
  # normalising values of samples with the reference set
  normalised_values <- lapply(samples,function(s1) as.numeric(s1)/as.numeric(reference_values))
  #axis range
  ymin <- if (logscale) min(unlist(normalised_values), na.rm = TRUE) else 0
  ymax <- max(unlist(normalised_values), na.rm = TRUE)
  # generate plot 
  plot(1:length(reference_values),normalised_values[[1]],main= "Spider Plot",type="b",xaxt="n",,xlab="Rare Earth Elements",ylab ="Concentration (PPM)",log = 'y',ylim = c(ymin, ymax),col=colors[1])
  axis(1,at=1:length(reference_values),labels = REE_symbols)
  grid(nx=NULL, ny=NULL,col="gray",lwd=1)
  
  abline(h=1,col="black",lwd=3)
  
  # generate other plots 
   if(length(normalised_values)>1)
   { 
       for(i in 2:length(normalised_values))
       {
         lines(1:length(reference_values),normalised_values[[i]],type="b",col=colors[i])
       }
   }
   legend("bottomleft", legend = rock_names_graph, col = colors[1:length(samples)], lty = 1, lwd = 2, bty = "n")
}

Week 2

Reading Notes

Research Paper- O’Neill(2016) :“The Smoothness and Shapes of Chondrite-normalized Rare Earth Element Patterns in Basalts”

1. Normalised REE patterns are important to study since they provide information about behaviour of incompatible trace elements.
2. Spider Plots are problematic since they ‘destroy’ the shape of the REE pattern.
3. Aim of this paper is to quantify REE shapes precisely with a large number of datasets.
4. Lanthanide Contraction: 5s and 5p shells contract towards the nucleus due increased nuclear charge as a result of the 4f electrons present in the shell before them.
5. The 3+ ionic state governs chemical properties of REE.
6. Z and ionic radius can both be used as independent variables to plot normalised REE data to produce smooth graphs.
7. Y/Ho values are almost constant due similar ionic radii.
8. Anomaly in smoothness is due to Eu2+ state which has a half filled subshell providing stability.
9. If Eu2+ /Eu3+ is measurably greater than zero in a petrogenetic process, the Eu partition coefficient for that process will deviate(due to Eu anomaly) from the smooth pattern estalished by the other REE

Research Paper- Coryell,Chase & Winchester(1963) :“A Procedure for Geochemical Interpretation of Terrestrial Rare-Earth Abundance Patterns”

1. The zig–zag pattern in REE concentration plots emerges from relative enrichment or depletion of REE during natural processes. This zig zag pattern can be removed by normalising the REE abundance with bronzite chondrite REE concentrations.
2. Eu abundance is considerably low for most minerals(by factor of 14 for samarskite, 3.7 for euxenite, 2.6 for monazite, 5.5 for xenotime) due to it being divalent at some stage.
3. The technique of extrapolation of the best fit line in a normalised graph is of great value as it helps deteremine values of other REE that might have not been measured in the sample. It also makes it possible to determine REE abundance in whole new rock.
   
4. Kirovograd Granite:( How do you interpret these graphs? )
   i)HREE very high in garnet and apatite.
   ii)Biotite, Monazite, Feldspar+Quartz decrease with increase in Z.
5. LREE are more enriched than HREE in basalt.

Research Paper- Anenburg(2020) :“Rare earth mineral diversity controlled by REE pattern shapes”

1. Aim: use lambda combinations to construct CI normalised REE plots in such a way that REE patterns for natural materials become predictable.
   
2. La3+: large ionic radii, Lu3+: quite small due to lanthanide contraction.
   
3. LREE: La to Sm.
   
4. MREE: Sm to Dy.
   
5. HREE: Eu to Lu.
   
6. Oddo–Harkins Rule: elements with even atomic numbers are more abundant than their neighbouring odd-numbered elements. 7.λ0 is the overall abundance of the REE.
   
7. λ1 is the linear slope (most instructive).
   
8. λ2 is the quadratic curvature (most instructive).
   
9. λ3 represents the inflections at the ends of patterns (sinusoidality).
   
10. λ4 term is required if any ‘W’ pattern is present.
    
11. LaN/LuN is a measure of the overall slope of the pattern.LaN/ SmN is a measure of the LREE enrichment of a pattern. These are not always useful to represent sinusoidality or curves which makes authors use differnet ratios leading to too much confusion and biases.
    
12. They have generated REE patterns for λ values..
    
13. Ce and Y dominate the REE pattern, just like the real world where they make up nearly 75.7% of REE minerals.Without Ce and Y other REE like La,Dy,Yb,Nd dominate. Yb is most dominant in these Ce and Y absent graphs.
    
14. Eu is present in high concentrations in Ca-rich minerals like plagioclase and epidote. Eu anomaly needs to be at a certain level to be dominant for different λ values. Eu rich mineral could be found in anorthosites or epidosites.
15. If Y is removed from the Levison suffixes, more REE minerals can be discovered.

Week 3

Notes

Using AlambdaR website to understand effects to change in λ values on the shape of the Normalised Values vs REE radii

We have taken a negative Eu anomaly and positive Ce anomaly for all cases

1. Changing λ0: When other coefficients are kept constant, changing λ0 changes the CI-normalised concentration by a factor and no change in the shape of the curve is observed. The entire curve goes up if λ0 is increased and goes down if λ0 is decreased.
   
2. Changing λ1: When other coefficients are kept constant, changing λ1 to a positive number or increasing it causes the concentration of HREEs to decrease while the LREEs are more abundant. When λ1 is decreased or changed to a negative number, the trend flips and HREEs become more abundant while LREEs have a lower concentration.
   
3. Changing λ2: When other coefficients are kept constant, increasing λ2 makes the Eu anomaly (at the vertex) is more prevelant as the curve becomes U-shaped withe the MREEs having low concentration while the LREEs and HREEs are more abundant. When λ2 is 0 or a small number (between -200 & +200) the curve is relatively straight and not parabolic. When λ2 is a large negative number, the MREEs have a higher abundance than LREEs and HREEs but negative Eu anomaly persists.
   
4. Changing λ3: adds to the cruvature, concentration doesn’t change much even if it. is changed from +1000 to -1000. However, the sinusodiality becomes prevelant if λ3 is a very big or small number. It is hard to predict as clearly as prevoius cases how λ3 affects concentration.
   
5. Changing λ4: when other λ are kept in a reasonable value range, λ4 produces a W-shape only when it is the order of 10^4 which seems like a very unreasonable value.

Data Simplifying Function

data_simplifying <- function(input_data)
{
  REE_list <- c("LATITUDE MIN","LONGITUDE MIN","LATITUDE MAX","LONGITUDE MAX","ROCK NAME","SIO2(WT%)","LA(PPM)","CE(PPM)","PR(PPM)","ND(PPM)","SM(PPM)","EU(PPM)","GD(PPM)","TB(PPM)","DY(PPM)","HO(PPM)","ER(PPM)","TM(PPM)","YB(PPM)","LU(PPM)")
  complete_row <- complete.cases(input_data[, REE_list])
  refined_REE_data <- input_data[complete_row,REE_list]
  
} 

Week 4

Simplification of all data

rock_files_csv <- list.files(path="~/Desktop/MSci Project/Rock Files CSV",full.names = TRUE)
rock_list_read <- lapply(rock_files_csv,read.csv,check.names=FALSE)
rock_files_names_extract <- basename(rock_files_csv)
rock_files_names <- sub(".csv","",rock_files_names_extract)
names(rock_list_read) <- rock_files_names
zero_containing_simplified_rocks <- lapply(rock_list_read,data_simplifying)
non_zero <- function(input)
{
  keep_row <- rowSums(input==0,na.rm = TRUE)==0
  return(input[keep_row, ])
}
simplified_rocks <- lapply(zero_containing_simplified_rocks,non_zero)

Saving simplified data in new folder

Removing bias for further study by taking random 100 samples

Week 5

Further Simplification of Data

output_directory <- "Simplified_Data"
dir.create(path = output_directory)

## Warning in dir.create(path = output_directory): 'Simplified_Data' already
## exists

set.seed(21)
for(rock_name in names(simplified_rocks))
{
  data_for_loop <- simplified_rocks[[rock_name]]
  nrows <-nrow(data_for_loop)
  selected_data <- NULL
  sample_size <- 100
  if (nrows >= sample_size)
  {
    random_rows <- sample(1:nrows,sample_size,replace = FALSE)
    selected_random_rows <- data_for_loop[random_rows, ]
    selected_data <- selected_random_rows
  }
  else
  {
    selected_data<- data_for_loop
  }
  # Averaging Lat and Lon 
   avg_lat <- (selected_data$"LATITUDE MIN"+selected_data$"LATITUDE MAX")/2
   lon_min_rad <- (selected_data$"LONGITUDE MIN" * pi/180)
   lon_max_rad <- (selected_data$"LONGITUDE MAX" * pi/180)
   x1 <- cos(lon_min_rad)
   x2 <-cos(lon_max_rad)
   y1 <-sin(lon_min_rad)
   y2 <-sin(lon_max_rad)
   x <- (x1 +x2)/2
   y <- (y1 +y2)/2
   avg_rad <- atan2(y,x)
   
   avg_long <- avg_rad*180/pi
   
   removed_col <- c("LATITUDE MIN","LONGITUDE MIN","LATITUDE MAX","LONGITUDE MAX")
   kept_col <- selected_data[, !colnames(selected_data) %in% removed_col]
   kept_col$AVERAGE_LAT <- avg_lat
   kept_col$AVERAGE_LON <- avg_long
   file_name <- paste0(rock_name,".csv")
   saving_path <- file.path(output_directory,file_name)
   write.csv(kept_col,file = saving_path,row.names = FALSE)
   
}

Reading Week

Saving Combined Rock File

new_output_directory <- "Combined_Data_File"
smaller_files <- list.files(path=output_directory,pattern = ".csv$", full.names = TRUE)
combined_list <- lapply(smaller_files,read.csv,encoding = "UTF-8")
combined_data <- do.call(rbind,combined_list)
combined_data$ROCK.NAME <- gsub(" *\\[.*?\\]", "", combined_data$ROCK.NAME, useBytes = TRUE)
saving_directory <- file.path(new_output_directory,"Combined_rock_data.csv")
write.csv(combined_data,file = saving_directory,row.names = FALSE)

Log-ratio Transformation

non_ree_col <- c("ROCK.NAME","AVERAGE_LAT","AVERAGE_LON","SIO2.WT..")
ree_data <- combined_data[, !colnames(combined_data) %in% non_ree_col]
non_ree_data <- combined_data[,non_ree_col]

#log transformation
log_ree_data_no_other <- log(ree_data)
log_mean_no_other <-  apply(log_ree_data_no_other,1,mean)
#ree_transformed <- matrix(data= NA,nrow=nrow(log_ree_data),ncol=ncol(log_ree_data))

# for(i in 1:nrow(log_ree_data))
# {
#   row_ree_data <- log_ree_data[i,]
#   row_ree_data_mean <- log_mean[i]
#   ree_transformed[i,] <- row_ree_data - row_ree_data_mean
# }

ree_transformed_no_other <- sweep(log_ree_data_no_other,1,log_mean_no_other,FUN = "-")

Week 6

PCA

pca_no_other <- prcomp(ree_transformed_no_other)
current_names_no_other <- row.names(pca_no_other$rotation)
clean_names_no_other <- sub("\\.[A-Z]+\\.?$", "", current_names_no_other)
row.names(pca_no_other$rotation) <- clean_names_no_other
biplot(pca_no_other,xlabs=rep('',nrow(ree_transformed_no_other)),main="REE Biplot")

plot(pca_no_other,main = "Scree Plot for REE Biplot")

Initial Interpretations of the Biplot

1. There is good correlation amongst the LREEs.
   
2. There is good correlation amongst the HREEs.
   
3. Eu concentration is not related to LREEs and other HREEs but is more related to MREEs. It has the most isolated concentration.
   
4. Correlation among MREE end members is weak.
   
5. The significance of the Oddo Harkins Rule is not very prevelant in the biplot. In traditional Spider plots, it is easy to see that odd atomic numbered elements have lower concentration than even numbered which is not the case here.
   
6. Negatively correlated REE pairs: La and Dy, Ce and Dy, Pr and Dy, Nd and Ho, Sm and Lu, Yb and Sm
7. Eu and Gd are better represented by PC2 while most other REEs are represented by PC1.
8. Gd, Sm and Tb are not well represented as compared to others due to their short arrow length

Initial Interpretations of the Screeplot

1. From the screeplot it is evident that PC1 is the most important Principal Component variable. Together with PC2, it covers xx% of the variance.
2. Beyond PC3, the variables do not explain the variance much and can be ignored.

Contribution of each PC in Screeplot 1

#round((pca_no_other$sdev^2/sum(pca_no_other$sdev^2))*100)

PC 1: 79% PC 2: 10% PC 3: 4% PC 4: 2% PC 5- PC 8: 1% Beyond PC 8: negligible

Week 7

PC vs Atomic Number Plots

PC_data_no_other <- pca_no_other$rotation
for (i in 1:5)
{
  pc_number_no_other <- paste0("PC",i)
  PC_no_other <- PC_data_no_other[,pc_number_no_other]
  plot(PC_no_other,ylab = pc_number_no_other,xlab = "Rare Earth Elements",xaxt="n")
  axis(side=1,at=1:length(PC_no_other),labels = clean_names_no_other,cex.axis=0.8)
}

Biplot without Eu

remove_Eu <- "EU.PPM."
ree_without_Eu_no_other <- ree_transformed_no_other[, !names(ree_transformed_no_other) %in% remove_Eu]
pca_no_Eu_no_other <- prcomp(ree_without_Eu_no_other)
biplot(pca_no_Eu_no_other,xlabs=rep('',nrow(ree_without_Eu_no_other)),main="Biplot 2")

plot(pca_no_Eu_no_other)

Contribution of each PC in Screeplot without Eu

#round((pca_no_Eu_no_other$sdev^2/sum(pca_no_Eu_no_other$sdev^2))*100)

PC 1: 85% PC 2: 7% PC 3: 3% PC 4: 2% PC 5- PC 8: 1% Beyond PC 8: negligible

Week 8

Replotting biplots with ‘Other’ column

#making 'other'column
total_ree_per_sample <- rowSums(ree_data)
other_col <- pmax(0,1000000 - total_ree_per_sample)
ree_data$Other <- other_col
combined_data$Other <- other_col
#log transformation
log_ree_data <- log(ree_data)
log_mean <-  apply(log_ree_data,1,mean)
ree_transformed <- sweep(log_ree_data,1,log_mean,FUN = "-")

pca <- prcomp(ree_transformed)
current_names <- row.names(pca$rotation)
clean_names <- sub("\\.[A-Z]+\\.?$", "", current_names)
row.names(pca$rotation) <- clean_names
biplot(pca,xlabs=rep('',nrow(ree_transformed)),main= "Biplot 3")

plot(pca)

Contribution of each PC in Screeplot 3

#round((pca$sdev^2/sum(pca$sdev^2))*100)

PC 1: 66% PC 2: 20% PC 3: 7% PC 4: 2% PC 5- PC 7: 1% Beyond PC 7: negligible ## PC Plots with Other

PC_data <- pca$rotation
for (i in 1:5)
{
  pc_number <- paste0("PC",i)
  PC <- PC_data[,pc_number]
  plot(PC,ylab = pc_number,xlab = "Rare Earth Elements",xaxt="n")
  axis(side=1,at=1:length(PC),labels = clean_names,cex.axis=0.8)
}

Biplot without Eu

remove_Eu <- "EU.PPM."
ree_withou_Eu <- ree_transformed[, !names(ree_transformed) %in% remove_Eu]
pca_no_Eu <- prcomp(ree_withou_Eu)
biplot(pca_no_Eu,xlabs=rep('',nrow(ree_withou_Eu)),main="Biplot 4")

plot(pca_no_Eu)

## Contribution of each PC in Screeplot 4 without Eu

#round((pca_no_Eu$sdev^2/sum(pca_no_Eu$sdev^2))*100)

PC 1: 69% PC 2: 21% PC 3: 5% PC 4: 2% PC 5- PC 7: 1% Beyond PC 7: negligible

PC Plots with Atomic Number as label

atomic_number <- c(57,58,59,60,62,63,64,65,66,67,68,69,70,71,"Other")
xcoordinate <- numeric(length(clean_names))
other_coordinate <- which(clean_names=="Other")

xcoordinate[clean_names != "Other"] <- atomic_number

## Warning in xcoordinate[clean_names != "Other"] <- atomic_number: number of
## items to replace is not a multiple of replacement length

xcoordinate[other_coordinate] <- 73

axis_ticks <- 57:71
axis_labels <- as.character(axis_ticks)
axis_labels[axis_ticks==61] <-"61"
ticks <- c(axis_ticks, 73) 
label <- c(axis_labels,"Other")

for (i in 1:5)
{
  pc_number <- paste0("PC",i)
  PC <- PC_data[,pc_number]
  plot(xcoordinate,PC, ylab = pc_number,xlab = "Atomic Number",xaxt="n")
  axis(side=1,at= ticks,labels= label,cex.axis=0.8)
}

PC Plots with Ionic radius as label

Ionic Radius values from Shannon-1976

ionic_radius8 <- c(1.160,1.143,1.126,1.109,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
xcoordinate <- numeric(length(clean_names))
other_coordinate <- which(clean_names=="Other","")

xcoordinate[clean_names != "Other"] <- ionic_radius8
xcoordinate[other_coordinate] <- 0.95

axis_labels_ionic8 <- as.character(round(ionic_radius8, 3))
ticks <- xcoordinate
label <- c(axis_labels_ionic8,"Other")

for (i in 1:5)
{
  pc_number <- paste0("PC",i)
  PC <- PC_data[,pc_number]
  par(las=2)
  plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=8) ",xaxt="n")
  axis(side=1,at= ticks,labels= label,cex.axis=0.8)
}

ionic_radius6 <- c(1.032,1.010,0.990,0.983,0.958,0.947,0.938,0.923,0.912,0.901,0.890,0.880,0.868,0.861)
xcoordinate <- numeric(length(clean_names))
other_coordinate <- which(clean_names=="Other")

xcoordinate[clean_names != "Other"] <- ionic_radius6
xcoordinate[other_coordinate] <- 0.85

axis_labels_ionic6 <- as.character(round(ionic_radius6, 3))
ticks <- xcoordinate
label <- c(axis_labels_ionic6,"Other")

for (i in 1:5)
{
  pc_number <- paste0("PC",i)
  PC <- PC_data[,pc_number]
  par(las=2)
  plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=6) ",xaxt="n")
  axis(side=1,at= ticks,labels= label,cex.axis=0.8)
}

PC2 vs PC3 Biplots

biplot(pca,choices=c(2,3),xlabs=rep('',nrow(ree_transformed)))

biplot(pca_no_Eu,choices=c(2,3),xlabs=rep('',nrow(ree_withou_Eu)))

Week 9

Switch Case

plots_switch_case <- function(input) {

  case_key <- as.character(input)

  result <- switch(case_key,

    "1" =
    {
      # Other + Eu
      pca_other_eu <- prcomp(ree_transformed)
      current_names <- row.names(pca$rotation)
      clean_names <- sub("\\.[A-Z]+\\.?$", "", current_names)
      row.names(pca$rotation) <- clean_names
      biplot(pca_other_eu,xlabs=rep('',nrow(ree_transformed)),main="Biplot PC1 vs PC2")
      biplot(pca_other_eu,choices = c(2,3),xlabs=rep('',nrow(ree_transformed)),main="Biplot PC2 vs PC3")
      plot(pca_other_eu)
      
      
      atomic_number <- c(57,58,59,60,62,63,64,65,66,67,68,69,70,71,"Other")
      xcoordinate <- numeric(length(clean_names))
      other_coordinate <- which(clean_names=="Other")
      xcoordinate[clean_names != "Other"] <- atomic_number
      xcoordinate[other_coordinate] <- 73
      axis_ticks <- 57:71
      axis_labels <- as.character(axis_ticks)
      axis_labels[axis_ticks==61] <-"61"
      ticks <- c(axis_ticks, 73) 
      label <- c(axis_labels,"Other")
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[,pc_number]
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Atomic Number",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }
      
      
      ionic_radius8 <- c(1.160,1.143,1.126,1.109,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
      xcoordinate <- numeric(length(clean_names))
      other_coordinate <- which(clean_names=="Other")
      xcoordinate[clean_names != "Other"] <- ionic_radius8
      xcoordinate[other_coordinate] <- 0.95
      axis_labels_ionic8 <- as.character(round(ionic_radius8, 3))
      ticks <- xcoordinate
      label <- c(axis_labels_ionic8,"Other")
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[,pc_number]
        par(las=2)
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=8) ",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }
      
      
      ionic_radius6 <- c(1.032,1.010,0.990,0.983,0.958,0.947,0.938,0.923,0.912,0.901,0.890,0.880,0.868,0.861)
      xcoordinate <- numeric(length(clean_names))
      other_coordinate <- which(clean_names=="Other")
      xcoordinate[clean_names != "Other"] <- ionic_radius6
      xcoordinate[other_coordinate] <- 0.85
      axis_labels_ionic6 <- as.character(round(ionic_radius6, 3))
      ticks <- xcoordinate
      label <- c(axis_labels_ionic6,"Other")
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[,pc_number]
        par(las=2)
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=6) ",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }

    },

    "2" = 
    {
    # No Other + Eu
      pca_no_other <- prcomp(ree_transformed_no_other)
      current_names_no_other <- row.names(pca_no_other$rotation)
      clean_names_no_other <- sub("\\.[A-Z]+\\.?$", "", current_names_no_other)
      row.names(pca_no_other$rotation) <- clean_names_no_other
      biplot(pca_no_other,xlabs=rep('',nrow(ree_transformed_no_other)),main="Biplot PC1 vs PC2")
      biplot(pca_no_other,choices = c(2,3),xlabs=rep('',nrow(ree_transformed_no_other)),main="Biplot PC2 vs PC3")
      plot(pca_no_other)
      
      PC_data_no_other <- data.frame(pca_no_other$rotation)
      
      atomic_number <- c(57,58,59,60,62,63,64,65,66,67,68,69,70,71)
      xcoordinate <- atomic_number
      axis_ticks <- 57:71
      axis_labels <- as.character(axis_ticks)
      axis_labels[axis_ticks==61] <-"61"
      ticks <- axis_ticks
      label <- axis_labels
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data_no_other[,pc_number]
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Atomic Number",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }
      
      
      ionic_radius8 <- c(1.160,1.143,1.126,1.109,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
      xcoordinate <- ionic_radius8
      axis_labels_ionic8 <- as.character(round(ionic_radius8, 3))
      ticks <- xcoordinate
      label <- axis_labels_ionic8
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data_no_other[,pc_number]
        par(las=2)
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=8) ",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }
      
      
      ionic_radius6 <- c(1.032,1.010,0.990,0.983,0.958,0.947,0.938,0.923,0.912,0.901,0.890,0.880,0.868,0.861)
      xcoordinate <- ionic_radius6
      axis_labels_ionic6 <- as.character(round(ionic_radius6, 3))
      ticks <- xcoordinate
      label <- axis_labels_ionic6
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data_no_other[,pc_number]
        par(las=2)
        plot(xcoordinate,PC, ylab = pc_number,xlab = "Ionic Radius (CN=6) ",xaxt="n")
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
}
    },

    "3" =
    {
      # Other + No Eu
      pca_other_no_eu <- prcomp(ree_withou_Eu)
      current_names_other_no_eu <- row.names(pca_other_no_eu$rotation)
      clean_names_other_no_eu <- sub("\\.[A-Z]+\\.?$", "", current_names_other_no_eu)
      row.names(pca_other_no_eu$rotation) <- clean_names_other_no_eu
      biplot(pca_other_no_eu,xlabs=rep('',nrow(ree_withou_Eu)),main="Biplot PC1 vs PC2")
      biplot(pca_other_no_eu,choices = c(2,3),xlabs=rep('',nrow(ree_withou_Eu)),main="Biplot PC2 vs PC3")
      plot(pca_other_no_eu)
      PC_data_other_no_eu <- data.frame(pca_other_no_eu$rotation)
      
      
      atomic_number <- c(57,58,59,60,62,64,65,66,67,68,69,70,71)
      xcoordinate <- numeric(length(clean_names_other_no_eu))
      other_coordinate <- which(clean_names_other_no_eu=="Other")
      xcoordinate[clean_names_other_no_eu != "Other"] <- atomic_number
      xcoordinate[other_coordinate] <- 73
      axis_ticks <- 57:71 
      axis_labels <- as.character(axis_ticks)
      axis_labels[axis_ticks==61] <-"61"
      ticks <- c(axis_ticks, 73) 
      label <- c(axis_labels,"Other")
       
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_other_no_eu[,pc_number]
          plot(xcoordinate,PC, ylab = pc_number,xlab = "Atomic Number",xaxt="n")
          axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }
      
      ionic_radius8 <- c(1.160,1.143,1.126,1.109,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
      ionic_radius8_plot <- ionic_radius8[c(1:5, 7:14)]
      other <- 1.20
      xcoordinate_ionic_radius <- numeric(length(clean_names_other_no_eu))
      other_coordinate_index <- which(clean_names_other_no_eu == "Other")
      xcoordinate_ionic_radius[clean_names_other_no_eu != "Other"] <- ionic_radius8_plot
      xcoordinate_ionic_radius[other_coordinate_index] <- other
      ticks<- c(ionic_radius8, other)
      label <- c(sprintf("%.3f",ionic_radius8), "Other")
      
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_other_no_eu[,pc_number]
          par(las=2)
          plot(xcoordinate_ionic_radius, PC, ylab = pc_number, xlab = "Ionic Radius (CN=8) ", xaxt="n")
          axis(side=1, at= ticks, labels= label, cex.axis=0.8)
      }
      
      
      ionic_radius6 <- c(1.032,1.010,0.990,0.983,0.958,0.947,0.938,0.923,0.912,0.901,0.890,0.880,0.868,0.861)
      ionic_radius6_plot <- ionic_radius6[c(1:5, 7:14)]
      other <- 1.05
      xcoordinate_ionic_radius <- numeric(length(clean_names_other_no_eu))
      other_coordinate_index <- which(clean_names_other_no_eu == "Other")
      xcoordinate_ionic_radius[clean_names_other_no_eu != "Other"] <- ionic_radius6_plot
      xcoordinate_ionic_radius[other_coordinate_index] <- other
      ticks<- c(ionic_radius6, other)
      label <- c(sprintf("%.3f",ionic_radius6), "Other")
      
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_other_no_eu[,pc_number]
          par(las=2)
          plot(xcoordinate_ionic_radius, PC, ylab = pc_number, xlab = "Ionic Radius (CN=6) ", xaxt="n")
          axis(side=1, at= ticks, labels= label, cex.axis=0.8)
      }

    },

    "4" = 
    {
      pca_no_eu_no_other <- prcomp(ree_without_Eu_no_other)
      current_names_no_eu_no_other <- row.names(pca_no_eu_no_other$rotation)
      clean_names_no_eu_no_other <- sub("\\.[A-Z]+\\.?$", "", current_names_no_eu_no_other)
      row.names(pca_no_eu_no_other$rotation) <- clean_names_no_eu_no_other
      num_ree <- nrow(ree_without_Eu_no_other)
      biplot(pca_no_eu_no_other, xlabs=rep('', num_ree), main="Biplot PC1 vs PC2 ")
      biplot(pca_no_eu_no_other, choices = c(2,3), xlabs=rep('', num_ree), main="Biplot PC2 vs PC3")
      plot(pca_no_eu_no_other)
      PC_data_no_eu_no_other <- data.frame(pca_no_eu_no_other$rotation)
      
      atomic_number <- c(57,58,59,60,62,64,65,66,67,68,69,70,71)
      xcoordinate <- atomic_number
      axis_ticks <- 57:71 
      axis_labels <- as.character(axis_ticks)
      axis_labels[axis_ticks == 61] <- "61"
      axis_labels[axis_ticks == 63] <- "63"
      ticks <- axis_ticks
      label <- axis_labels
      
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_no_eu_no_other[,pc_number] 
          plot(xcoordinate, PC, ylab = pc_number, xlab = "Atomic Number", xaxt="n", xlim=c(56.5, 71.5))
          axis(side=1, at= ticks, labels= label, cex.axis=0.8)
      }
      
      ionic_radius8 <- c(1.160, 1.143, 1.126, 1.109, 1.079, 1.066, 1.053, 1.040, 1.027, 1.015, 1.004, 0.994, 0.985, 0.977)
      xcoordinate_ionic_radius8 <- ionic_radius8[c(1:5, 7:14)] 
      ticks <- ionic_radius8
      label <- c(sprintf("%.3f", ionic_radius8))
      
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_no_eu_no_other[,pc_number]
          par(las=2)
          plot(xcoordinate_ionic_radius8, PC, ylab = pc_number, xlab = "Ionic Radius (CN=8) ", xaxt="n", xlim=c(0.97, 1.17))
          axis(side=1, at= ticks, labels= label, cex.axis=0.8)
      }
      
      ionic_radius6 <- c(1.032, 1.010, 0.990, 0.983, 0.958, 0.947, 0.938, 0.923, 0.912, 0.901, 0.890, 0.880, 0.868, 0.861)
      xcoordinate_ionic_radius6 <- ionic_radius6[c(1:5, 7:14)] 
      ticks <- ionic_radius6
      label <- c(sprintf("%.3f", ionic_radius6))
      
      for (i in 1:5)
      {
          pc_number <- paste0("PC",i)
          PC <- PC_data_no_eu_no_other[,pc_number]
          par(las=2)
          plot(xcoordinate_ionic_radius6, PC, ylab = pc_number, xlab = "Ionic Radius (CN=6) ", xaxt="n", xlim=c(0.85, 1.05))
          axis(side=1, at= ticks, labels= label, cex.axis=0.8)
      }
      
      
    },

    {
      stop(paste("Invalid case number specified:", input_number, ". Please use 1, 2, 3, or 4."))
    }
  )

  return(result)
}

input_number <- readline(prompt = "Enter a number: 1 for Other+Eu, 2 for No Other+Eu, 3 for Other+No Eu and 4 for No Other+ No Eu ")

case_number <- plots_switch_case(input_number)

print(case_number)

Christmas Break

# Other + Eu
      pca_other_eu <- prcomp(ree_transformed)
      current_names <- row.names(pca_other_eu$rotation)
      clean_names <- sub("\\.[A-Z]+\\.?$", "", current_names)
      row.names(pca_other_eu$rotation) <- clean_names
      biplot(pca_other_eu,xlabs=rep('',nrow(ree_transformed)),main="Biplot PC1 vs PC2")

      biplot(pca_other_eu,choices = c(2,3),xlabs=rep('',nrow(ree_transformed)),main="Biplot PC2 vs PC3")

      plot(pca_other_eu,main = "Scree Plot")

      atomic_number <- c(57,58,59,60,62,63,64,65,66,67,68,69,70,71)
      xcoordinate <- atomic_number
      axis_ticks <- 57:71
      axis_labels <- as.character(axis_ticks)
      axis_labels[axis_ticks==61] <-"61"
      ticks <- axis_ticks 
      label <- axis_labels
      element_symbol <- clean_names[clean_names != "Other"]
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[clean_names != "Other",pc_number]
        point_colour_coding <- ifelse(xcoordinate %% 2 == 0, "blue","black")
        plot(xcoordinate,PC,col=point_colour_coding, ylab = pc_number,xlab = "Atomic Number",xaxt="n")
        text(xcoordinate,PC,labels = element_symbol, pos=4, col="red", cex=0.5)
        legend("topleft",legend= c("Even","Odd"),col=c("blue","black"),pch=16,title="Atomic Number",bty="o",cex=0.5)
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }

      ionic_radius8 <- c(1.160,1.143,1.126,1.109,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
      xcoordinate <- ionic_radius8
      axis_labels_ionic8 <- as.character(round(ionic_radius8, 3))
      ticks <- xcoordinate
      label <- axis_labels_ionic8
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[clean_names != "Other",pc_number]
        par(las=2)
        point_colour_coding <- ifelse(atomic_number %% 2 == 0, "blue","black")
        plot(xcoordinate,PC,col=point_colour_coding, ylab = pc_number,xlab = "Ionic Radius (CN=8) ",xaxt="n")
        text(xcoordinate,PC,labels = element_symbol, pos=4, col="red", cex=0.5)
        legend("topright",legend= c("Even","Odd"),col=c("blue","black"),pch=16,title="Atomic Number",bty="o",cex=0.5)
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }

      ionic_radius6 <- c(1.032,1.010,0.990,0.983,0.958,0.947,0.938,0.923,0.912,0.901,0.890,0.880,0.868,0.861)
      xcoordinate <- ionic_radius6
      axis_labels_ionic6 <- as.character(round(ionic_radius6, 3))
      ticks <- xcoordinate
      label <- axis_labels_ionic6
      for (i in 1:5)
      {
        pc_number <- paste0("PC",i)
        PC <- PC_data[clean_names != "Other",pc_number]
        par(las=2)
        point_colour_coding <- ifelse(atomic_number %% 2 == 0, "blue","black")
        plot(xcoordinate,PC, col=point_colour_coding,ylab = pc_number, xlab = "Ionic Radius (CN=6) ",xaxt="n")
        text(xcoordinate,PC,labels = element_symbol, pos=4, col="red", cex=0.5)
        legend("topright",legend= c("Even","Odd"),col=c("blue","black"),pch=16,title="Atomic Number",bty="o",cex=0.5)
        axis(side=1,at= ticks,labels= label,cex.axis=0.8)
      }

Week 11

rock_row_ranges <- list(
  "Adakite" = 1:100, "Alkali basalt" = 101:200, "Andesite" = 201:300, 
  "Ankaramite" = 301:389, "Basalt" = 389:488, "Basaltic-andesite" = 489:588, 
  "Basaltic-trachyandesite" = 589:688, "Basanite" = 689:788, "Benmorite" = 789:888, 
  "Boninite" = 889:988, "Calc-alkali basalt" = 989:1088, "Carbonatite" = 1089:1188, 
  "Comenditic Rhyolite" = 1189:1288, "Dacite" = 1289:1388, "Diorite" = 1389:1488, 
  "Dolerite" = 1489:1588, "Dunite" = 1589:1688, "Foidite" = 1689:1788, 
  "Gabbro" = 1789:1888, "Granite" = 1889:1988, "Granodiorite" = 1989:2088, 
  "Harzburgite" = 2089:2188, "Hawaiite" = 2189:2288, "Kimberlite" = 2289:2388, 
  "Komatiite" = 2389:2488, "Lamproite" = 2489:2588, "Lamprophyre" = 2589:2688, 
  "Latite" = 2689:2788, "Leucitite" = 2789:2888, "Lherzolite" = 2889:2988, 
  "Melilitite" = 2989:3088, "Mugearite" = 3089:3188, "Nephelinite" = 3189:3288, 
  "Oceanite" = 3289:3314, "Ore" = 3315:3355, "Orthopyroxenite" = 3356:3455, 
  "Pantelleritic Rhyolite" = 3456:3555, "Peridotite" = 3556:3655, "Phonolite" = 3656:3755, 
  "Phonolitic-Tephrite" = 3756:3855, "Phoscorite" = 3856:3917, "Picrite" = 3918:4017, 
  "Pyroxenite" = 4018:4117, "Rhyodacite" = 4118:4217, "Rhyolite" = 4218:4317, 
  "Shoshonite" = 4318:4417, "Subalkali Basalt" = 4418:4517, "Tephrite" = 4518:4617, 
  "Tephritic-Phonolite" = 4618:4717, "Tholeiitic Andesite" = 4718:4742, 
  "Tholeiitic Basalt" = 4743:4842, "Tonalite" = 4843:4942, "Trachyandesite" = 4943:5042, 
  "Trachybasalt" = 5043:5142, "Trachydacite" = 5143:5242, "Trachyte" = 5243:5342, 
  "Transitional Basalt" = 5343:5442, "Wehrlite" = 5443:5542)

for (i in 1:15) # DO WE NEED A X-Label on plots as the rocks are already the labels ? 
{
  
  rock_PC_averages <- sapply(rock_row_ranges, function(row) mean(pca_other_eu$x[row, i], na.rm = TRUE))
  par(mar = c(10, 4, 4, 2)) 
  plot(rock_PC_averages, xaxt = "n", main = paste("Average PC", i),
       ylab = paste("Mean PC", i, "Value"),col = "red")
  axis(1, at = 1:58, labels = names(rock_PC_averages), las = 2, cex.axis = 0.7)
  grid()
}

for (i in 1:15) 
{
  rock_PC_sd <- sapply(rock_row_ranges, function(row) sd(pca_other_eu$x[row, i], na.rm = TRUE))
  par(mar = c(10, 4, 4, 2)) 
  plot(rock_PC_sd, xaxt = "n", main = paste("Standard Deviation of PC", i), 
       ylab = paste("SD of PC", i),col = "blue")
  axis(1, at = 1:58, labels = names(rock_PC_sd), las = 2, cex.axis = 0.7)
  grid()
}

# set.seed(1)
# data_for_pca <- combined_data[, sapply(combined_data, is.numeric)]
# pca_other_eu <- prcomp(data_for_pca,scale. = TRUE)
# 
# num_rocks <- 58
# num_PCs <- 15
# num_iterations <- 500
# 
# mean_val_500 <- array(NA, dim = c(num_rocks, num_PCs, num_iterations))
# sd_val_500 <- array(NA, dim = c(num_rocks, num_PCs,num_iterations))
# PC_loadings_list <- list() 
# 
# rock_row_ranges <- list(
#   "Adakite" = 1:100, "Alkali basalt" = 101:200, "Andesite" = 201:300, 
#   "Ankaramite" = 301:389, "Basalt" = 389:488, "Basaltic-andesite" = 489:588, 
#   "Basaltic-trachyandesite" = 589:688, "Basanite" = 689:788, "Benmorite" = 789:888, 
#   "Boninite" = 889:988, "Calc-alkali basalt" = 989:1088, "Carbonatite" = 1089:1188, 
#   "Comenditic Rhyolite" = 1189:1288, "Dacite" = 1289:1388, "Diorite" = 1389:1488, 
#   "Dolerite" = 1489:1588, "Dunite" = 1589:1688, "Foidite" = 1689:1788, 
#   "Gabbro" = 1789:1888, "Granite" = 1889:1988, "Granodiorite" = 1989:2088, 
#   "Harzburgite" = 2089:2188, "Hawaiite" = 2189:2288, "Kimberlite" = 2289:2388, 
#   "Komatiite" = 2389:2488, "Lamproite" = 2489:2588, "Lamprophyre" = 2589:2688, 
#   "Latite" = 2689:2788, "Leucitite" = 2789:2888, "Lherzolite" = 2889:2988, 
#   "Melilitite" = 2989:3088, "Mugearite" = 3089:3188, "Nephelinite" = 3189:3288, 
#   "Oceanite" = 3289:3314, "Ore" = 3315:3355, "Orthopyroxenite" = 3356:3455, 
#   "Pantelleritic Rhyolite" = 3456:3555, "Peridotite" = 3556:3655, "Phonolite" = 3656:3755, 
#   "Phonolitic-Tephrite" = 3756:3855, "Phoscorite" = 3856:3917, "Picrite" = 3918:4017, 
#   "Pyroxenite" = 4018:4117, "Rhyodacite" = 4118:4217, "Rhyolite" = 4218:4317, 
#   "Shoshonite" = 4318:4417, "Subalkali Basalt" = 4418:4517, "Tephrite" = 4518:4617, 
#   "Tephritic-Phonolite" = 4618:4717, "Tholeiitic Andesite" = 4718:4742, 
#   "Tholeiitic Basalt" = 4743:4842, "Tonalite" = 4843:4942, "Trachyandesite" = 4943:5042, 
#   "Trachybasalt" = 5043:5142, "Trachydacite" = 5143:5242, "Trachyte" = 5243:5342, 
#   "Transitional Basalt" = 5343:5442, "Wehrlite" = 5443:5542)
# 
# for (i in 1:num_iterations) 
# {
#   row_num <- unlist(lapply(rock_row_ranges, function(row) sample(row, length(row), replace = TRUE)))
#   pca_500 <- prcomp(data_for_pca[row_num, ], scale. = TRUE)
#   
#   #should we compare with Eu specifically ? 
#   for (j in 1:num_PCs) 
#   {
#     if (cor(pca_other_eu$rotation[,j], pca_500$rotation[,j]) < 0) 
#     {
#       pca_500$x[,j] <- -pca_500$x[,j]
#     }
#     mean_val_500[, j, i] <- sapply(rock_row_ranges, function(row) mean(pca_500$x[row, j]))
#     sd_val_500[, j, i] <- sapply(rock_row_ranges, function(row) sd(pca_500$x[row, j]))
#   }
#   
#   PC_loadings_list[[i]] <- pca_500$rotation[, 1:num_PCs]
# }
# 
# mean_low_lim <- apply(mean_val_500, c(1, 2), quantile, 0.025)
# mean_high_lim <- apply(mean_val_500, c(1, 2), quantile, 0.975)
# sd_low_lim <- apply(sd_val_500, c(1, 2), quantile, 0.025)
# sd_high_lim <- apply(sd_val_500, c(1, 2), quantile, 0.975)
# 
# 
# for (a in 1:num_PCs) 
# {
#   rock_mean <- sapply(rock_row_ranges, function(row) mean(pca_other_eu$x[row, a], na.rm = TRUE))
#   par(mar = c(10, 4, 4, 2)) 
#   plot(rock_mean, xaxt = "n", col = "red",ylim = range(c(mean_low_lim[,a], mean_high_lim[,a])),
#        main = paste("Average PC", a), ylab = paste("Mean PC", a))
#   
#   arrows(1:num_rocks, mean_low_lim[,a], 1:num_rocks, mean_high_lim[,a], length=0.03, angle=90, code=3)
#   axis(1, at = 1:num_rocks, labels = names(rock_mean), las = 2, cex.axis = 0.7)
#   grid()
# }
# 
# for (b in 1:num_PCs) 
# {
#   rock_sd <- sapply(rock_row_ranges, function(row) sd(pca_other_eu$x[row, b], na.rm = TRUE))
#   par(mar = c(10, 4, 4, 2)) 
#   plot(rock_sd, xaxt = "n", col = "blue",ylim = range(c(sd_low_lim[,b], sd_high_lim[,b])),
#        main = paste("Standard Deviation of PC", b), ylab = paste("SD of PC", b))
#   
#   arrows(1:num_rocks, sd_low_lim[,b], 1:num_rocks, sd_high_lim[,b], length=0.03, angle=90, code=3)
#   axis(1, at = 1:num_rocks, labels = names(rock_sd), las = 2, cex.axis = 0.7)
#   grid()
# }

Week 12

Rock Scatter Plot

rock_name_files <- list.files(output_directory, pattern = "\\.csv$", full.names = TRUE)
rock_data_list <- lapply(rock_name_files, function(read) read.csv(read, fileEncoding = "latin1"))
rock_names <- toupper(sub("\\.csv$", "", basename(rock_name_files), ignore.case = TRUE))
rock_rows <- unlist(lapply(rock_data_list, nrow))
rock_row_ranges <- split(1:sum(rock_rows), rep(rock_names, times = rock_rows))
num_rock_groups <- length(rock_row_ranges)

for (i in 1:15) 
{
  rock_mean_PC <- sapply(rock_row_ranges, function(row) mean(pca_other_eu$x[row, i], na.rm = TRUE))
  rock_sd_PC   <- sapply(rock_row_ranges, function(row) sd(pca_other_eu$x[row, i], na.rm = TRUE))
  
  sorted_indices <- order(rock_mean_PC)
  sorted_means   <- rock_mean_PC[sorted_indices]
  sorted_sds     <- rock_sd_PC[sorted_indices]
  sorted_names   <- names(sorted_means) 
  
  par(mar = c(10, 4, 4, 2)) 
  y_lims <- range(c(sorted_means - sorted_sds, sorted_means + sorted_sds), na.rm = TRUE)
  
  plot(sorted_means,xaxt = "n",main = paste("Average PC",i, "values with sd"),ylab = paste("Mean PC",i,"with sd"), col = "red",pch = 19, ylim = y_lims,xlim = c(1, num_rock_groups), xlab = "")

  arrows(1:num_rock_groups, sorted_means - sorted_sds,1:num_rock_groups, sorted_means + sorted_sds,code = 3, angle = 90, length = 0.05, col = "blue")
  
  axis(1,at = 1:num_rock_groups,labels = sorted_names,las = 2,cex.axis = 0.5,gap.axis = -1)
  grid()
}

Flipping Problem and PC Plots

atomic_number <- c(57,58,59,60,62,63,64,65,66,67,68,69,70,71)
ionic_radius8 <- c(1.160, 1.143, 1.126, 1.109, 1.079, 1.066, 1.053, 1.040, 1.027, 1.015, 1.004, 0.994, 0.985, 0.977)
ionic_radius6 <- c(1.032, 1.010, 0.990, 0.983, 0.958, 0.947, 0.938, 0.923, 0.912, 0.901, 0.890, 0.880, 0.868, 0.861)
pca_list <- list()
num_iterations <- 500
clean_names <- sub("\\.[A-Z]+\\.?$", "", row.names(pca_other_eu$rotation))
ree_indices <- which(clean_names != "Other")
element_symbol <- clean_names_no_other
random_pca_indices <- sample(1:num_iterations, 4)
flipped_examples <- list()

for (i in 1:num_iterations) 
{
  samples <- ree_transformed_no_other[sample(nrow(ree_transformed_no_other), replace = TRUE), ]
  pca_500_samples <- prcomp(samples)
  rotation <- pca_500_samples$rotation[ree_indices, 1:5]
  if (i %in% random_pca_indices) 
  {
    flipped_examples[[as.character(i)]] <- rotation
  }
  # check that Lu loading on all principle components is +ve
  flip <- ifelse(rotation["LU.PPM.",] < 0,-1,1)
  rotation <- sweep(rotation,2,flip,"*")
  pca_list[[i]] <- rotation
  
}

for (pc_no in 1:4)
{
   # One big plot per PC with 4 subplots
  par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
  for (iteration in names(flipped_examples)) 
  {
    pc_val <- flipped_examples[[iteration]][, pc_no]
    plot(atomic_number, pc_val, col = ifelse(atomic_number %% 2 == 0, "blue", "black"), 
         pch = 16, main = paste("PC", pc_no, "Iteration", iteration), 
         xlab = "Atomic Number", ylab = paste("Mean PC", pc_no))
    text(atomic_number, pc_val, labels = element_symbol, pos = 4, col = "red", cex = 0.6)
  }
}

pca_array <- array(unlist(pca_list), dim = c(length(ree_indices), 5, num_iterations))
pc_mean <- apply(pca_array, c(1, 2), mean)
pc_sd <- apply(pca_array, c(1, 2), sd)
x_axis_list <- list(atomic_number, ionic_radius8, ionic_radius6)
x_axis_labels <- c("Atomic Number", "Ionic Radius (CN=8)", "Ionic Radius (CN=6)")

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

for (graph_type in 1:3) 
{
  current_graph <- x_axis_list[[graph_type]]
  label_name <- x_axis_labels[graph_type]
  
  for (pc_num in 1:5) 
    {
    y_mean <- pc_mean[, pc_num]
    y_sd <- pc_sd[, pc_num]
    point_colour_coding <- ifelse(atomic_number %% 2 == 0, "blue", "black")

    plot(current_graph, y_mean, col = point_colour_coding, pch = 16,main = paste("", label_name, "vs PC", pc_num),xlab = label_name, ylab = paste("Mean PC", pc_num),ylim = range(c(y_mean - y_sd,y_mean + y_sd)),xaxt = "n")
    
    arrows(x0 = current_graph, y0 = y_mean - y_sd,x1 = current_graph, y1= y_mean + y_sd,code = 3,angle = 90, length = 0.05, col = "darkorange")
    
    text(current_graph, y_mean, labels = element_symbol, pos = 4, col = "red", cex = 0.6)
    axis(side = 1, at = current_graph, labels = round(current_graph, 3), las = 2, cex.axis = 0.7)
    legend("bottomright", legend = c("Even", "Odd"), col = c("blue", "black"),pch = 16, title = "Atomic Num", cex = 0.7)
  }
}

Week 13

Scatter Plot vs Silica Unsorted

silica_column_name <- "SIO2.WT.." 

rock_silica_mean <- sapply(rock_row_ranges, function(row) mean(combined_data[row, silica_column_name], na.rm = TRUE))

for (i in 1:15) 
{
  rock_mean_PC <- sapply(rock_row_ranges, function(row) mean(pca_other_eu$x[row, i], na.rm = TRUE))
  rock_sd_PC   <- sapply(rock_row_ranges, function(row) sd(pca_other_eu$x[row, i], na.rm = TRUE))
  
  sorted_indices <- order(rock_mean_PC)
  sorted_means   <- rock_mean_PC[sorted_indices]
  sorted_sds     <- rock_sd_PC[sorted_indices]
  sorted_names   <- names(sorted_means) 
  sorted_silica  <- rock_silica_mean[sorted_indices]
  
  par(mar = c(10, 4, 4, 2)) 
  y_lims <- range(c(sorted_means - sorted_sds, sorted_means + sorted_sds), na.rm = TRUE)
  
  plot(sorted_means, xaxt = "n",ylab = paste("Mean PC", i, "with sd"), col = "red", pch = 19, ylim = y_lims, xlim = c(1, num_rock_groups), xlab = "")

  arrows(1:num_rock_groups, sorted_means - sorted_sds, 1:num_rock_groups,sorted_means + sorted_sds, code = 3, angle = 90, length = 0.05, col = "blue")

  #rock name axis
  axis(1, at = 1:num_rock_groups, labels = sorted_names, las = 2, cex.axis = 0.5, gap.axis = -1)
  
  #silica percentage axis
  axis(side = 3, at = 1:num_rock_groups, labels = round(sorted_silica, 1),las = 2, cex.axis = 0.5, col.axis = "darkorange")
  mtext("Average SiO2 (WT%)", side = 3, line = 3, cex = 0.8, col = "black")

  grid()
}

Week 14

Scatter Plot vs Silica Sorted

silica_column_name <- "SIO2.WT.." 

rock_silica_mean <- sapply(rock_row_ranges, function(row) mean(combined_data[row, silica_column_name], na.rm = TRUE))

for (i in 1:15) 
{
  rock_mean_PC <- sapply(rock_row_ranges, function(row) mean(pca_other_eu$x[row, i], na.rm = TRUE))
  rock_sd_PC   <- sapply(rock_row_ranges, function(row) sd(pca_other_eu$x[row, i], na.rm = TRUE))
  
  sorted_indices <- order(rock_silica_mean)
  sorted_means   <- rock_mean_PC[sorted_indices]
  sorted_sds     <- rock_sd_PC[sorted_indices]
  sorted_names   <- names(sorted_means) 
  sorted_silica  <- rock_silica_mean[sorted_indices]
  
  par(mar = c(10, 4, 4, 2)) 
  y_lims <- range(c(sorted_means - sorted_sds, sorted_means + sorted_sds), na.rm = TRUE)
  
  plot(sorted_means, xaxt = "n",ylab = paste("Mean PC", i, "with sd"), col = "red", pch = 19, ylim = y_lims, xlim = c(1, num_rock_groups), xlab = "")

  arrows(1:num_rock_groups, sorted_means - sorted_sds, 1:num_rock_groups,sorted_means + sorted_sds, code = 3, angle = 90, length = 0.05, col = "blue")

  #rock name axis
  axis(1, at = 1:num_rock_groups, labels = sorted_names, las = 2, cex.axis = 0.5, gap.axis = -1)
  
  #silica percentage axis
  axis(side = 3, at = 1:num_rock_groups, labels = round(sorted_silica, 1),las = 2, cex.axis = 0.5, col.axis = "darkorange")
  mtext("Average SiO2 (WT%)", side = 3, line = 3, cex = 0.8, col = "black")

  grid()
}

Week 15

Calculate 4 Vectors

symbol_header <- c("La", "Ce", "Pr", "Nd", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Other")
n<-length(symbol_header)
ree_mid_point <- (max(atomic_number)+min(atomic_number)) /2 
standardize_vectors <- function(v)
{
  v_centered <- v - mean(v)
  v_final <- v_centered / sqrt(sum(v_centered^2))
  return(v_final)
}
# Vector 1
# Other = 1, Eu=0, REE(x13)=0
v1_initial <- c(rep(-1,n-1),1)
vector_1 <- standardize_vectors(v1_initial) * -1

# Vector 2
# Other = 0, Eu=1, REE(x13)=0

v2_initial <- c(rep(0,n))
v2_initial[symbol_header == "Eu"] <- 1
v2_no_EU <- (symbol_header != "Eu" & symbol_header != "Other")
v2_initial[v2_no_EU]<- -1/(n-2)
vector_2 <- standardize_vectors(v2_initial)

# Vector 3
# LREE vs HREE represents PC1 

v3_initial <- rep(0,n)
v3_initial[1:14] <- atomic_number
v3_initial[15]<- mean(v3_initial[1:14])
vector_3 <- standardize_vectors(v3_initial)

# Vector 4
# Quadratic LREE vs HREE represents PC2 

v4_initial <- rep(0,n)
v4_initial[1:14] <- (atomic_number-ree_mid_point)^2
v4_initial[15] <- mean(v4_initial[1:14])
vector_4 <- standardize_vectors(v4_initial)

combined_vectors <- cbind(vector_1,vector_2,vector_3,vector_4)
rownames(combined_vectors) <- symbol_header
colnames(combined_vectors) <- c("Vector_1","Eu_Anomaly","Linear_LREEvsHREE","Quad_LREEvsHREE")
#print(combined_vectors)

Solving Vectors

x_project <- function(vectors,transformeddata)
{
  A <- as.matrix(vectors)
  B <- as.matrix(transformeddata)
  scores <- B %*% A
  return(as.data.frame(scores))
}

projected_scores <- x_project(combined_vectors,ree_transformed)
colnames(projected_scores) <- c("Total_REE", "Eu_Anomaly", "Slope", "Quad")
#print(projected_scores)
reconstructed_data <- as.matrix(projected_scores) %*% MASS::ginv(combined_vectors)
residual_epsilon <- reconstructed_data - as.matrix(ree_transformed)

Reading Week 2

Adding rock name col

cleanest_rock_names <- as.character(non_ree_data$ROCK.NAME)
cleanest_rock_names <- sapply(strsplit(cleanest_rock_names, "[/,]"), `[`, 1)
cleanest_rock_names <- trimws(toupper(cleanest_rock_names))
cleanest_rock_names <- gsub(" ", "_",cleanest_rock_names)                

projected_scores$Rock_Type <- cleanest_rock_names
clean_sorted_names <- gsub(" ", "_", trimws(toupper(sorted_names)))

projected_scores_filtered <- projected_scores[projected_scores$Rock_Type %in% clean_sorted_names, ]

rock_mean_scores <- aggregate(cbind(Total_REE, Eu_Anomaly, Slope, Quad) ~ Rock_Type,data = projected_scores_filtered,FUN = mean)
rock_sd_scores <- aggregate(cbind(Total_REE, Eu_Anomaly, Slope, Quad) ~ Rock_Type, 
                            data = projected_scores_filtered, FUN = sd)

Rock scatter Plots

axes_pair <- list(c("Total_REE", "Eu_Anomaly"),c("Total_REE", "Slope"),c("Total_REE", "Quad"),c("Eu_Anomaly", "Slope"),c("Eu_Anomaly", "Quad"),c("Slope", "Quad"))



for(i in 1:length(axes_pair)) 
{
  x_n <- axes_pair[[i]][1]
  y_n <- axes_pair[[i]][2]
  
  x_val <- rock_mean_scores[[x_n]]
  y_val <- rock_mean_scores[[y_n]]
  
  #x_limit <- 
  #y_limit <- 
  plot(x_val, y_val,xlab = x_n, ylab = y_n,main = paste(x_n, "vs", y_n),pch = 19, cex = 0.5)
  text(x_val, y_val,labels = rock_mean_scores$Rock_Type,pos = 3, offset =0.1,cex = 0.5)
}

Week 16

chosen rocks

chosen_rocks <- data.frame(
Rock_Type = 
c(
    "ADAKITE", "ALKALI_BASALT", "ANDESITE", "BASALT", "CARBONATITE", 
    "DIORITE", "DOLERITE", "DUNITE", "GABBRO", "GRANITE", 
    "HARZBURGITE", "KIMBERLITE", "KOMATIITE", "LAMPROPHYRE", "LATITE", 
    "ORE", "ORTHOPYROXENITE", "PERIDOTITE", "PICRITE", "PYROXENITE", 
    "RHYOLITE", "TRACHYTE", "TRANSITIONAL_BASALT", "WEHRLITE"),
  Short_Name = c(
    "Ada", "AlB", "And", "Bas", "Carb", 
    "Dio", "Dol", "Dun", "Gab", "Grn", 
    "Hzb", "Kim", "Kom", "Lmp", "Lat", 
    "Ore", "Opx", "Per", "Pic", "Pxn", 
    "Rhy", "Tra", "TBas", "Whr")
)

rock_mean_scores_filtered <- merge(rock_mean_scores, chosen_rocks, by = "Rock_Type")
rock_sd_filtered <- merge(rock_sd_scores, chosen_rocks, by = "Rock_Type")

plotting chosen rocks

for(i in 1:length(axes_pair)) 
{
  x_n <- axes_pair[[i]][1]
  y_n <- axes_pair[[i]][2]
  
  x_val <- rock_mean_scores_filtered[[x_n]]
  y_val <- rock_mean_scores_filtered[[y_n]]
  
  x_sd <- rock_sd_filtered[[x_n]]
  y_sd <- rock_sd_filtered[[y_n]]
  
  x_limit <- c(min(x_val - x_sd, na.rm=T), max(x_val + x_sd, na.rm=T))
  y_limit <- c(min(y_val - y_sd, na.rm=T), max(y_val + y_sd, na.rm=T))
  
  plot(x_val, y_val,xlab = x_n, ylab = y_n,main = paste(x_n, "vs", y_n),pch = 19, cex = 0.5, xlim = x_limit, ylim = y_limit)
  
  #error bars
  #vertical
  arrows(x_val, y_val - y_sd, x_val, y_val + y_sd, length = 0.05, angle = 90, code = 3, col = "gray")
  #horizontal
  arrows(x_val - x_sd, y_val, x_val + x_sd, y_val, length = 0.05, angle = 90, code = 3, col = "gray")
  
  text(x_val, y_val,labels = rock_mean_scores_filtered$Short_Name,pos = 3, offset =0.1,cex = 0.5, col = "red")
}

residual plot

plot(as.matrix(ree_transformed), reconstructed_data,xlab="True Values (ree_transformed)", ylab="Predicted Values (reconstructed)",main="Fit vs True Values", pch=19, cex=0.5, col=rgb(0,0,0,0.2))

abline(0, 1, col="green", lwd=2, lty=2)

Week 17

CHUR

# CHUR values from McDonough, W. F., & Sun, S. S. (1995). The composition of the Earth.
chur_ppm <- c(0.237, 0.613, 0.092, 0.457, 0.148, 0.056, 0.199, 0.036, 0.249, 0.054, 0.160, 0.025, 0.161, 0.024)

chur_with_other <- c(chur_ppm, 1e6 - sum(chur_ppm))
# clr transformation 
chur_clr <- log(chur_with_other / exp(mean(log(chur_with_other))))
ree_minus_chur <- sweep(as.matrix(ree_transformed), 2, chur_clr, "-")
# projected scores with chur
projected_scores_chur <- as.data.frame(ree_minus_chur %*% combined_vectors)
colnames(projected_scores_chur) <- c("Total_REE", "Eu_Anomaly", "Slope", "Quad")
# make total_ree +ve 
projected_scores_chur$Total_REE <- projected_scores_chur$Total_REE 

projected_scores_chur$Rock_Type <- cleanest_rock_names
projected_scores_chur_filtered <- projected_scores_chur[projected_scores_chur$Rock_Type %in% chosen_rocks$Rock_Type, ]

# mean and sd
rock_mean_chur <- aggregate(cbind(Total_REE, Eu_Anomaly, Slope, Quad) ~ Rock_Type, data = projected_scores_chur_filtered, FUN = mean)
rock_sd_chur <- aggregate(cbind(Total_REE, Eu_Anomaly, Slope, Quad) ~ Rock_Type, data = projected_scores_chur_filtered, FUN = sd)

plot_data_chur <- merge(rock_mean_chur, chosen_rocks, by = "Rock_Type")
sd_data_chur <- merge(rock_sd_chur, chosen_rocks, by = "Rock_Type")

# chur as origin (0,0)
chur_ref <- data.frame(Rock_Type="CHUR", Total_REE=0, Eu_Anomaly=0, Slope=0, Quad=0, Short_Name="CHUR")
plot_data_chur <- rbind(plot_data_chur, chur_ref)
sd_data_chur <- rbind(sd_data_chur, chur_ref)

for(i in 1:length(axes_pair)) 
{
  x_n <- axes_pair[[i]][1]
  y_n <- axes_pair[[i]][2]
  
  x_val <- plot_data_chur[[x_n]]
  y_val <- plot_data_chur[[y_n]]
  x_sd  <- sd_data_chur[[x_n]]
  y_sd  <- sd_data_chur[[y_n]]
  
  x_limit <- range(c(x_val - x_sd, x_val + x_sd, 0), na.rm=T)
  y_limit <- range(c(y_val - y_sd, y_val + y_sd, 0), na.rm=T)
  
  plot(x_val, y_val, xlab = paste(x_n, "(vs CHUR)"), ylab = paste(y_n, "(vs CHUR)"),main = paste(x_n, "vs", y_n, "Normalized to CHUR"),pch = 19, cex = 0.6, xlim = x_limit, ylim = y_limit,col = ifelse(plot_data_chur$Short_Name == "CHUR", "blue", "black"), asp=1)
  
  abline(h = 0, v = 0, lty = 3, col = "gray60")
  
  # error bars showing warnings ?
suppressWarnings(arrows(x_val, y_val - y_sd, x_val, y_val + y_sd, length = 0.04, angle = 90, code = 3, col = "gray70"))
suppressWarnings(arrows(x_val - x_sd, y_val, x_val + x_sd, y_val, length = 0.04, angle = 90, code = 3, col = "gray70"))
  
  text(x_val, y_val, labels = plot_data_chur$Short_Name, pos = 3, offset = 0.2, cex = 0.6, col = ifelse(plot_data_chur$Short_Name == "CHUR", "blue", "red"))
}

Chur residual

reconstructed_data_chur <- as.matrix(projected_scores_chur[, 1:4]) %*% MASS::ginv(combined_vectors)

plot(as.matrix(ree_minus_chur), reconstructed_data_chur,xlab="True Values (ree_minus_chur)", ylab="Predicted Values (reconstructed_data-chur)",main="Fit vs True Values", pch=19, cex=0.5, col=rgb(0,0,0,0.2))

abline(0, 1, col="green", lwd=2, lty=2)

Week 18

Chur Residual PCA

res_mat <- as.matrix(ree_minus_chur[, 1:14]) - reconstructed_data_chur[, 1:14]
colnames(res_mat) <- clean_names[1:14]
pca_chur_residual <- prcomp(res_mat, scale. = TRUE)
biplot(pca_chur_residual, xlabs = rep(".", nrow(pca_chur_residual$x)))

PC1 : Odd ve even

chur_residual_matrix <- as.matrix(ree_minus_chur[, 1:14]) - reconstructed_data_chur[, 1:14]

par(mfrow = c(4, 4), mar = c(2, 2, 2, 1))
invisible(sapply(1:14, function(i) hist(chur_residual_matrix[, i], main = colnames(chur_residual_matrix)[i],col = "blue", breaks = 30, xlab = "")))

Week 19

CHUR residual scatter plots

point_col <- ifelse(atomic_number %% 2 == 0, "blue", "black")
par(mar = c(6, 4, 4, 2))
for (i in 1:3) 
{
  plot(x_axis_list[[i]], pca_chur_residual$rotation[,1], col = point_col, pch = 16, xlab = x_axis_labels[i], ylab = "Residual PC1 Loadings",main = paste("Residual PC1 vs", x_axis_labels[i]),xaxt = "n")
  text(x_axis_list[[i]], pca_chur_residual$rotation[,1], labels = element_symbol, pos = 4, col = "red", cex = 0.7)
  axis(side = 1, at = x_axis_list[[i]], labels = round(x_axis_list[[i]], 3), las = 2, cex.axis = 0.7)
  abline(h = 0, lty = 2, col = "gray60")
  legend("bottomright", legend = c("Even", "Odd"), col = c("blue", "black"), pch = 16, cex = 0.4)
}

More ref points

ref_ree_list <- list(
  CHUR_1= c(0.237,0.612,0.095,0.4670,0.1530,0.0580,0.2055,0.0374,0.2540,0.0566,0.1655,0.0255,0.17,0.0254), # avg CI Chondrite Sun and McDonough-1989
  CHUR_2= c(0.2469,0.6321,0.0959,0.4864,0.1556,0.0599,0.2003,0.0378,0.2577,0.0554,0.1667,0.0261,0.1694,0.0256), # avg CI Chondrite Pourmand et al 2012
  CHUR_3  = c(0.2530,0.6,0.0910,0.4640,0.1530,0.0586,0.2060,0.0375,0.2540,0.0566,0.166,0.0262,0.168,0.0246), # avg CI Chondrite Barrat et al 2012
  CHUR_4 = c(0.3670,0.9570,0.1370,0.711,0.231,0.087,0.306,0.058,0.381,0.0851,0.2490,0.0356,0.248,0.0381), # CI Chondrite volatile free Taylor and Mclennan 2009
  CHUR_5 = c(0.3191,0.8201,0.1211,0.6151,0.2,0.0761,0.2673,0.0494,0.33,0.756,0.216,0.0329,0.2209,0.033), # CI Chondrite volatile free korotev 1996
  NASC = c(32,73,7.9,33,5.7,1.24,5.2,0.85,5.8,1.04,3.4,0.5,3.1,0.48),# Mclennan 1989
  UCC = c(31,63,7.1,27,4.7,1,4,0.7,3.9,0.83,2.3,0.3,2,0.31) # rudnick and gao 2014
)

ref_matrix <- do.call(rbind, ref_ree_list)
ref_matrix_all <- cbind(ref_matrix, 1e6 - rowSums(ref_matrix))

#  CLR and Normalize to original CHUR
ref_clr_mat <- log(ref_matrix_all / exp(rowMeans(log(ref_matrix_all))))
ref_shifted <- sweep(ref_clr_mat, 2, chur_clr, "-") # chur_clr is the origin

ref_coords <- as.data.frame(ref_shifted %*% combined_vectors)
colnames(ref_coords) <- c("Total_REE", "Eu_Anomaly", "Slope", "Quad")
ref_coords$Short_Name <- names(ref_ree_list)

#Merge the new references

ref_coords$Rock_Type <- ref_coords$Short_Name 
plot_data_chur <- rbind(plot_data_chur, ref_coords)

#no sd in ref points
ref_sd_null <- ref_coords
ref_sd_null[, 1:4] <- 0
sd_data_chur <- rbind(sd_data_chur, ref_sd_null)

#point names
all_refs <- c("CHUR", names(ref_ree_list))

#updated loop
for(i in 1:length(axes_pair)) 
{
  x_n <- axes_pair[[i]][1]
  y_n <- axes_pair[[i]][2]
  
  x_val <- plot_data_chur[[x_n]]
  y_val <- plot_data_chur[[y_n]]
  x_sd  <- sd_data_chur[[x_n]]
  y_sd  <- sd_data_chur[[y_n]]
  
  x_limit <- range(c(x_val - x_sd, x_val + x_sd, 0), na.rm=T)
  y_limit <- range(c(y_val - y_sd, y_val + y_sd, 0), na.rm=T)
  
  # makes ref points blue
  plot(x_val, y_val, 
       xlab = paste(x_n, "(vs CHUR)"), ylab = paste(y_n, "(vs CHUR)"),
       main = paste(x_n, "vs", y_n, "Normalized to CHUR"),
       pch = 19, cex = 0.6, xlim = x_limit, ylim = y_limit,
       col = ifelse(plot_data_chur$Short_Name %in% all_refs, "blue", "black"), 
       asp=1)
  
  abline(h = 0, v = 0, lty = 3, col = "gray60")
  
  
  suppressWarnings(arrows(x_val, y_val - y_sd, x_val, y_val + y_sd, length = 0.04, angle = 90, code = 3, col = "gray70"))
  suppressWarnings(arrows(x_val - x_sd, y_val, x_val + x_sd, y_val, length = 0.04, angle = 90, code = 3, col = "gray70"))
  
  text(x_val, y_val, labels = plot_data_chur$Short_Name, pos = 3, offset = 0.2, cex = 0.6, 
       col = ifelse(plot_data_chur$Short_Name %in% all_refs, "blue", "red"))
}

spider_samples <- list(ree_data[249,1:14],ree_data[434,1:14],ree_data[1158,1:14],ree_data[1973,1:14],ree_data[3616,1:14],ree_data[4278,1:14])
spider_ref <- c(0.237,0.612,0.095,0.4670,0.1530,0.0580,0.2055,0.0374,0.2540,0.0566,0.1655,0.0255,0.17,0.0254)
rock_names_graph <- c("Andesite","Basalt","Carbonatite","Granite","Peridotite","Rhyolite")

plot_spider(spider_samples,spider_ref)

# 1. Create the plot but suppress the X axis (xaxt = "n")
plot(1:14, as.numeric(ree_data[1, 1:14]), 
     type = "b",          # Connects points with lines
     pch = 19,            # Solid circles
     col = "black", 
     main = "REE Concentration in Adakite Sample",
     xlab = "Rare Earth Elements",   # Overall title for X axis
     ylab = "Concentration (PPM)",
     xaxt = "n")          # "n" tells R not to draw the X axis yet

# 2. Add the custom labels (clean_names_no_other) to the X axis
axis(side = 1,            # 1 = bottom axis
     at = 1:14,           # Where to put the labels
     labels = clean_names_no_other, # The actual element names
     las = 1,             # Optional: makes labels vertical if they overlap
     cex.axis = 0.8)      # Optional: slightly smaller text if names are crowded