[ncl-talk] Spatial plot for Significance of Linear Trend

[Adam Phillips]

Hi Sri,
First, I highly recommend you speak with a statistician or your advisor
about the issues we have been discussing. Also review the documentation of
the NCL functions regCoef and dtrend. It is important that you fully
understand the particulars about these mathematical functions and to know
exactly what your code is doing.

To answer the questions from your first reply:

> Just to clarify here the *regression coefficient* is synonymous with the
slope or *trend*?
> The main aim is to plot out the trend and overlay with p values.

Did you try testing it? Here's a simple test:
a = (/4.5,3.2,-1.5,0.5,1.7,2.9/)
rc = regCoef(ispan(0,dimsizes(a)-1,1),a)
tt = dtrend(a,True)
print(rc+" "+tt at slope)

Again, please review the documentation for these functions, and try them
out using simple tests as is shown above.

> Is there another way to calculate the significance apart from:
> finobs_rtt = (/ (1.-betainc(df/(df+tval^2), df/2.0, beta_b))*100. /)  ;
probability in %?

I do not know of one, but you should consult a statistician for this.

> Is it possible to do ttest on the trend analysis?

You can do a t-test if you instead choose to do an epoch difference
calculation instead of a trend calculation. In an epoch difference
calculation you (for example) would subtract the average of the first 15
years of your data from the average of the last 15 years of data. Ask a
colleague or your advisor on whether this would be appropriate for your
situation.
---------------------
With regards to the code you attached in your last email, on first glance I
do not see anything amiss with your script. However, as I noted in a
previous reply, the slope attribute returned from dtrend is in units of per
time step. For your script, as you have 80 years of annual 2m Temperature
data,  you will want to multiply the returned slope from dtrend by 80.
finobs1 = finobs1*80.
Adam

On Wed, Mar 7, 2018 at 4:00 AM, Sri Nandini <snandini at marum.de> wrote:

>
> Hello
>  Using the given code, i attach the 2 pltos i get from doing trend
> analysis on T2m and getting the significance. IT seems really small and i
> am not convinced these are the right values.
> Could someone please have a look at the script below which i used and
> confirm i use the right scripting for getting getting significance of and
> linear trends for plotting.
> Much appreciated
>
> Here is the part of my code that works:
>
> yrStrt = 2020
> yrLast = 2099
> ;===================================================================
> f= addfile("T2M_45.nc", "r") ;
> time   = f->time
> YYYY   = cd_calendar(time,-1)/100
> ; entire file
> iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
> T2    = f->TREFHT(iYYYY,:,:)
> printVarSummary(T2)
> ; (time, lat,lon)
> T2=T2-273.15
> T2 at _FillValue = -9.96921e+36
> T2 at units= "degC"
> printVarSummary(T2)
> ; (time, lat,lon)
>
> ;=======================================================================
> ;Very important to give the year dimension its cordinates
> ;=======================================================================
>   yyyymm = cd_calendar(T2&time, -1)
>    yyyy   = yyyymm/100
>
>    dimx = dimsizes(T2)
>    ntim = dimx(0)                ; all years and months
>    nlat = dimx(1)
>    mlon = dimx(2)
>
>    year  = ispan(yyyy(0), yyyy(ntim-1), 1)
>    nyrs  = dimsizes(year)
>
> ;=========================================================================
>  T2mm    = month_to_annual(T2, 1)
>  T2mm&year  = year
>  printVarSummary(T2mm)
> ;=============================================================
> ; Calculate the regression coefficients (slopes)
> ;=============================================================
>
> tarr1 = T2mm(lat|:,lon|:,year|:)
>      printVarSummary(tarr1)
>                                                          ; (time,
> lat,lon)
> tttt1 = dtrend_msg(ispan(0,dimsizes(tarr1&year)-1,1),tarr1,True,True)
>         printVarSummary(tttt1)
> finobs1 = tarr1(:,:,0)
> finobs_rtt1 = finobs1
> finobs1 = (/ onedtond(tttt1 at slope, (/dimsizes(tarr1&lat),dimsizes(tarr1&lon)/)
> ) /)   ; reform the slope attribute to 2D (lat x lon) array. This will be
> used for plotting the trend
>        printVarSummary(finobs1)
>
> tval1 = tarr1(:,:,0)
> df1 = tarr1(:,:,0)
>
> rc1 = regcoef(ispan(0,dimsizes(tarr1&year)-1,1),tarr1,tval1,df1)
> ; tval and df are filled in by regcoef
>
> df1 = equiv_sample_size(tarr1,0.05,0)
> ; Optional! over rule df returned from regcoef.
>
> df1 = df1-2
>                                       ;regcoef/equiv_sample_size return
> N, need N-2
> beta_b1 = tval1
> beta_b1 = 0.5
>
> printVarSummary(df1)
> ;lat,lon
> printVarSummary(tval1)
> ;lat,lon
>
>  finobs_rtt1 = (/ (1.-betainc(df1/(df1+tval1^2), df1/2.0, beta_b1))*100.
> /) ; significance of trends expressed from 0 to 100%
> printVarSummary(finobs_rtt1)
> ;print(finobs_rtt1)                                        ; (time,
> lat,lon)
>
> _______________________________________________
> 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/20180307/366de48e/attachment.html>