R: [R] slope estimations of teeth like data

[Gabor Grothendieck]

I am not entirely sure what is essential and what is inessential in 
this problem.  In the yy data shown, the peaks occur in one step and
all the peaks are equal.  Just plotting

plot(diff(yy))

shows that there is a huge region of potential cutoffs for separating
out the large single step increases from the other differences.

If it continues to be true that the increases are large relative 
to everything else but the cutoff is not quite so obvious as 
here, we can calculate a cutoff by just sorting the unique
differences and partitioning them in every possible
way into small and large differences and then choosing the 
partition whose residual variance is least.   

Its actually pretty simple to do that from first pinciples, as
shown here, but there might be a function in R that does it
even more directly.

diffyy <- diff(yy)
yy. <- unique(sort(diffyy)) # unique diffs sorted
yy. <- yy.[-1] - diff(yy.)/2  # midpoints

# For each 2 class partition into small & large, get residual var
res <- sapply(yy., function(x) 
                     var(resid(lm(diffyy ~ factor(diffyy < x)))) )

cutoff <- yy.[which.min(res)]

# indices of points involved in rise
which(diffyy > cutoff) + 1

---

Petr Pikal <petr.pikal at precheza.cz>

 
On 15 Jun 2004 at 13:52, Vito Muggeo wrote: 

> Dear Petr, 
> Probably I don't understand exactly what you are looking for. 
> 
> However your "plot(x,c(y,z))" suggests a broken-line model for the 
> response "c(y,x)" versus the variables x. Therefore you could estimate 
> a segmented model to obtain (different) slope (and breakpoint) 
> estimates. See the package segmented. 

Thank you Vito, but it is not what I want. plot(x,c(y,z)) shows only one "spike" 
and I have many such spikes in actual data. 

My actual data look like those 

set.seed(1) 
y <- 0.03*x[1:100]+rnorm(100, mean=.001, sd=.03) 
z <- 3-rep(seq(1,100,10),each=10)*.03+rnorm(100,mean=.001, sd=.03) 
yy <- NULL 
for( i in 1:10) yy <- c(yy,c(y,z)[1:floor(runif(1)*200)]) 
y.l <- length(yy) 
plot(1:y.l, yy) 

x axis is actually a time and y is a weight of gradually filled conteiner, which is 
irregularly emptied. I want to do an hourly and/or daily averages of increases in 
weight (it can by done by aggregate) 

myfac <- gl(y.l/12,12,length=1271) #hopefully length is ok 

y.agg <- aggregate(diff(yy), list(myfac), mean) 
## there will be list(hod=cut(time.axis,"hour")) construction actually 

0.03 can be expected average result and some aggregated values ar OK but some 
are wrong as they include values from emptying time. 

*** This*** is probably what I need, I need to set some logical vector which will 
be TRUE when there was a filling time and FALSE during other times. And I 
need to specify it according a data I have available. 

Best what I was able to do was to consider filling time as a time when let say 

diff(yy) >= 0 

was between prespecified limits, but you know how it is with real life and 
prespecified limits. 

Or I can plot my data against time, manually find out regions which are correct 
and make a aggregation only with correct data. But there are 24*60*3 values 
each day so I prefer not to do it manually. 

Or finally I can throw away any hourly average which is not in set limits, but I 
prefer to throw away as little data as possible. 

I hope I was able to clarify the issue a bit. 

Thank you 
Best regards 
Petr 


> 
> best, 
> vito 
> 
> 
> 
> ----- Original Message ----- 
> From: Petr Pikal <petr.pikal at precheza.cz> 
> To: <r-help at stat.math.ethz.ch> 
> Sent: Tuesday, June 15, 2004 1:11 PM 
> Subject: [R] slope estimations of teeth like data 
> 
> 
> > Dear all 
> > 
> > Suppose I have teeth like data similar like 
> > 
> > x <- 1:200 
> > y <- 0.03*x[1:100]+rnorm(100, mean=.001, sd=.03) 
> > z <- 3-rep(seq(1,100,10),each=10)*.03+rnorm(100,mean=.001, sd=.03) 
> > plot(x,c(y,z)) 
> > 
> > and I want to have a gradient estimations for some values from 
> > increasing 
> part of 
> > data 
> > 
> > like 
> > 
> > y.agg <- aggregate(diff(c(y,z)), 
> > list(rep(seq(1,200,10),each=10)[1:199]), 
> mean) 
> > 
> > y.agg[1:10,] ##is OK, I want that 
> > y.agg[11:20,] ##is not OK, I do not want that 
> > 
> > actual data are similar but more irregular and have subsequent 
> > gradual 
> increases 
> > and decreases, more like 
> > 
> > set.seed(1) 
> > yy<-NULL 
> > for( i in 1:10) yy <- c(yy,c(y,z)[1:floor(runif(1)*200)]) 
> > length(yy) 
> > [1] 1098 
> > 
> > plot(1:1098,yy) 
> > 
> > Is there anybody who has some experience with such data, mainly how 
> > to 
> extract 
> > only increasing portions or to filter values of "yy" such as only 
> aggregated slopes 
> > from increasing parts are computed and other parts are set to NA. 
> Sometimes 
> > actual data have so long parts of steady or even slightly increasing 
> values at 
> > decreasing part that aggregated values are slightly positive 
> > although they 
> are 
> > actually from decreasing portion of data. 
> > 
> > Thank you 
> > Petr Pikal 
> > petr.pikal at precheza.cz 
> >