On Fri, May 14, 2010 at 1:47 PM, Thomas Levine <thomas.levine@gmail.com>
wrote:
> Creating the 5 indicator variables will be easy if you post your code
> and sample data. This may also allow people to help with the first
> problem you were having.

Here you go.

Fragment of data file with 4000 or so pts:

RawTime Time OxpcA THCA OxyA DeOxyA Aux1 Aux2 Aux3
9213290       0 71.2   51.47   36.64   14.82   3.88  3.88  0.35
9213390     0.1 66.1   46.26   30.59   15.67   3.88  3.88  0.35
9213490     0.2 81.8   54.07   44.25   9.815   3.88  3.88  0.35
9213590     0.3 58.1    56.7   32.96   23.74   3.88  3.88  0.35
9213690     0.4 50.2   48.55   24.39   24.16   3.88  3.88  0.35
9213790     0.5 59.7   47.02   28.08   18.94   3.88  3.88  0.35
9213890     0.6 68.0   52.78    35.9   16.88  3.88  3.88  0.35
9213990     0.7 60.7   46.62   28.28   18.34   3.88  3.88  0.35
9214090     0.8 57.5   47.81   27.51   20.31   3.88  3.88  0.35
9214190     0.9 60.4   52.13   31.49   20.64   3.88  3.88  0.35
9214290       1 25.5   48.21   12.31    35.9   3.88  3.88  0.35
9214390     1.1 58.3   49.67   28.95   20.72   3.88  3.88  0.35
9214490     1.2 65.5   57.28   37.52   19.76   3.88  3.88  0.35
9214590     1.3 75.1   53.84   40.44    13.4   3.88  3.88  0.35
9214690     1.4 58.9   40.41   23.79   16.62   3.88  3.88  0.35
9214790     1.5 76.6    50.1   38.38   11.73  3.88  3.88  0.35
9214890     1.6 72.4   52.78   38.21   14.57   3.88  3.88  0.35
9214990     1.7 66.5   41.27   27.43   13.84   3.88  3.88  0.35
9215090     1.8 62.4    52.4   32.68   19.72   3.88  3.88  0.35
9215190     1.9 43.9   45.11    19.8   25.31   3.88  3.88  0.35
9215290       2 47.1   46.17   21.75   24.42   3.88  3.88  0.35

#convert binary Aux1--Aux3 into decimal
x<- (d$Aux1>mean(range(d$Aux1)))*1 + (d$Aux2>mean(range(d$Aux2)))*2
+(d$Aux3>mean(range(d$Aux3)))*4
x<-abs(x-3)  #do this because resp to +contrast is same as to -contrast

hbr<-d$DeOxyA

#detrend() is my own function which detrends :-)
hbr<-detrend(hbr)

#lm.lag is my own function which shifts x in time so it produces the best
fit with y
#then fits lm(y~(shifted x))
#This treats x as continuous (which it is)
fit1<-lm.lag(x,hbr,lag=c(50,200)); summary(fit1$fit); fit1$lag

xs<-shift(x,lag=97)
plot(0:3, aggregate(cbind(hbr,xs), list(xs), mean)[,2]);
abline(a=18.00991,b= -.18183)

#this is the simple indicator var approach that compares 0 vs 1, 0 vs 2, 0
vs 3
# It is NOT what I want
fit2<-lm(hbr~factor(xs),contrasts=contr.treatment(3))
summary(fit2)


I am wondering, first of all, if an lm() approach can tell me what I want to
know. Namely, the coefficients of the linear fit for y~x, and also some
coefficients which indicate what kind of nonlinearities I have and how
strong they are. I was asking how to do it and threw out lm() as a
possibility. I am open to anything.

Thanks for any help
Bill

	[[alternative HTML version deleted]]

