[R] Peak Over Threshold values

Vito Muggeo (UniPa) vito.muggeo at unipa.it
Wed May 26 11:07:00 CEST 2010


dear Tonja,
By plotting your data

plot(df)

it seems to me that you are looking for a piecewise linear 
relationships. If this is the case, have a look to the package 
segmented. You have to specify or not the number and the starting values 
for the breakpoints

library(segmented)
olm<-lm(walevel~day)

#specify number and starting values for the breakpoints..
oseg<-segmented(olm, seg.Z=~day, psi=c(20,50,80))
plot(oseg,add=TRUE,col=2)
oseg$psi

#or you can avoid to specify starting values for psi
oseg<-segmented(olm, seg.Z=~day, 
psi=NA,control=seg.control(stop.if.error=FALSE,K=30))

plot(oseg,add=TRUE,col=2)
oseg$psi


best,
vito


Tonja Krueger ha scritto:
>    Dear List
> 
>    I hope you can help me: I’ve got a dataframe (df) within which I am looking
>    for Peak Over Threshold values as well as the length of the events. An event
>    starts when walevel equals 5.8 and it should end when walevel equals the
>    lower threshold value (5.35).
> 
>    I tried “clusters (…)” from “evd package”, and varied r (see example) but it
>    did not work for all events (again see example).
> 
>    walevel <- c(5.75, 5.75, 5.75, 5.76, 5.80, 5.82, 5.85, 5.86, 5.91, 5.91,
>    5.88, 5.92, 5.99, 6.02, 6.00, 6.08, 6.11, 6.10, 6.10, 6.11, 6.14, 6.12,
>    6.15, 6.17, 6.15, 6.08, 6.01, 5.95, 5.89, 5.83, 5.77, 5.72, 5.70, 5.65,
>    5.59, 5.51, 5.43, 5.33, 5.30, 5.25, 5.22, 5.21, 5.19, 5.19, 5.19, 5.13,
>    5.15, 5.07, 5.04, 5.02, 4.99, 5.05, 5.07, 5.12, 5.17, 5.22, 5.32, 5.29,
>    5.33, 5.33, 5.36, 5.37, 5.42, 5.43, 5.48, 5.55, 5.57, 5.66, 5.68, 5.72,
>    5.77, 5.79, 5.81, 5.80, 5.80, 5.77, 5.72, 5.70, 5.66, 5.62, 5.62, 5.61,
>    5.59, 5.57, 5.51, 5.47, 5.41, 5.39, 5.40, 5.40, 5.42, 5.45, 5.47, 5.50,
>    5.50, 5.49, 5.43, 5.39, 5.33, 5.26)
> 
>    day <- c(1:100)
> 
>    df <- data.frame(day,walevel)
> 
>    library(evd)
>    clusters(df$walevel, u = 5.80, r = 1, ulow = 5.35, cmax = T, plot = T)
>    clusters(df$walevel, u = 5.80, r = 50, ulow = 5.35, cmax = T, plot = T)
> 
>    What have I done wrong?
> 
>    Tonja
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

-- 
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 23895240
fax: 091 485726/485612
http://dssm.unipa.it/vmuggeo



More information about the R-help mailing list