[ncl-talk] Spatial plot for Significance of Linear Trend
Adam Phillips
asphilli at ucar.edu
Tue Mar 6 14:30:30 MST 2018
Hi Sri,
I've again pasted my code below but added numerous comments to answer your
questions. I had a mistake in my original code, and changed this:
finobs = tarr(0,:,:)
to this:
finobs = tarr(:,:,0)
;======================================================
tarr = obs_ann(:,:,5::12) ; tarr will be dimensioned lat x lon x time
tttt = dtrend_msg(ispan(0,dimsizes(tarr&time)-1,1),tarr,True,True) ; I
am using dtrend_msg to calculate the linear trends of the tarr array. This
has nothing to do with the significance calculation.
finobs = tarr(:,:,0) ; preallocate space for finobs (lat x lon)
finobs_rtt = finobs ; preallocate space for finobs_rtt (lat x lon)
finobs = (/ onedtond(tttt at slope,
(/dimsizes(obs_seas&lat),dimsizes(obs_seas&lon)/)
) /) ; reform the slope attribute to 2D (lat x lon) array.
tval = tarr(:,:,0) ; preallocate space for tval (lat x lon)
df = new((/dimsizes(tarr&lat),dimsizes(tarr&lon)/),"integer") ;
preallocate space for df
rc = regcoef(ispan(0,dimsizes(tarr&time)-1,1),tarr,tval,df) ; tval and
df are filled in by regcoef
df = equiv_sample_size(tarr,0.05,0) ; Optional! over rule df returned
from regcoef. Up to you whether you want to use equiv_sample_size or not
df = df-2 ;regcoef/equiv_sample_size return N, need N-2
beta_b = tval
beta_b = 0.5
finobs_rtt = (/ (1.-betainc(df/(df+tval^2), df/2.0, beta_b))*100. /) ;
probability in %
In the above finobs will have units of (for SST) degrees C per year.
Note that there are multiple ways to compute the linear trend. I choose to
use dtrend (or dtrend_msg or dtrend_msg_n). You can regress your data
against a straight line (as your regress_4 example does), or you can use
the slopes returned from dtrend/dtrend_msg as is illustrated in example #2
here:
https://www.ncl.ucar.edu/Document/Functions/Built-in/dtrend.shtml
Note that on the same page example #4 shows one using dtrend (for the
trends) and then regcoef/betainc to calculate the significance of the
trends.
Hope that makes sense. If not, please respond to the ncl-talk email list.
Adam
On Tue, Mar 6, 2018 at 1:18 AM, Sri Nandini <snandini at marum.de> wrote:
> Hello
>
> Indeed i was regressing a EOF principal component timeseries on a
> detrended time x lat x lon array (temperature).
>
> However, i managed to make that work by setting the right metadata.
>
>
> Now is for the calculation of significance of trends.
>
> I have some questions regarding the code below, posted next to the line in
> the code:
>
>
> Here is some code that I use to calculate the significance of trends
> (using the betainc function as you were doing):
>
> tarr = obs_ann(:,:,5::12)
> =====>what is the data format? [lat, lon, time?]
> tttt = dtrend_msg(ispan(0,dimsizes(tarr&time)-1,1),tarr,True,True)
> ====>Are you using the dtrend_msg function because it returns attributes
> like @slope which is what i need for graphical purposes? why not just use
> dtrend function?
>
>
> finobs = tarr(0,:,:)
> ====>have you rearranged the inital data here?
> finobs_rtt = finobs
> finobs = (/ onedtond(tttt at slope, (/dimsizes(obs_seas&lat),dimsizes(obs_seas&lon)/)
> ) /) ===>i understand we get the slope from the detrend function and not
> from regcoef?
>
> tval = tarr(:,:,0) ====>is this still arranged as lat, lon,
> time?
> df = new((/dimsizes(tarr&lat),dimsizes(tarr&lon)/),"integer")
> =====>defining new variable===degrees of freedom
>
>
>
>
> rc = regcoef(ispan(0,dimsizes(tarr&time)-1,1),tarr,tval,df) ======> i
> dont understand what is happening here? from NCL online, if you want to
> calculate trends https://www.ncl.ucar.edu/Applications/Scripts/regress_
> 4.ncl shows to use
>
> ;************************************************
> ; Calculate the regression coefficients (slopes)
> ;************************************************
> rc = regCoef(year,xann(lat|:,lon|:,year|:))
>
> rc at long_name = "Trend"
> rc at units = xann at units+"/year"
> copy_VarCoords(xann(0,:,:), rc) ; copy lat,lon coords
>
> printVarSummary(rc)
>
>
> df = equiv_sample_size(tarr,0.05,0) ; Optional! over rule df returned
> from regcoef. Up to you whether you want to use equiv_sample_size or not
> df = df-2 ;regcoef/equiv_sample_size return N, need N-2
> beta_b = tval
> beta_b = 0.5
> finobs_rtt = (/ (1.-betainc(df/(df+tval^2), df/2.0, beta_b))*100. /)
> ====> is this ttest? with a 2 tailed probability?
>
> Would appreciate is someone would help guide me on this.
>
> The idea is to just to do a simple trend analysis on temperature from
> 2020-2100 (80years) and see if its significant or not?
>
>
>
>
>
>
>
> On Mon, Mar 5, 2018 at 2:36 AM, Sri Nandini <snandini at marum.de> wrote:
>
> Hello
>
> I am trying to *calculate statistical significance* of a *linear* trend
> (Regression analysis).
> The regression analysis simple enough.
>
> I have troubles understanding how the probabilities are computed based on
> the trend i calculate?
> What i actually wanted to do is calculate the 2 tailed ttest probabilities
> but it doesn't seem possible with NCL, unless someone else knows of a way?
> Or does
>
> *using the NCL function "betainc" to get the p-value for a student test is
> the same thing?* *>>>*My code is below:
>
>
> f = addfile ("air.mon.mean.nc", "r")
> latS = 29.3
> latN = 67.26
> lonL = 19.
> lonR = 60.
>
> TIME1 = f->time
> YYYY1 = cd_calendar(TIME1,-1)/100 ; entire file
> iYYYY1 = ind(YYYY1.ge.yrStrt .and. YYYY1.le.yrLast)
> t2m = f->air(iYYYY1,:,:)
> printVarSummary(t2m)
>
> T2M = month_to_season (t2m, season)
> printVarSummary(T2M)
> nyrs = dimsizes(T2M&time)
>
> t2mw = T2M({lat|latS:latN},{lon|lonL:lonR},time|:)
> printVarSummary(t2mw)
>
> t2mw = dtrend(t2mw,False)
> printVarSummary(t2mw)
>
> eofT = eofunc_Wrap(t2mw, neof, optEOF)
> printVarSummary(eofT)
>
> ;===========================================================
> ===============
> eof_regres= eofT ;create an array with meta
> data
> do ne=0,neof-1
> eof_regres(ne,:,:)= regCoef(eof_ts(ne,:), t2mw(lat|:,lon|:,time|:));
> end do
> printVarSummary(eof_regres)
>
> ;===========================================================
> ===============
> ;Calculate significance at 95% confidence
> ;The regCoef function returns the following attributes: yintercept
> (y-intercept), tval (t-statistic), rstd (standard error of the estimated
> regression coefficient) and nptxy (number of elements used).
> ;===========================================================
> ===============
>
> df = onedtond(eof_regres at nptxy,dimsizes(eof_regres))-2 ;
> degrees of freedom
> tval = onedtond(eof_regres at tval, dimsizes(eof_regres)) ;
> t-statistic
>
> b = tval ; b must be same size as
> tval (and df)
> b = 0.5
> prob = betainc(df/(df+tval^2),df/2.0,b)
> printVarSummary(prob) ; prob(nlat,nlon)
>
> ;copy_VarCoords(eof_regres,prob)
> prob at long_name = "probability"
> ;prob at units = "fraction: [0,1]"
>
> prob_95=where (prob.ge.0.95,1,0) ; **** note 0.95 ****
> copy_VarCoords(prob, prob_95)
> printVarSummary(prob_95)
>
> eof_regres at long_name = "regression coefficient"
> prob at long_name = "probability"
>
> ;==================================The error message i get is:
>
> warning:VarVarWrite: rhs has no dimension name or coordinate variable,
> deleting name of lhs dimension number (1) and destroying coordinate var,
> use "(/../)" if this is not desired outcome
> warning:VarVarWrite: rhs has no dimension name or coordinate variable,
> deleting name of lhs dimension number (2) and destroying coordinate var,
> use "(/../)" if this is not desired outcome
> warning:["Execute.c":8640]:Execute: Error occurred at or near line 135 in
> file NAOregreT2M_obs.ncl
>
>
> Variable: eof_regres
> Type: float
> Total Size: 49856 bytes
> 12464 values
> Number of Dimensions: 3
> Dimensions and sizes: [evn | 2] x [76] x [82]
> Coordinates:
> evn: [1..2]
> Number Of Attributes: 11
> tval : <ARRAY of 6232 elements>
> yintercept : <ARRAY of 6232 elements>
> rstd : <ARRAY of 6232 elements>
> nptxy : <ARRAY of 6232 elements>
> long_name : EOF: DJF: Monthly mean of surface temperature
> _FillValue : -9.96921e+36
> method : transpose
> matrix : covariance
> pcvar : ( 59.16475, 20.32462 )
> eval : ( 18982.65, 6521.03 )
> eval_transpose : ( 169.488, 58.22348 )
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
>
> --
> Adam Phillips
> Associate Scientist, Climate and Global Dynamics Laboratory, NCAR
> www.cgd.ucar.edu/staff/asphilli/ 303-497-1726 <(303)%20497-1726>
>
> <http://www.cgd.ucar.edu/staff/asphilli>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
--
Adam Phillips
Associate Scientist, Climate and Global Dynamics Laboratory, NCAR
www.cgd.ucar.edu/staff/asphilli/ 303-497-1726
<http://www.cgd.ucar.edu/staff/asphilli>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180306/dd25aff1/attachment.html>
More information about the ncl-talk
mailing list