NCL Home > Documentation > Functions > General applied math, Statistics

regline

Calculates the linear regression coefficient between two series.

Prototype

	function regline (
		x [*] : numeric,  
		y [*] : numeric   
	)

	return_val [1] :  float or double

Arguments

x
y

One-dimensional arrays of the same length. Missing values should be indicated by x@_FillValue and y@_FillValue. If x@_FillValue or y@_FillValue are not set, then the NCL default (appropriate to the type of x and y) will be assumed.

Return value

The return value is a scalar of type double if either x or y are double, and float otherwise. Some attributes are returned as well. See the description below.

Description

regline computes the information needed to construct a regression line: regression coefficient (trend, slope,...) and the average of the x and y values. regline is designed to work with one-dimensional x and y arrays. Missing data are allowed.

regline also returns the following attributes:

xave (scalar, float or double, depending on x and y)
average of x
yave (scalar, float or double, depending on x and y)
average of y
tval (scalar, float or double, depending on x and y)
t-statistic (assuming null-hypothesis)
rstd (scalar, float or double, depending on x and y)
standard error of the regression coefficient
yintercept (scalar, float or double, depending on x and y)
y-intercept at x=0
nptxy (scalar, integer)
number of points used

If the regression coefficients for multi-dimensional arrays are needed, use regCoef.

If the series is autocorrelated the returned degrees of freedom must be recalculated. See Example 2.

          Taking Serial Correlation into Account in Tests of the Mean
          F. W. Zwiers and H. von Storch, J. Climate 1995, pp336-

Examples

Example 1

The following example was taken from:

    Brownlee
    Statistical Theory and Methodology
    J Wiley 1965   pgs: 342-346     QA276  .B77

The regression line information for the example below is: (a) rc=0.9746, (b) tval=38.7, (c) nptxy=18 which yields 16 degrees of freedom (df=nptxy-2). To test the null hypothesis (i.e., rc=0) at the two-tailed 95% level, we note that t(16) is 2.120 (table look-up: 0.975). Clearly, the calculated t-statistic greatly exceeds 2.120 so the null hypothesis is rejected at the 5% level.

Rather than a table lookup, the following could be used to calculate the actual significance level.

       alpha = betainc(df/(df+rc@tval^2), df/2.0, 0.5)

or, alternatively,

       prob = 1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5)

The example series are:

      x    = (/ 1190.,1455.,1550.,1730.,1745.,1770. \
             ,  1900.,1920.,1960.,2295.,2335.,2490. \
             ,  2720.,2710.,2530.,2900.,2760.,3010. /)
      
      y    = (/ 1115.,1425.,1515.,1795.,1715.,1710. \
             ,  1830.,1920.,1970.,2300.,2280.,2520. \
             ,  2630.,2740.,2390.,2800.,2630.,2970. /)
      
      rc   = regline (x,y)
     
                                    ; Note use of attributes
      df   = rc@nptxy-2
      prob = (1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5) )
     
      yReg = rc*x + rc@yintercept   ; NCL array notation 
                                    ; yReg is same length as x, y
      print(rc)
     
      print("prob="+prob)
     
      print(yReg)

The "print(rc)" statement will yield:

Variable: rc
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes:   [1]
Coordinates: 
Number Of Attributes: 7
  _FillValue :  -999
  yintercept :  15.35229
  yave       :  2125.278
  xave       :  2165
  nptxy      :  18
  rstd       :  0.025155
  tval       :  38.74286

(0)     0.9745615

(0)     prob=1

The "print(yReg)" statement will yield:

Variable: yReg
Type: float
Total Size: 72 bytes
            18 values
Number of Dimensions: 1
Dimensions and sizes:   [18]
Coordinates: 
Number Of Attributes: 1
  _FillValue :  -999
(0)     1175.08
(1)     1433.339
(2)     1525.923
(3)     1701.344
(4)     1715.962
(5)     1740.326
(6)     1867.019
(7)     1886.51
(8)     1925.493
(9)     2251.971
(10)    2290.953
(11)    2442.01
(12)    2666.159
(13)    2656.414
(14)    2480.993
(15)    2841.581
(16)    2705.142
(17)    2948.782

Note 1: The above assumes that all the points are independent. If this is not the case, then the number used to test for significance should be less.

Note 2: To construct 95% confidence limits for the hypothesis that the regression coefficient is one (i.e., rc=1) :

As noted above, the t for 0.975 and 16 degrees of freedom is 2.120 [table look-up].
rc@rstd * 2.12 = 0.053. This yields 95% confidence limits of (0.97-0.053) < 0.97 < (0.97+0.053) or (0.92 to 1.03). Thus, the hypothesis that rc=1 can not be rejected.

Example 2

Often geophysical time series are correlated. If the auto-correlation is significant [specified, a priori, by the user], then the degrees of freedom returned should be adjusted. The adjusted dof should be used in the calculation of the regression coefficient significance. Example 1 used:

 
      rc   = regline (x,y)
      df   = rc@nptxy-2
      prob = (1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5) )

A-priori, let's set the lag-1 autocorrelation significance level at (say) 0.10.

        N    = rc@nptxy                         ; convenience/clarity
        acr  = esacr(y,2)
        if (acr(1).gt.0.0) then
            pr1     = rtest(acr(1), N, 0)
            rsiglvl = 0.10
            if (pr1.lt.rsiglvl) then
                df   = N*(1.0-acr(1))/(1.0+acr(1))
                prob = (1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5) )
            end if
        end if

Example 3

Use the 'Old Faithful' data set used by several R regression examples Here, we use NCL's regline and the ANOVA version regline_stats.

Simplefied R code:

    eruption.lm = lm(eruptions ~ waiting, data=faithful)
    summary(eruption.lm) 

Summary output:

    Call: 
    lm(formula = eruptions ~ waiting, data = faithful) 
     
    Residuals: 
        Min      1Q  Median      3Q     Max 
    -1.2992 -0.3769  0.0351  0.3491  1.1933 
     
    Coefficients: 
                Estimate Std. Error t value Pr(>|t|) 
    (Intercept) -1.87402    0.16014   -11.7   <2e-16 *** 
    waiting      0.07563    0.00222    34.1   <2e-16 *** 
    --- 
    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 
     
    Residual standard error: 0.497 on 270 degrees of freedom 
    Multiple R-squared: 0.811,      Adjusted R-squared: 0.811

The NCL code is:


   diri = "./"
   fili = "faithful_R.dat"
   pthi = diri + fili

   ncol = 3
   data = readAsciiTable(pthi, ncol, "float", 1)

   x    = data(:,2)   ; eruption interval (minutes); independent
   y    = data(:,1)   ; eruption duration (minutes); dependent
   nxy  = dimsizes(y)

   x@long_name = "Eruption Interval (min)"
   y@long_name = "Eruption Duration (min)"

  ;print("============ OLD Faithful================")
  ;print(sprintf("%4.0f", x)+"  "+sprintf("%4.0f", y))
  ;print("=========================================")

;---Assorted statistics and regression line fit
  ;rxy = escorc(x,y)         ;  cross-correlation
  ;print(rxy)
  ;print("=========================================")

   rcl  = regline(x, y)      ; basic linear regression
   print(rcl)

   df   = rcl@nptxy-2
   prob = (1 - betainc(df/(df+rcl@tval^2), df/2.0, 0.5) )
   print("=========================================")

   rcl_anova =  regline_stats(x,y) ; linear regression: ANOVA
   print(rcl_anova)
   print("=========================================")

The (slightly) edited output is: