R: [R] slope estimations of teeth like data

Petr Pikal petr.pikal at precheza.cz
Tue Jun 15 15:21:46 CEST 2004


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 
> > 
> > ______________________________________________ 
> > R-help at stat.math.ethz.ch mailing list 
> > https://www.stat.math.ethz.ch/mailman/listinfo/r-help 
> > PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html 

Petr Pikal
petr.pikal at precheza.cz




More information about the R-help mailing list