Commit d82ffcf8 authored by Boulangeat Isabelle's avatar Boulangeat Isabelle
Browse files

analyse hmean model

parent 8ab10cda
...@@ -159,3 +159,46 @@ ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = ea ...@@ -159,3 +159,46 @@ ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = ea
labs(y="easting", x="Type vegetation") labs(y="easting", x="Type vegetation")
``` ```
### Ensemble de variables abiotiques
```{r ,fig=TRUE}
releves_all = readRDS("releves_all_plus.rds")
library(ade4)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900", "frost_severe.gdd300", "frost_severe.gdd600",
"frost_severe.gdd900", "frost.gdd300", "frost.gdd600", "frost.gdd900",
"gddspeed.gdd300", "gddspeed.gdd600", "gddspeed.gdd900", "radiations.gdd300",
"radiations.gdd600", "radiations.gdd900", "rainfall.gdd300", "rainfall.gdd600",
"rainfall.gdd900", "snowdays.gdd300", "snowdays.gdd600", "snowdays.gdd900" )
dat = na.omit(releves_all[,vars])
dim(dat)
pca = dudi.pca(dat, scan=FALSE, nf = 4)
inertia.dudi(pca, row=FALSE, col =TRUE)
par(mfrow=c(1,2))
s.corcircle(pca$co)
s.corcircle(pca$co, xax=3, yax=4)
```
### Modélisation de la hauteur moyenne
```{r ,fig=TRUE}
dat = na.omit(releves_all[,c(vars, "hmean")])
dat$lf = as.factor(dat$lf)
dim(dat)
library(randomForest)
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd300 + rainfall.gdd900 + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, ntree = 150, mtry=16)
mod
varImpPlot(mod)
```
```{r ,fig=TRUE}
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd900 + slope + northing + easting + lf, data = dat, ntree = 100, mtry=16)
mod
library(plotmo)
plotmo(mod, ylim = NA)
```
This source diff could not be displayed because it is too large. You can view the blob instead.
...@@ -230,3 +230,317 @@ ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = ea ...@@ -230,3 +230,317 @@ ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = ea
![](README_files/figure-html/unnamed-chunk-14-1.png)<!-- --> ![](README_files/figure-html/unnamed-chunk-14-1.png)<!-- -->
### Ensemble de variables abiotiques
```r
releves_all = readRDS("releves_all_plus.rds")
```
```
## Warning in readRDS("releves_all_plus.rds"): strings not representable in native
## encoding will be translated to UTF-8
```
```r
library(ade4)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900", "frost_severe.gdd300", "frost_severe.gdd600",
"frost_severe.gdd900", "frost.gdd300", "frost.gdd600", "frost.gdd900",
"gddspeed.gdd300", "gddspeed.gdd600", "gddspeed.gdd900", "radiations.gdd300",
"radiations.gdd600", "radiations.gdd900", "rainfall.gdd300", "rainfall.gdd600",
"rainfall.gdd900", "snowdays.gdd300", "snowdays.gdd600", "snowdays.gdd900" )
dat = na.omit(releves_all[,vars])
dim(dat)
```
```
## [1] 326 29
```
```r
pca = dudi.pca(dat, scan=FALSE, nf = 4)
inertia.dudi(pca, row=FALSE, col =TRUE)
```
```
## Inertia information:
## Call: inertia.dudi(x = pca, row.inertia = FALSE, col.inertia = TRUE)
##
## Decomposition of total inertia:
## inertia cum cum(%)
## Ax1 1.214e+01 12.14 41.87
## Ax2 3.815e+00 15.96 55.03
## Ax3 2.938e+00 18.90 65.16
## Ax4 2.495e+00 21.39 73.76
## Ax5 2.155e+00 23.55 81.19
## Ax6 1.097e+00 24.64 84.98
## Ax7 9.678e-01 25.61 88.31
## Ax8 8.821e-01 26.49 91.36
## Ax9 6.313e-01 27.12 93.53
## Ax10 5.531e-01 27.68 95.44
## Ax11 3.860e-01 28.06 96.77
## Ax12 2.496e-01 28.31 97.63
## Ax13 2.234e-01 28.54 98.40
## Ax14 1.103e-01 28.65 98.78
## Ax15 1.040e-01 28.75 99.14
## Ax16 8.821e-02 28.84 99.45
## Ax17 4.356e-02 28.88 99.60
## Ax18 3.133e-02 28.91 99.70
## Ax19 2.890e-02 28.94 99.80
## Ax20 2.208e-02 28.97 99.88
## Ax21 1.484e-02 28.98 99.93
## Ax22 1.126e-02 28.99 99.97
## Ax23 5.520e-03 29.00 99.99
## Ax24 3.240e-03 29.00 100.00
## Ax25 1.676e-05 29.00 100.00
##
## Column contributions (%):
## alti slope lf northing
## 3.448 3.448 3.448 3.448
## easting snowfree_time cumgdd0 cumgdd60d
## 3.448 3.448 3.448 3.448
## albedo.gdd300 albedo.gdd600 albedo.gdd900 frost_severe.gdd300
## 3.448 3.448 3.448 3.448
## frost_severe.gdd600 frost_severe.gdd900 frost.gdd300 frost.gdd600
## 3.448 3.448 3.448 3.448
## frost.gdd900 gddspeed.gdd300 gddspeed.gdd600 gddspeed.gdd900
## 3.448 3.448 3.448 3.448
## radiations.gdd300 radiations.gdd600 radiations.gdd900 rainfall.gdd300
## 3.448 3.448 3.448 3.448
## rainfall.gdd600 rainfall.gdd900 snowdays.gdd300 snowdays.gdd600
## 3.448 3.448 3.448 3.448
## snowdays.gdd900
## 3.448
##
## Column absolute contributions (%):
## Axis1 Axis2 Axis3 Axis4
## alti 2.35689 9.71773 0.78391 1.266e+00
## slope 0.57743 3.75004 0.84670 5.896e-02
## lf 0.01174 2.01464 1.31009 1.490e-01
## northing 1.10364 0.33057 2.65148 3.062e-01
## easting 0.76428 1.59800 0.03336 7.096e-02
## snowfree_time 0.59547 9.49618 11.35769 1.095e+00
## cumgdd0 0.76251 15.99983 5.03075 6.603e-01
## cumgdd60d 0.67773 17.29617 4.14299 4.977e-01
## albedo.gdd300 7.49195 0.17806 0.71629 8.400e-01
## albedo.gdd600 7.63105 0.19890 0.07115 6.415e-01
## albedo.gdd900 5.70258 1.73350 0.56764 2.602e-01
## frost_severe.gdd300 1.64300 0.04100 1.73959 2.947e+01
## frost_severe.gdd600 1.64300 0.04100 1.73959 2.947e+01
## frost_severe.gdd900 1.59739 0.01727 1.62883 2.936e+01
## frost.gdd300 6.40042 0.68087 0.82578 4.488e-01
## frost.gdd600 6.40042 0.68087 0.82578 4.488e-01
## frost.gdd900 5.91168 1.46054 0.65286 3.809e-01
## gddspeed.gdd300 7.19982 0.38288 0.10875 4.519e-01
## gddspeed.gdd600 6.82718 0.40637 1.89308 1.875e-01
## gddspeed.gdd900 3.57737 0.21330 7.70715 4.415e-04
## radiations.gdd300 5.53633 0.21733 2.59806 8.290e-03
## radiations.gdd600 3.99008 0.74914 4.79785 1.218e-01
## radiations.gdd900 2.48107 1.44946 7.83808 6.309e-01
## rainfall.gdd300 2.03240 8.53988 5.09403 5.156e-03
## rainfall.gdd600 0.86467 7.75374 11.52810 3.083e-01
## rainfall.gdd900 0.45027 6.73033 12.71229 1.262e-01
## snowdays.gdd300 5.39061 2.36479 3.81811 9.447e-01
## snowdays.gdd600 5.39061 2.36479 3.81811 9.447e-01
## snowdays.gdd900 4.98841 3.59283 3.16190 8.390e-01
##
## Signed column relative contributions:
## Axis1 Axis2 Axis3 Axis4
## alti -28.6203 -37.07608 2.3028 3.157825
## slope -7.0119 -14.30755 -2.4872 -0.147095
## lf 0.1425 -7.68647 3.8485 -0.371675
## northing -13.4019 -1.26123 -7.7889 -0.763863
## easting 9.2808 6.09687 -0.0980 0.177033
## snowfree_time -7.2309 -36.23082 33.3638 2.732736
## cumgdd0 -9.2594 -61.04422 14.7781 1.647412
## cumgdd60d -8.2298 -65.99013 12.1703 1.241669
## albedo.gdd300 90.9769 -0.67935 -2.1041 -2.095804
## albedo.gdd600 92.6660 -0.75887 -0.2090 -1.600464
## albedo.gdd900 69.2480 -6.61381 1.6675 -0.649227
## frost_severe.gdd300 19.9514 0.15641 -5.1101 73.533219
## frost_severe.gdd600 19.9514 0.15641 -5.1101 73.533219
## frost_severe.gdd900 19.3976 0.06589 -4.7848 73.254688
## frost.gdd300 77.7221 -2.59773 -2.4258 -1.119687
## frost.gdd600 77.7221 -2.59773 -2.4258 -1.119687
## frost.gdd900 71.7872 -5.57239 -1.9178 -0.950205
## gddspeed.gdd300 87.4295 1.46080 0.3195 -1.127433
## gddspeed.gdd600 82.9044 1.55044 5.5610 -0.467794
## gddspeed.gdd900 43.4410 -0.81380 22.6402 0.001101
## radiations.gdd300 67.2292 -0.82917 7.6319 -0.020684
## radiations.gdd600 48.4527 -2.85821 14.0939 0.303766
## radiations.gdd900 30.1283 -5.53013 23.0248 1.574111
## rainfall.gdd300 24.6799 32.58222 14.9640 -0.012863
## rainfall.gdd600 10.5000 29.58286 33.8644 0.769237
## rainfall.gdd900 5.4678 25.67824 37.3431 0.314976
## snowdays.gdd300 65.4597 -9.02240 -11.2159 -2.357037
## snowdays.gdd600 65.4597 -9.02240 -11.2159 -2.357037
## snowdays.gdd900 60.5756 -13.70774 -9.2883 -2.093141
##
## Cumulative sum of column relative contributions (%):
## Axis1 Axis1:2 Axis1:3 Axis1:4 Axis5:25
## alti 28.6203 65.696 68.00 71.16 28.843
## slope 7.0119 21.319 23.81 23.95 76.046
## lf 0.1425 7.829 11.68 12.05 87.951
## northing 13.4019 14.663 22.45 23.22 76.784
## easting 9.2808 15.378 15.48 15.65 84.347
## snowfree_time 7.2309 43.462 76.83 79.56 20.442
## cumgdd0 9.2594 70.304 85.08 86.73 13.271
## cumgdd60d 8.2298 74.220 86.39 87.63 12.368
## albedo.gdd300 90.9769 91.656 93.76 95.86 4.144
## albedo.gdd600 92.6660 93.425 93.63 95.23 4.766
## albedo.gdd900 69.2480 75.862 77.53 78.18 21.821
## frost_severe.gdd300 19.9514 20.108 25.22 98.75 1.249
## frost_severe.gdd600 19.9514 20.108 25.22 98.75 1.249
## frost_severe.gdd900 19.3976 19.463 24.25 97.50 2.497
## frost.gdd300 77.7221 80.320 82.75 83.87 16.135
## frost.gdd600 77.7221 80.320 82.75 83.87 16.135
## frost.gdd900 71.7872 77.360 79.28 80.23 19.772
## gddspeed.gdd300 87.4295 88.890 89.21 90.34 9.663
## gddspeed.gdd600 82.9044 84.455 90.02 90.48 9.516
## gddspeed.gdd900 43.4410 44.255 66.89 66.90 33.104
## radiations.gdd300 67.2292 68.058 75.69 75.71 24.289
## radiations.gdd600 48.4527 51.311 65.40 65.71 34.291
## radiations.gdd900 30.1283 35.658 58.68 60.26 39.743
## rainfall.gdd300 24.6799 57.262 72.23 72.24 27.761
## rainfall.gdd600 10.5000 40.083 73.95 74.72 25.283
## rainfall.gdd900 5.4678 31.146 68.49 68.80 31.196
## snowdays.gdd300 65.4597 74.482 85.70 88.06 11.945
## snowdays.gdd600 65.4597 74.482 85.70 88.06 11.945
## snowdays.gdd900 60.5756 74.283 83.57 85.66 14.335
```
```r
par(mfrow=c(1,2))
s.corcircle(pca$co)
s.corcircle(pca$co, xax=3, yax=4)
```
![](README_files/figure-html/unnamed-chunk-15-1.png)<!-- -->
### Modélisation de la hauteur moyenne
```r
dat = na.omit(releves_all[,c(vars, "hmean")])
dat$lf = as.factor(dat$lf)
dim(dat)
```
```
## [1] 222 30
```
```r
library(randomForest)
```
```
## randomForest 4.6-14
```
```
## Type rfNews() to see new features/changes/bug fixes.
```
```
##
## Attaching package: 'randomForest'
```
```
## The following object is masked from 'package:ggplot2':
##
## margin
```
```
## The following object is masked from 'package:dplyr':
##
## combine
```
```r
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd300 + rainfall.gdd900 + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, ntree = 150, mtry=16)
```
```
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
```
```r
mod
```
```
##
## Call:
## randomForest(formula = hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd300 + rainfall.gdd900 + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, ntree = 150, mtry = 16)
## Type of random forest: regression
## Number of trees: 150
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 13.07629
## % Var explained: 55.78
```
```r
varImpPlot(mod)
```
![](README_files/figure-html/unnamed-chunk-16-1.png)<!-- -->
```r
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd900 + slope + northing + easting + lf, data = dat, ntree = 100, mtry=16)
```
```
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
```
```r
mod
```
```
##
## Call:
## randomForest(formula = hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd900 + slope + northing + easting + lf, data = dat, ntree = 100, mtry = 16)
## Type of random forest: regression
## Number of trees: 100
## No. of variables tried at each split: 8
##
## Mean of squared residuals: 13.10248
## % Var explained: 55.69
```
```r
library(plotmo)
```
```
## Loading required package: Formula
```
```
## Loading required package: plotrix
```
```
## Loading required package: TeachingDemos
```
```r
plotmo(mod, ylim = NA)
```
```
## plotmo grid: gddspeed.gdd300 radiations.gdd300 cumgdd60d rainfall.gdd900
## 39 4890.969 456.8318 2.860266
## slope northing easting lf
## 12.70785 -0.6139406 -0.1065332 5
```
![](README_files/figure-html/unnamed-chunk-17-1.png)<!-- -->
...@@ -177,8 +177,8 @@ head(sites_all) ...@@ -177,8 +177,8 @@ head(sites_all)
reltot = merge(sites_all@data[,c("id_site", "alti", "slope", "lf", "northing", "easting")], releves_all, by.x = "id_site", by.y="ref_site") reltot = merge(sites_all@data[,c("id_site", "alti", "slope", "lf", "northing", "easting")], releves_all, by.x = "id_site", by.y="ref_site")
head(reltot) head(reltot)
saveRDS("reltot", file = "releves_all_plus.rds") saveRDS(reltot, file = "releves_all_plus.rds")
saveRDS("sites_all", file = "sites_all_plus.rds") saveRDS(sites_all, file = "sites_all_plus.rds")
###### analyse sites - type - topo ###### analyse sites - type - topo
library(ggplot2) library(ggplot2)
......
releves_all = readRDS("releves_all.rds") releves_all = readRDS("releves_all_plus.rds")
head(releves_all) head(releves_all)
dim(releves_all)
length(unique(releves_all$releve)) length(unique(releves_all$releve))
length(unique(releves_all$ref_site)) length(unique(releves_all$id_site))
head(releves_all[which(is.na(releves_all$hmean)),])
library(ade4)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900", "frost_severe.gdd300", "frost_severe.gdd600",
"frost_severe.gdd900", "frost.gdd300", "frost.gdd600", "frost.gdd900",
"gddspeed.gdd300", "gddspeed.gdd600", "gddspeed.gdd900", "radiations.gdd300",
"radiations.gdd600", "radiations.gdd900", "rainfall.gdd300", "rainfall.gdd600",
"rainfall.gdd900", "snowdays.gdd300", "snowdays.gdd600", "snowdays.gdd900" )
dat = na.omit(releves_all[,vars])
dim(dat)
pca = dudi.pca(dat, scan=FALSE, nf = 4)
inertia.dudi(pca, row=FALSE, col =TRUE)
par(mfrow=c(1,2))
s.corcircle(pca$co)
s.corcircle(pca$co, xax=3, yax=4)
## ax1 41% albedo (frost?, gddspeed, radiations300, snowdays?)
## ax2 14% cumgdd (snowfree_time, alti, opp. rainfall300-600)
## ax3 10% snowfree_time, rainfall600-900,
## ax4 8% frost severe
# other variables not represented : slope, north/east,
### ==== MODELE
dat = na.omit(releves_all[,c(vars, "hmean")])
dat$lf = as.factor(dat$lf)
dim(dat)
library(randomForest) library(randomForest)
sub = releves_all[releves_all$ref_typoveg=="SUB" & !(is.na(releves_all$hmean)),] mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd300 + rainfall.gdd900 + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, ntree = 150, mtry=16)
que = releves_all[releves_all$ref_typoveg=="QUE" & !(is.na(releves_all$hmean)),] # plot(mod$rsq, type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")
prod = releves_all[releves_all$ref_typoveg=="PROD" & !(is.na(releves_all$hmean)),]
alp = releves_all[releves_all$ref_typoveg=="ALP" & !(is.na(releves_all$hmean)),] # library(caret)
# mod = caret::train(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd0 + rainfall.gdd300 + rainfall.gdd900 + snowfree_time + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, method = "rf", ntree=150)
# print(mod)
# print(mod$finalModel)
mod
varImpPlot(mod)
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd900 + slope + northing + easting + lf, data = dat, ntree = 100, mtry=16)
mod = randomForest(hmean ~ cumgdd.90 + frost.90 + rainfall.90, data = alp)
mod mod
varImpPlot(mod)
plot(hmean~cumgdd.90, data = releves_all) library(plotmo)
plotmo(mod, ylim = NA)
library(ggplot2) # ###
dat = ggplot(dat, aes(x=value, y=hmean)) + # library(randomForest)
geom_point() + # sub = releves_all[releves_all$ref_typoveg=="SUB" & !(is.na(releves_all$hmean)),]
facet_wrap(~period, ncol = 2, scales = "free") # que = releves_all[releves_all$ref_typoveg=="QUE" & !(is.na(releves_all$hmean)),]
# prod = releves_all[releves_all$ref_typoveg=="PROD" & !(is.na(releves_all$hmean)),]
# alp = releves_all[releves_all$ref_typoveg=="ALP" & !(is.na(releves_all$hmean)),]
#
# mod = randomForest(hmean ~ snowfree_time + cumgdd0 + cumgdd60d + , data = alp)
# mod
#
# plot(hmean~cumgdd.90, data = releves_all)
#
# library(ggplot2)
# dat = ggplot(dat, aes(x=value, y=hmean)) +
# geom_point() +
# facet_wrap(~period, ncol = 2, scales = "free")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment