0
$\begingroup$

So I am simply trying to plot observed data and the Theil-Sen trend estimate on one plot, but the estimates have been strange but consistent across methods/packages. As you can see in the figure, I have a time series of averaged yearly minimum temperatures across over 30 stations over 50 years ranging between 5 and 10 degrees C: Station avereaged TMin across valley

When I use different packages in R (trend, mlblm, and this one:https://rdrr.io/github/jentjr/gwstats/man/get_theilsen.html)

I get statistically significant slopes of around 0.0148, but the intercept is always [-28,-21]. I don't understand why I am getting this result, not only from a plotting standpoint (as it would not line up with the data at all) but also given what I know about the median method used to calculate this intercept. Can someone help me understand this result and help me figure out how to plot the trend and observed data together?

Best,

Jacob

edit:

I used Tom's advice and got the following: enter image description here

The intercept seems high, and while I know it is based on the median of pairs rather than the average, that trend line just does not sit right with me.

$\endgroup$
  • $\begingroup$ Some generic discussion in Sections 1 and 2 of stata-journal.com/sjpdf.html?articlenum=st0394 I'd add that re-expressing the slope in terms of deg C/century makes it easier to think about. (To check whether the intercept is too high I think we need access to the data.) $\endgroup$ – Nick Cox Jul 3 '18 at 18:02
  • $\begingroup$ Can you give the data and explain how the procedure you used says it calculates intercept? The line is shifted high as you suggest and doesn't seem to be consistent with any of the Theil-Sen intercept estimators I have seen implemented. By eye it looks like it should be sitting about 0.5 lower (roughly) $\endgroup$ – Glen_b -Reinstate Monica Jul 4 '18 at 0:53
1
$\begingroup$

I am assuming that the solution would be to set 1950 to be 0. As you do not have data between 0 and 1950. As you can see your origin will be way to the left of your x-axis on the plot.

$\endgroup$
  • $\begingroup$ Tom, thank you very much, that makes a whole lot of sense and I fixed the issue and got the plot seen above in the edit. The trend line still looks wrong (intercept too high) but I know it is because that is based off the median of pairs of data. Does that explain why it looks slightly off to me? $\endgroup$ – Jay Coulee Jul 3 '18 at 17:05
  • $\begingroup$ I am not familiar with the method you are using. Maybe it is aiming to stay within some small chosen distance to the largest possible number of points? Also, perhaps the value at 2015 is an outlier? $\endgroup$ – Tom Jul 3 '18 at 17:18
  • $\begingroup$ The theil-sen method basically takes every pair of points and takes the median slope of the lines formed from these combinations and by that same logic, the median intercept from the lines formed. That is why I guess it is strange to me that 2015 would have such a large effect? It's supposed to be robust against this kind of thing but perhaps that is just a coincidence it looks like it is being swayed by 2015 $\endgroup$ – Jay Coulee Jul 3 '18 at 17:27
1
$\begingroup$

The advice by @Tom helps: set 1950 as Year 0, and the results are much more reasonable. (Shown below, blue line). I don't know why this would be, but I also don't know how mblm calculates the intercept. As a note, this problem does not occur with the quantile regression shown below with red line.

Data are approximate.

if(!require(mblm)){install.packages("mblm")}
if(!require(ggplot2)){install.packages("ggplot2")}

Data = read.table(header=T, text="
Year               MinTemp
1950.0382043935053 5.519238641015146
1950.9878564606358 5.918338108882523
1952.116250511666  6.569279574293901
1953.0945558739256 6.527834629553828
1954.2161277118298 7.283667621776505
1955.0211488606906 6.906467458043391
1956.202756174103  6.739255014326648
1956.7649065356802 7.096193205075728
1957.8455450948288 8.48137535816619
1958.6410151453133 8.251023331968891
1959.8376313276028 7.853049529267295
1961.038340837768  7.392140810478921
1961.8024287078729 7.644289807613591
1962.7616318733799 7.8965411379451504
1963.7208350388867 8.14879246827671
1964.748260335653  7.3521285304952935
1965.9394187474418 7.038067949242736
1966.894528585073  7.353254195661073
1968.0652203574841 7.353868194842407
1968.7815527357075 8.34025787965616
1970.0109155410016 7.4388047482603366
1970.9714831491337 7.670077773229637
1971.9102196752629 8.23700368399509
1973.0863692181745 8.153704461727385
1974.0769545640608 7.923454768726976
1974.8342202210397 8.280495292672944
1976.021285304953  8.029369627507165
1976.9791240278348 8.302599263200984
1977.9437849638423 8.470937372083505
1978.9384636376042 8.177752762996317
1979.701187065084  8.450880065493248
1981.2593805430483 8.493655341792879
1982.0412061672807 8.473086369218175
1983.0126893164143 8.536532951289399
1983.8095238095239 8.28520261972984
1984.8042024832857 7.992018010642654
1985.7715923045437 8.118399508800655
1986.726702142175  8.433585755218994
1988.1129758493655 8.119627507163324
1988.9411925228544 7.385796152271798
1989.9249556556147 7.260437986082687
1991.032883067267  8.226054031927958
1991.9948151180245 8.436348751534998
1992.9826715786603 8.24805566925911
1993.8272615636513 7.262484650020468
1994.785100286533  7.5357142857142865
1995.9953608950746 6.927957429390096
1996.7075999454223 7.9772820302906275
1997.6627097830537 8.292468276708965
1998.8402237685907 8.188190749079002
1999.817164688225  8.167724109701188
2000.8336744439898 7.538886614817848
2002.0193750852777 7.308739255014327
2002.756174102879  7.980454359394189
2003.7440305635148 7.792161277118298
2004.8833401555466 8.275276299631601
2005.6815390912814 8.002967662709784
2006.6707599945423 7.7936962750716345
2008.0556692591076 7.500716332378224
2009.4815118024287 9.578387228817029
2010.7313412471005 8.362259516987312
2011.7055532814845 8.383749488334017
2012.660663119116  8.698935734752354
")

Data$Year = round(Data$Year)

Data$Year = Data$Year - 1950


library(mblm)

model= mblm(MinTemp ~ Year, data=Data)

summary(model)

   ### Coefficients:
   ###             Estimate      MAD V value Pr(>|V|)    
   ### (Intercept) 7.820528 0.499264    2016 5.29e-12 ***
   ### Year        0.009301 0.012388    1758 2.88e-07 ***

Sum = summary(model)$coefficients

library(ggplot2)

ggplot(Data, aes(x=Year, y=MinTemp)) + 
  geom_point() +
  geom_abline(intercept = Sum[1], slope = Sum[2], color="blue", size=1.2) +
  labs(x = "Years after 1950")

enter image description here

### Optional quantile regression follows

if(!require(quantreg)){install.packages("quantreg")}

library(quantreg)

model.q = rq(MinTemp ~ Year, data = Data, tau = 0.5)

summary(model.q)

   ### Coefficients:
   ###             coefficients lower bd upper bd
   ### (Intercept) 7.46789      6.97021  7.78631 
   ### Year        0.01470      0.00656  0.02508 

library(ggplot2)

model.null = rq(MinTemp ~ 1, data = Data, tau = 0.5)

anova(model.q, model.null)

   ### Quantile Regression Analysis of Deviance Table
   ### 
   ### Model 1: MinTemp ~ Year
   ### Model 2: MinTemp ~ 1
   ###   Df Resid Df F value  Pr(>F)  
   ### 1  1       61  5.7424 0.01964 *

Sumq = summary(model.q)$coefficients

ggplot(Data, aes(x=Year, y=MinTemp)) + 
  geom_point() +
  geom_abline(intercept = Sumq[1], slope = Sumq[2], color="red", size=1.2) +
  labs(x = "Years after 1950")

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.