[Biomod-commits] Need help with some basic questions
Wilfried Thuiller
wilfried.thuiller at ujf-grenoble.fr
Thu Feb 23 08:21:03 CET 2012
>
> So I have averaged the binary prop.mean.weighted projections (shown below) and next I want to average the thresholds (from consensus_cccm80a2a_results) correct? How do I navigate to these thresholds. And once I have taken the average can I visualize the averaged output using level.plot?
>
> MyAverage=apply(cbind(Total_consensus_cccm80a2a_Bin[,1,"prob.mean.weighted"], Total_consensus_had80a2a_Bin[,1,"prob.mean.weighted"]),1,mean)
The total consensus projections are simply the average of the weighted probabilities means of all repetitions and pseudo-absence projections.
Just do the same for the threshold:
MyTre = mean(consensus_cccm80a2a_results $Cflorida$Cflorida$thresholds$prob.mean.weighted)
then, simply run the BinaryTransformation function
MyAverageBin = BinaryTransformation(MyAverage, MyTre)
Something you might want to do as well is to remove the projections from the models made on 100% of the data. Here when you run the Ensemble Forecasting, you have the projections from the model trained on 100% of the data and and the ones trained on 70% of the data. Better to avoid mixing the two on the final outcome you want to look at.
Just add one option into the Ensemble Forecasting option:
Ensemble.Forecasting(Proj.name="cccm80a2a", weight.method='Roc', PCA.media=F, binary=T, bin.method='Roc', Test=T, decay="proportional", repetition.models=T, final.model.out=TRUE)
### note there was a mistake in the help file, it was said that final.model.out=TRUE keeps the final model in. This is the way around. I have just updated the help file.
Additionnaly, I will also recommend to use the decay = "proportional". The weights are awarded for each method proportionaly to their evaluation scores. The advantage is that the discrimination is more fair than with the decay. In the latter case, close scores can strongly diverge in the weights they are awarded, when the proportional method will consider them as being fairly similar in prediction quality and award them a similar weight.
For the level.plot, once you have made your averaging, you should have a vector of probability or binary values
You should also have your coordinates from somewhere
level.plot(MyAverage, MyCoordinates)
Cheers,
Wilfried
>
> Thanks so much for your help!
>
>
>
>
>
>
>
> On Wed, Feb 22, 2012 at 12:36 PM, Wilfried Thuiller <wilfried.thuiller at ujf-grenoble.fr> wrote:
> Dear Ashley,
>
> The weights, thresholds etc are the same because they are estimated from the models, not from the projections. The Ensemble.Forecasting function does not re-run the models, it just pick up the results from every model run in Models, and then estimate the weights and threshold. These weights are then used to make the ensemble forecast. Obviously, the weighting scheme is not going to change in function of the climate change scenario you use because it only relates to the models themselves that are independent of the climate change scenarios.
>
> Concerning your second question, of course it could be done easily. BIOMOD does not do the for you but this is rather easy to compute.
>
> Averaging needs to be done with the probability values (not the binary ones).
>
> For instance, if you want to average your two projections for Cflorida:
>
> load the Total_consensus projections for each GCMs. They should be named along these lines
>
> load("MYPATH/Total_consensus_cccm80a2a")
>
> Let's say your second GCM is called hadcm3
>
> load("MYPATH/Total_consensus_hadcm3a2a")
>
> then let's say you want to average the prob.mean.weighted projections.
>
> Each Total_consensus file is a 3-D array. In the second dimension, you should just have one column (one species). The third dimension corresponds to the different ensemble forecast strategies.
>
> We are going to use the mean function between the two datasets, using the apply function.
>
> MyAverage = apply(cbind(Total_consensus_cccm80a2a[,1,"prob.mean.weighted"], Total_consensus_hadcm3a2a[,1,"prob.mean.weighted"]), 1, mean)
>
> For the threshold, just do the same and take the average.
>
> Does it help?
>
> Wilfried
>
>
>
>
>
>
>
> Le 22 févr. 2012 à 18:18, Ashley Brooks a écrit :
>
>> Hi All!
>>
>> I have a few questions that I'm a little confused about and need some clarification:
>>
>> 1) I have run a few projections and then ensemble forecasts using different sets of future climate data. However the weights, thresholds, etc for the different runs have the same exact output. Is this normal or have I messed something up somewhere? See below:
>>
>> > Initial.State(Response=my.data[1], Explanatory=my.data[,4:13], IndependentResponse=my.data[1], IndependentExplanatory=my.data[,4:13])
>>
>> > Models(GLM=T, TypeGLM="poly", Test="AIC", GBM=T, No.trees=2000, GAM=T, Spline=3, CTA=T, CV.tree=50, ANN=T, CV.ann=2, SRE=T,quant=0.05, FDA=T, MARS=T, RF=T, NbRunEval=3, DataSplit=70, Yweights=NULL, Roc=T, Optimized.Threshold.Roc=T, Kappa=T, TSS=T, KeepPredIndependent=T, VarImport=5, NbRepPA=2, strategy="circles", coor=Coor, distance=2, nb.absences=1000)
>>
>> > Projection(Proj=cccm80a2a[,3:12], Proj.name='cccm80a2a', GLM=T, GBM=T, GAM=T, CTA=T, ANN=T, SRE=T, quant=0.05, FDA=F, MARS=T, RF=T, BinRoc=T, BinKappa=F, BinTSS=F, FiltRoc=T, FiltKappa=F, FiltTSS=F, repetition.models=T)
>>
>> > Ensemble.Forecasting(Proj.name="cccm80a2a", weight.method='Roc', PCA.media=F, binary=T, bin.method='Roc', Test=T, decay=1.6, repetition.models=T)
>> Cflorida
>>
>> consensus_cccm80a2a_results
>> $Cflorida
>> $Cflorida$weights
>> ANN CTA GAM GBM GLM MARS FDA RF SRE
>> PA1 0.0371 0.0232 0.1522 0.2434 0.0951 0.0594 0 0.3895 0
>> PA1_rep1 0.0371 0.0232 0.1522 0.0594 0.0951 0.3895 0 0.2434 0
>> PA1_rep2 0.0594 0.0232 0.1522 0.3895 0.0951 0.0371 0 0.2434 0
>> PA1_rep3 0.0232 0.0483 0.0951 0.2434 0.0483 0.1522 0 0.3895 0
>> PA2 0.0371 0.0232 0.0951 0.1522 0.0594 0.3165 0 0.3165 0
>> PA2_rep1 0.1022 0.0232 0.1022 0.2434 0.0371 0.3895 0 0.1022 0
>> PA2_rep2 0.0371 0.0232 0.0951 0.1522 0.0594 0.2434 0 0.3895 0
>> PA2_rep3 0.0302 0.0302 0.0951 0.1522 0.0594 0.3165 0 0.3165 0
>>
>> $Cflorida$PCA.median
>> model.selected
>> PA1 NA
>> PA1_rep1 NA
>> PA1_rep2 NA
>> PA1_rep3 NA
>> PA2 NA
>> PA2_rep1 NA
>> PA2_rep2 NA
>> PA2_rep3 NA
>>
>> $Cflorida$thresholds
>> PA1 PA1_rep1 PA1_rep2 PA1_rep3 PA2 PA2_rep1 PA2_rep2 PA2_rep3
>> prob.mean 552.4331 476.5176 476.5653 459.1243 485.8907 485.8569 428.8626 292.2260
>> prob.mean.weighted 567.7707 305.4160 432.0515 423.9614 447.9690 351.9328 297.7658 237.1679
>> median 599.1970 529.4070 492.7660 582.6600 594.6120 520.6430 497.5000 301.3960
>> Roc.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>> Kappa.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>> TSS.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>>
>> $Cflorida$test.results
>> PA1 PA1_rep1 PA1_rep2 PA1_rep3 PA2 PA2_rep1 PA2_rep2
>> prob.mean 0.9986525 0.9958785 0.9951836 0.9978842 0.9985311 0.9991299 0.9975960
>> prob.mean.weighted 0.9995932 0.9970734 0.9957514 0.9991808 0.9993503 0.9987062 0.9987006
>> median 0.9982740 0.9945085 0.9942486 0.9974350 0.9977938 0.9975819 0.9970339
>> Roc.mean 0.9979492 0.9911780 0.9908023 0.9960763 0.9982514 0.9969859 0.9920226
>> Kappa.mean 0.9991215 0.9951836 0.9917288 0.9977910 0.9988672 0.9991667 0.9926667
>> TSS.mean 0.9987797 0.9946638 0.9928955 0.9974322 0.9986610 0.9993701 0.9954548
>> PA2_rep3
>> prob.mean 0.9971582
>> prob.mean.weighted 0.9979944
>> median 0.9972599
>> Roc.mean 0.9948672
>> Kappa.mean 0.9934689
>> TSS.mean 0.9931836
>>
>> And then I get the same exact results for other runs as well:
>>
>> > Projection(Proj=cccm80b2b[,3:12], Proj.name='cccm80b2b', GLM=T, GBM=T, GAM=T, CTA=T, ANN=T, SRE=T, quant=0.05, FDA=F, MARS=T, RF=T, BinRoc=T, BinKappa=F, BinTSS=F, FiltRoc=T, FiltKappa=F, FiltTSS=F, repetition.models=T)
>>
>> > Ensemble.Forecasting(Proj.name="cccm80b2b", weight.method='Roc', PCA.media=F, binary=T, bin.method='Roc', Test=T, decay=1.6, repetition.models=T)
>> Cflorida
>>
>> consensus_cccm80b2b_results
>> $Cflorida
>> $Cflorida$weights
>> ANN CTA GAM GBM GLM MARS FDA RF SRE
>> PA1 0.0371 0.0232 0.1522 0.2434 0.0951 0.0594 0 0.3895 0
>> PA1_rep1 0.0371 0.0232 0.1522 0.0594 0.0951 0.3895 0 0.2434 0
>> PA1_rep2 0.0594 0.0232 0.1522 0.3895 0.0951 0.0371 0 0.2434 0
>> PA1_rep3 0.0232 0.0483 0.0951 0.2434 0.0483 0.1522 0 0.3895 0
>> PA2 0.0371 0.0232 0.0951 0.1522 0.0594 0.3165 0 0.3165 0
>> PA2_rep1 0.1022 0.0232 0.1022 0.2434 0.0371 0.3895 0 0.1022 0
>> PA2_rep2 0.0371 0.0232 0.0951 0.1522 0.0594 0.2434 0 0.3895 0
>> PA2_rep3 0.0302 0.0302 0.0951 0.1522 0.0594 0.3165 0 0.3165 0
>>
>> $Cflorida$PCA.median
>> model.selected
>> PA1 NA
>> PA1_rep1 NA
>> PA1_rep2 NA
>> PA1_rep3 NA
>> PA2 NA
>> PA2_rep1 NA
>> PA2_rep2 NA
>> PA2_rep3 NA
>>
>> $Cflorida$thresholds
>> PA1 PA1_rep1 PA1_rep2 PA1_rep3 PA2 PA2_rep1 PA2_rep2 PA2_rep3
>> prob.mean 552.4331 476.5176 476.5653 459.1243 485.8907 485.8569 428.8626 292.2260
>> prob.mean.weighted 567.7707 305.4160 432.0515 423.9614 447.9690 351.9328 297.7658 237.1679
>> median 599.1970 529.4070 492.7660 582.6600 594.6120 520.6430 497.5000 301.3960
>> Roc.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>> Kappa.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>> TSS.mean 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
>>
>> $Cflorida$test.results
>> PA1 PA1_rep1 PA1_rep2 PA1_rep3 PA2 PA2_rep1 PA2_rep2 PA2_rep3
>> prob.mean 0.9986525 0.9958785 0.9951836 0.9978842 0.9985311 0.9991299 0.9975960 0.9971582
>> prob.mean.weighted 0.9995932 0.9970734 0.9957514 0.9991808 0.9993503 0.9987062 0.9987006 0.9979944
>> median 0.9982740 0.9945085 0.9942486 0.9974350 0.9977938 0.9975819 0.9970339 0.9972599
>> Roc.mean 0.9979492 0.9911780 0.9908023 0.9960763 0.9982514 0.9969859 0.9920226 0.9948672
>> Kappa.mean 0.9991215 0.9951836 0.9917288 0.9977910 0.9988672 0.9991667 0.9926667 0.9934689
>> TSS.mean 0.9987797 0.9946638 0.9928955 0.9974322 0.9986610 0.9993701 0.9954548 0.9931836
>>
>>
>>
>> 2) Also I am using two GCM (cccma and had) with two strylines each (A2A and B2B) and would like to average the two GCM's to have just two predictions (one for the A2A storyline and one for the B2B storyline). My question is how do I average the two GCM's? Can I do this in Biomod or do I do it in another program? What commands would I use?
>>
>> Any advice or help would be most appreciated!
>> Thanks,
>> Ashley
>>
>> _______________________________________________
>> Biomod-commits mailing list
>> Biomod-commits at lists.r-forge.r-project.org
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/biomod-commits
>
> --------------------------
> Dr. Wilfried Thuiller
> Laboratoire d'Ecologie Alpine, UMR CNRS 5553
> Université Joseph Fourier
> BP53, 38041 Grenoble cedex 9, France
> tel: +33 (0)4 76 51 44 97
> fax: +33 (0)4 76 51 42 79
>
> Email: wilfried.thuiller at ujf-grenoble.fr
> Personal website: http://www.will.chez-alice.fr
> Team website: http://www-leca.ujf-grenoble.fr/equipes/emabio.htm
>
> ERC Starting Grant TEEMBIO project: http://www.will.chez-alice.fr/Research.html
> FP6 European EcoChange project: http://www.ecochange-project.eu
>
>
>
>
>
>
>
>
> _______________________________________________
> Biomod-commits mailing list
> Biomod-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/biomod-commits
--------------------------
Dr. Wilfried Thuiller
Laboratoire d'Ecologie Alpine, UMR CNRS 5553
Université Joseph Fourier
BP53, 38041 Grenoble cedex 9, France
tel: +33 (0)4 76 51 44 97
fax: +33 (0)4 76 51 42 79
Email: wilfried.thuiller at ujf-grenoble.fr
Personal website: http://www.will.chez-alice.fr
Team website: http://www-leca.ujf-grenoble.fr/equipes/emabio.htm
ERC Starting Grant TEEMBIO project: http://www.will.chez-alice.fr/Research.html
FP6 European EcoChange project: http://www.ecochange-project.eu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/biomod-commits/attachments/20120223/5936dea6/attachment-0001.html>
More information about the Biomod-commits
mailing list