Une mise-à-jour de Gitlab vers la version 14.0 est prévue le 24 juin entre 14:00 et 15:00. Le service sera inaccessible ou instable pendant cette période. Merci de votre compréhension.

prg1.R 3.55 KB
Newer Older
Dumoulin Nicolas's avatar
Dumoulin Nicolas committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98


setwd("D://workspace/StockedDiad/R")
library(pls)
library(plsdepot)

mydata <- read.table("D://workspace/StockedDiad/data/output/stockedDiadOutput5final.txt", sep=",", header=TRUE)

mydata$recovery = mydata$wildActiveFemaleSpawner>5
mydata1=subset(mydata, mydata$lastFailInReproduction<49.25)

# output="geoMeanEffective"
output="lastFailInReproduction"
 # output="meanWildActiveFemaleSpawner"
#output = "recovery"
boxplot(mydata[, output] ~ mydata$probabilityToBirthHabitat, xlab='prob to birth habitat', ylab=output)
boxplot(mydata[, output] ~ mydata$probabilityToBirthRiver, xlab='prob to birth river')
boxplot(mydata1[, output] ~ mydata1$releaseAmount, xlab='release amount')
mydata$Mafter
table(mydata[, c("recovery","femaleMeanAgeAtMaturation" )])
table(mydata[, c("recovery","Mafter")])
# ======================================================================================
# pls avec pls
summary(glm(as.formula(paste(output,"~ as.factor(probabilityToBirthHabitat) + as.factor(probabilityToBirthRiver)", sep="")), data=mydata))

pls1=plsr(lastFailInReproduction ~ as.factor(M0_3months)+as.factor(M3_12months)+ 
            as.factor(M1year)+ as.factor(Mafter)+
       as.factor(femaleMeanAgeAtMaturation) + as.factor(femaleSdAgeAtMaturation) +
       as.factor(fecundity)+
       as.factor(femaleMeanInterSpawningInterval) + 
       as.factor(maleMeanAgeAtMaturation) + as.factor(maleSdAgeAtMaturation) +
       as.factor(maleMeanInterSpawningInterval)+
       as.factor(probabilityToBirthRiver)+ as.factor(probabilityToBirthHabitat) +
       as.factor(releaseAmount) + as.factor(sigmaRecruitment),     
     ncomp = 10, data = mydata1, method = "oscorespls")

summary(pls1)

vip = VIP(pls1)
barplot(vip[1,order(vip[1,], decreasing = TRUE)])
plot(pls1)
rm(pls2)
pls2=plsr(lastFailInReproduction ~ M0_3months + M3_12months + 
            M1year+ Mafter +
            femaleMeanAgeAtMaturation + femaleSdAgeAtMaturation +
            fecundity +
            femaleMeanInterSpawningInterval + 
            maleMeanAgeAtMaturation + maleSdAgeAtMaturation +
            maleMeanInterSpawningInterval +
            probabilityToBirthRiver + probabilityToBirthHabitat +
            releaseAmount + sigmaRecruitment,     
          ncomp = 10, data = mydata1, method = "oscorespls", model=TRUE, validation ="CV")
pls2$model
summary(pls2)
plot(RMSEP(pls2))
explvar(pls2)
RMSEP(pls2)
pls2$loading.weights
pls2$Yloadings
VIP(pls2)

VIPjh(pls2,1,4)

names(vip2[1,])
names(vip2[1,])
barplot(vip2[1, order(vip2[1,], decreasing = TRUE)])
names(vip2[1, order(vip2[1,], decreasing = TRUE)])

vip2[1,]
vip2[1,order(vip2[1,])]
# ===============================================================================
# pls avec plsdepot


sum(mydata$lastFailInReproduction<49.25)
hist(mydata1$lastFailInReproduction)

Xvariable = c("M0_3months", "M3_12months", "M1year", "Mafter", 
              "femaleMeanAgeAtMaturation", "femaleSdAgeAtMaturation", "fecundity",
              "femaleMeanInterSpawningInterval",
              "maleMeanAgeAtMaturation", "maleSdAgeAtMaturation", "maleMeanInterSpawningInterval",
              "probabilityToBirthRiver", "probabilityToBirthHabitat",
              "releaseAmount", "sigmaRecruitment" )
myplsreg1 = plsreg1(mydata1[,Xvariable], mydata1$lastFailInReproduction, comps=3)

myplsreg1$x.loads
myplsreg1$R2Xy

myplsreg2$
myplsreg1$R2
myplsreg1$R2Xy
head(myplsreg1$T2)
summary(myplsreg1)
plot(myplsreg1)
plot(mydata1$lastFailInReproduction, myplsreg1$y.pred, xlab = "Original",
     ylab = "Predicted")

data(oliveoil)