In martinEcarnot/signe: Useful functions to handle spectra

Load libraries

library(R.matlab)
library(prospectr)
library(nirsextra)
library(MetStaT)
source("Script_R_2020/sp2dfclo.R")

Load globalmatix

# brb3="/home/ecarnot/Documents/INRA/Projets/SIGNE/2019/globalmatrix"
brb3="globalmatrix"
load(file=brb3)
sp=globalmatrix

Set parameters

# Savitsky-Golay parameters (p=degree of polynom, n= window size, m=derivation order)
p=2
n=11
m=1

Pretreatments

## Remove gaps from 3 sub-spectrometers (Montpellier: 1000nm (#651) et 1800nm (#1451))
sp_pre=adj_asd(sp,c(276,676)) #Retrouvé empiriquement, ne correspondait pas à ce qui etait indiqué precedemment.

## SNV
sp_pre=t(scale(t(sp_pre)))

##Savitsky Golay Derivation
sp=savitzkyGolay(sp_pre, m = m, p = p, w = n)

Creaction of the factors

class=as.factor(substr(rownames(sp),9,9)) #9,9
classclo=as.factor(substr(rownames(sp),9,13)) #9,13

## Variable qui mesure le nombre de classes
c=length(levels(class))
cclo=length(levels(classclo))


# On créé un facteur datclone qui groupe un clone à 1 date
datclone=as.factor(substr(rownames(sp),1,13))
ndc=length(unique(datclone))
# On créé un facteur souche qui groupe les 6 spectres de chaque souche
numsp=as.numeric(substr(rownames(sp),15,16))
souche=cut(numsp, breaks = c(0,6,12,18),labels=c("s1","s2","s3"))
endroit=as.numeric(substr(rownames(sp),15,16))%%3
rownames(sp)=paste(rownames(sp),souche,endroit)

sp=sp2dfclo(sp,class,classclo)
sp=cbind(sp,datclone,souche)

colnames(sp)[c(1,2)]=c("cep","clo")
sp$annee=as.factor(substr(rownames(sp),1,4))
sp$parcelle=as.factor(substr(rownames(sp),18,18))
sp$position=as.factor(substr(rownames(sp),23,23))
sp$jour=as.factor(substr(rownames(sp),1,8))
sp$mois=as.factor(substr(rownames(sp),5,6))

Filter data + Save data set for Matlab

set_name=c("all","parcG","cepG","parcGcepG","parcGcepC","parcGcepS")

for (i in 1:length(set_name)) {
  if (set_name[i]=="all") {
    ired=sample(1:nrow(sp),nrow(sp)/2)  # Too big : reduse by 2
    spf=sp[ired,]
    name_fact=c("cep","clo","annee","parcelle","mois")
  }
  else if (set_name[i]=="parcG") {
    spf=droplevels(sp[which(sp$parcelle=="G"),])
    name_fact=c("cep","clo","annee","mois")
  }
  else if (set_name[i]=="cepG") {
    spf=droplevels(sp[which(sp$cep=="G"),])
    name_fact=c("clo","annee","parcelle","mois")
  }
  else if (set_name[i]=="parcGcepG") {
    spf=droplevels(sp[which(sp$parcelle=="G" & sp$cep=="G"),])
    name_fact=c("clo","annee","mois")
  }
  else if (set_name[i]=="parcGcepC") {
    spf=droplevels(sp[which(sp$parcelle=="G" & sp$cep=="C"),])
    name_fact=c("clo","annee","mois")
  }
  else if (set_name[i]=="parcGcepS") {
    spf=droplevels(sp[which(sp$parcelle=="G" & sp$cep=="S"),])
    name_fact=c("clo","annee","mois")
  }

  # Factors to numeric (for ASCA function)
  id_x=match(name_fact, names(spf))
  fact=as.matrix(sapply(spf[,id_x], as.numeric))

  writeMat(paste0("asca_save/SIGNE_asca_",set_name[i],".mat"), x=spf$x, fact=fact, name_fact=name_fact)
}

Load Matlab result for plot

ascaMat=readMat('asca_save/SIGNE_asca_ev.mat')

#All
# dat=data.frame(t(ascaMat$ev.all[[1]][[1]]))
# colnames(dat)=c("Residal",unlist(ascaMat$fact.all[[1]][[1]]))
# barplot(as.matrix(dat[,-1]),main=ascaMat$set.name[[1]][[1]])

# Parcelle G (Grau du Roi, parcelle 96), all cep.
dat=data.frame(t(ascaMat$ev.all[[2]][[1]]))
colnames(dat)=c("Residal",unlist(ascaMat$fact.all[[2]][[1]]))
barplot(as.matrix(dat[,-1]),main=ascaMat$set.name[[2]][[1]])

# Parcelle G (Grau du Roi, parcelle 96), intra-cep variance
dat=t(data.frame(ascaMat$ev.all[[4]][[1]],ascaMat$ev.all[[5]][[1]],ascaMat$ev.all[[6]][[1]]))
rownames(dat)=substr(unlist(lapply(ascaMat$set.name,function (x) {x[[]][[1]]}))[4:6],6,9)
colnames(dat)=c("Residal",unlist(ascaMat$fact.all[[5]][[1]]))
barplot(dat[,-1], beside = TRUE, legend=T,main="Parcelle G")

ASCA MetStat # Donne des resultats differents de ASCA Matlab (version )

spf_xr=spf$x[seq(1,nrow(spf$x),2),seq(1,ncol(spf$x),1)]
factr=fact[seq(1,nrow(fact),2),]

res=ASCA.Calculate(spf_xr, factr) # Donne des resultats differents de ASCA Matlab (version )
#ASCA.DoPermutationTest(res)


opt=asca("options")
opt$permtest="off"
r=asca(spf_xr, factr,opt)

ev=c(r$XA$EffectExplVar,r$XB$EffectExplVar,r$XC$EffectExplVar,r$XD$EffectExplVar,r$XE$EffectExplVar)
#p=c(r$XA$EffectSignif$p,r$XB$EffectSignif$p,r$XAB$EffectSignif$p)
disp_fac=paste( r$TermLabels[-1],c(colnames(factr), paste(colnames(factr),collapse="x")),sep = "=")
cat(disp_fac)
print(ev)
#print(p)