Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Sine curve fit using lm and nls in R

I am a beginner in curve fitting and several posts on Stackoverflow really helped me.

I tried to fit a sine curve to my data using lm and nls but both methods show a strange fit as shown below. Could anyone point out where I went wrong. I would suspect something to do with time but could not get it right. My data can be accessed from here. plot

data <- read.table(file="900days.txt", header=TRUE, sep="")
time<-data$time
temperature<-data$temperature

#lm fitting
xc<-cos(2*pi*time/366)
xs<-sin(2*pi*time/366)
fit.lm<-lm(temperature~xc+xs)
summary(fit.lm)
plot(temp~time, data=data, xlim=c(1, 900))
par(new=TRUE)
plot(fit.lm$fitted, type="l", col="red", xlim=c(1, 900), pch=19, ann=FALSE, xaxt="n",
yaxt="n")

#nls fitting
fit.nls<-nls(temp~C+alpha*sin(W*time+phi),
   start=list(C=27.63415, alpha=27.886, W=0.0652, phi=14.9286))
summary(fit.nls)
plot(fit.nls$fitted, type="l", col="red", xlim=c(1, 900), pch=19, ann=FALSE, xaxt="n", 
axt="n")
like image 889
Eddie Avatar asked Nov 20 '13 19:11

Eddie


1 Answers

This is because the NA values are removed from the data to be fit (and your data has quite a few of them); hence, when you plot fit.lm$fitted the plot method is interpreting the index of that series as the 'x' values to plot it against.

Try this [note how I've changed variable names to prevent conflicts with the functions time and data (read this post)]:

Data <- read.table(file="900days.txt", header=TRUE, sep="")
Time <- Data$time 
temperature <- Data$temperature

xc<-cos(2*pi*Time/366)
xs<-sin(2*pi*Time/366)
fit.lm <- lm(temperature~xc+xs)

# access the fitted series (for plotting)
fit <- fitted(fit.lm)  

# find predictions for original time series
pred <- predict(fit.lm, newdata=data.frame(Time=Time))    

plot(temperature ~ Time, data= Data, xlim=c(1, 900))
lines(fit, col="red")
lines(Time, pred, col="blue")

This gives me:

enter image description here

Which is probably what you were hoping for.

like image 53
Andy Barbour Avatar answered Oct 25 '22 05:10

Andy Barbour