R help - Jan 2003 - (v2) quadratic trends and changes in slopes (R-help digest, Vol 1 #52 - 16 msgs)

-----Original Message-----
 Message: 6
Date: Mon, 20 Jan 2003 01:11:24 +0100
From: Martin Michlmayr <tbm at cyrius.com>
To: r-help at stat.math.ethz.ch
Subject: [R] quadratic trends and changes in slopes

I'd like to use linear and quadratic trend analysis in order to find
out a change in slope.  Basically, I need to solve a similar problem as
discussed in
http://www.gseis.ucla.edu/courses/ed230bc1/cnotes4/trend1.html

....

RESPONSE: This message updates my earlier response. I just started
programming with R on Friday and I wanted to make my response to you
more general... before starting on a more difficult problem of the same
type at work. The two improvements I've added are: (1) a more general
procedure for more quickly and accurately creating the contrast vectors
and (2) direct program interaction with results in output tables. Access
to results in output tables is so much easier in R than what I've had to
do in SAS that it's incredible. An R program to work the example you
cited is appended. I hope this meets your needs.

Chuck White

# R program for working example at:
# http://www.gseis.ucla.edu/courses/ed230bc1/cnotes4/trend2.html
#
# Copy and paste the example data from the web to a plain text editor
# and save as 'example.txt'. If this program is copied and passted into
# the R Console then the program is expected to run without further
# operator intervention.

#Read the data as follows:

example<-read.table("example.txt", header = FALSE, 
col.names=c('y','grp','o1','o2','o3'))
example

# Make variable names available to session

attach(example)

# Convert grp from numeric format to factor format for ANOVA. If you
don't 
# do this then you get the wrong test (with only 1 degree of freedom).

fgrp<-factor(grp)
fgrp

# Conduct ANOVA on grp

GRP<-lm(y~fgrp)
anova(GRP)

# Conduct ANOVA on contrasts

## The example aready contains the contrast vectors but you'll have to 
## create them yourself when you use real data. A quick and accurate way
to ## set them up is to use the repeat command (rep) as follows:

Linear<-c(rep(-3,8),rep(-1,8),rep(1,8),rep(3,8))
Quadratic<-c(rep(1,8),rep(-1,16),rep(1,8))
Cubic<-c(rep(-1,8),rep(3,8),rep(-3,8),rep(1,8))

Contrasts<-lm(y~Linear+Quadratic+Cubic)
Contrasts.anova<-anova(Contrasts)
Contrasts.anova

# Calculate the F-test and P-value for Nonlinear

SS.Quadratic<-Contrasts.anova$"Sum Sq"[2]
SS.Cubic<-Contrasts.anova$"Sum Sq"[3]
SS.Nonlinear<-SS.Quadratic+SS.Cubic

DF.Quadric<-Contrasts.anova$"Df"[2]
DF.Cubic<-Contrasts.anova$"Df"[3]
DF.Nonlinear<-DF.Quadric+DF.Cubic

SS.Residuals<-Contrasts.anova$"Sum Sq"[4]
DF.Residuals<-Contrasts.anova$"Df"[4]

Nonlinear.test<-(SS.Nonlinear/DF.Nonlinear)/(SS.Residuals/DF.Residuals)

Nonlinear.test

Nonlinear.pvalue<-1-pf(Nonlinear.test,DF.Nonlinear,DF.Residuals)

Nonlinear.pvalue

On 20 Jan 2003 at 21:49, Chuck White wrote:
> 
> I'd like to use linear and quadratic trend analysis in order to find
> out a change in slope.  Basically, I need to solve a similar problem as
> discussed in
> http://www.gseis.ucla.edu/courses/ed230bc1/cnotes4/trend1.html
> 
This response show how to do the test of non-linearity in a 
complicated way, all can be done much easier in R, start to 
look at poly() and contr.poly() (and summary.aov with the argument 
split=). But that is not the point. The original poster did'nt want 
to test for nonlinearity, he assumed there is nonlinearity and wanted
to estimate the change point. He also said that the usual procedure
to do that in his field is to estimate cuadratic models for data
1, 1:2, 1:3, ..., 1:9 (or some similar number) and take the change-
point as the value of i above (in 1:i) where the quadratic term
first is significant. That cannot be sound, as you obviously must go
 somewhat past the changepoint before the quadratic term can become 
significant! So this method cannot possibly give an consistent 
estimator of the change-point. He should use some other method, like 
building a model with an explicit change-point and estimate that.

Kjetil Halvorsen