[ncl-talk] Spatial plot for Significance of Linear Trend
Adam Phillips
asphilli at ucar.edu
Mon Mar 5 16:23:45 MST 2018
Hi Sri,
I see a couple of things in your script that you should check:
1 - You have the following:
do ne=0,neof-1
eof_regres(ne,:,:)= regCoef(eof_ts(ne,:), t2mw(lat|:,lon|:,time|:));
end do
...
df = onedtond(eof_regres at nptxy,dimsizes(eof_regres))-2 ; degrees of
freedom
tval = onedtond(eof_regres at tval, dimsizes(eof_regres)) ; t-statistic
Note that the @tval and @nptxy attributes get overwritten each time through
the do loop. If neof = 1, the above will work, otherwise the coding should
be changed.
2 - With regards to this error message:
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
I do not know what line 135 is, but the error message above is telling you
how to fix it. The error message is generated because the metadata on the
left hand side of the equal sign is getting overwritten by the metadata on
the right hand side of the equal sign. Note: This error message does not
mean that your calculation of the significance of your trend is incorrect.
However, I do not see where you are calculating a trend, as you are
seemingly regressing a EOF principal component timeseries on a time x lat x
lon array.
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)
tttt = dtrend_msg(ispan(0,dimsizes(tarr&time)-1,1),tarr,True,True)
finobs = tarr(0,:,:)
finobs_rtt = finobs
finobs = (/ onedtond(tttt at slope,
(/dimsizes(obs_seas&lat),dimsizes(obs_seas&lon)/) ) /)
tval = tarr(:,:,0)
df = new((/dimsizes(tarr&lat),dimsizes(tarr&lon)/),"integer")
rc = regcoef(ispan(0,dimsizes(tarr&time)-1,1),tarr,tval,df)
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 %
Hope that helps. If you have any further queries please respond to ncl-talk.
Adam
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
<http://www.cgd.ucar.edu/staff/asphilli>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180305/366c2cbd/attachment.html>
More information about the ncl-talk
mailing list