[R] Predict with loess
Prof Brian Ripley
ripley at stats.ox.ac.uk
Tue Jun 13 10:58:44 CEST 2006
On Tue, 13 Jun 2006, Jouanin Celine wrote:
> I want to do a nonparametric regression. I’m using the function loess.
> The variable are the year from 1968 to 1977 and the dependant variable is a proportion P. The dependant variable have missing value (NA).
> The script is :
>
> year <- 1969:2002
> length(year)
> [1] 34
>
> P <- c(NA,0.1,0.56,NA,NA,0.5,0.4,0.75,0.9,
> 0.98,0.2,0.56,0.7,0.89,0.3,0.1,0.45,0.46,0.49,0.78,
> 0.25,0.79,0.23,0.26,0.46,0.12,0.56,0.8,0.55,0.41,
> 0.36,0.9,0.22,0.1)
> length(P)
> [1] 34
>
> lo1 <- loess(P~year,span=0.3,degree=1)
> summary(lo1)
> yearCo <- 1969:2002
> year_lo <- data.frame(year = yearCo )
> length(year_lo)
> [1] 34
I get 1 here, and so should you.
> mlo <- predict(loess(P~year,span=0.3,degree=1),new.data=year_lo,se=T)
It should be newdata, not new.data
> mlo$fit
> mlo$se.fit
Notice that these are of length 31, not 34
You are trying to predict at the values used for fitting (possibly not
what you intended), so you don't actually need this. Try
lo1 <- loess(P~year,span=0.3,degree=1, na.action=na.exclude)
fitted(lo1)
plot(year,P,type='o')
lines(year, fitted(lo1))
Or if you want to try interpolation
lines(year, predict(lo1, newdata=year_lo))
This will not extrapolate to 1969, and as far as I recall the version of
loess in R does not allow extrapolation.
> plot(year,P,type='o')
> lines(year,predict(loess(P~year,span=0.15,degree=1),new.data=year_lo,
> se=T,na.action=na.omit)$fit,col='blue',type='l')
>
> The message error indicates that x and y don’t have the same length.
> In fact in m$fit and m$se.fit there are 3 values who don’t have a
> fitted value.
Correct, and that's because you used na.action=na.omit and did not specify
newdata.
[...]
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list