[R] trigonometric regression

William Dunlap wdunlap at tibco.com
Thu Jun 17 18:12:07 CEST 2010


> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of Duncan Murdoch
> Sent: Thursday, June 17, 2010 3:19 AM
> To: William Simpson
> Cc: r-help at r-project.org
> Subject: Re: [R] trigonometric regression
> 
> William Simpson wrote:
> > Suppose I do a trigonometric regression
> > fit<-lm(y~ cf + sf)
> > where cf and sf are the cos and sine components.
> >
> > b<-coef(fit)
> > I have the fitted sine component b[2] and the cos component b[3].
> > Doing summary(fit) gives me the p-values and SEs for b[2] and b[3].
> >
> > But I want the amplitude of the fitted waveform
> > amp<-sqrt(b[2]^2+b[3]^2)
> >
> > Can someone please tell me how to get the p-value for amp?
> >   
> 
> "the p-value for amp" is ambiguous; p-values refer to tests, not 
> functions.  But let's assume you want to test whether amp = 0.  Then 
> this is achieved by an F test comparing the model with cf and 
> sf versus 
> one without it.  You'll see it in summary(fit) at the bottom of the 
> display.  If you want to include other covariates in the 
> model, you can 
> use anova, e.g.
> 
> anova(lm(y ~ other), lm(y ~ cf + sf + other))

You can also define a function that keeps the cos
and sin terms together so anova(fit) shows
one entry for the (cos,sin) pair.  E.g., define
the following function
 cs <- function(x, freq)cbind(cos=cos(x*freq), sin=sin(x*freq))
and use it as in
 > time<- sort(runif(30,0,20))
 >
y<-sin(1.7*time+.1)+2*sin(1.1*time+.8)+rnorm(length(time),mean=2,sd=0.3)
 > fit <- lm(y~cs(time,1.7)+cs(time,1.1)+cs(time,1.3))
 > anova(fit)
 Analysis of Variance Table 

 Response: y
               Df Sum Sq Mean Sq  F value    Pr(>F)    
 cs(time, 1.7)  2 12.000  5.9998  50.7267 3.692e-09 ***
 cs(time, 1.1)  2 57.355 28.6776 242.4609 3.493e-16 ***
 cs(time, 1.3)  2  0.780  0.3902   3.2988   0.05501 .  
 Residuals     23  2.720  0.1183                       
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 

> 
> Duncan Murdoch
> 
> ______________________________________________
> 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.
> 



More information about the R-help mailing list