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

Sri Nandini snandini at marum.de
Tue Mar 6 23:25:05 MST 2018 

Thank you

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.

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 %?

Is it possible to do ttest on the trend analysis?



On Mar 6, 2018 10:30:30 PM, Adam Phillips wrote:

> 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. U> p 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. U> > > p 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> > >  
> > > 

> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > > > > 
> > > 
> > 

> > _______________________________________________
> > 
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 
> 

> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> > 
> 



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180307/fffebcb4/attachment.html>

More information about the ncl-talk mailing list