Nonlinear least squares regression of skewed normal distribution in R (or any language)

user_13282905 0
Viewed 413

First time poster. Apologize in advance if I use improper etiquette or vocabulary.

I've time series data of chemical concentration (y) vs time (x) from a USGS river survey. It exhibits a skew normal distribution that I would like to model via non-linear least squares regression. I'm able to fit a normal distribution curve to the data, but can't seem to incorporate "skewness" into the model.

I arrived at the my normal distribution fit from the answer given by Whuber here... Linear regression best polynomial (or better approach to use)?

my data and code...

y <- c(0.532431978850729, 0.609737363640599, 0.651964078008195, 0.657368066358271, 
0.741496240155044, 0.565435828629966, 0.703655525439792, 0.718855614453251, 
0.838983191559565, 0.743767469276213, 0.860155614137561, 0.81923941209205, 
1.07899884812998, 0.950877380129941, 1.01284743983765, 1.11717867112622, 
1.08452873942528, 1.14640319037414, 1.35601176845714, 1.55587090166098, 
1.81936731953165, 1.79952819117948, 2.27965075864338, 2.92158756334143, 
3.28092981974249, 1.09884083379528, 4.52126319475028, 5.50589160306292, 
6.48951979830975, 7.61196542128105, 9.56700470248019, 11.0814901164772, 
13.3072954022821, 13.8519364143597, 11.4108376964234, 8.72143939873907, 
5.12221325838613, 2.58106436004881, 1.0642701141608, 0.44945378376047, 
0.474569233285229, 0.128299654944011, 0.432876244482592, 0.445456125461339, 
0.435530646939433, 0.337503495863836, 0.456525976632425, 0.35851011819921, 
0.525854215793115, 0.381206935673774, 0.548351975353343, 0.365384673834335, 
0.418990479166088, 0.50039125911365, 0.490696977485334, 0.376809405620949, 
0.484559448760701, 0.569111550743562, 0.439671715276438, 0.353621820313257, 
0.444241243031233, 0.415197754444015, 0.474852839357701, 0.462144150397257, 
0.535339727332139, 0.480714031175711)

#creating an arbitrary vector to represent time
x <- seq(1,length(y), by=1)

#model of normal distribution 
f <- function(x, theta)  { 
  m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4];
  a*exp(-0.5*((x-m)/s)^2) + b
}

# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y))

# Do the fit.  (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0))

# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]

par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
     xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)

So, any suggestions on how to adjust the model to accommodate skewness?

Cheers, Jamie

skew least-squares normal-distribution regression r
3 Answers

Data are concentration versus time of certain compound in water samples from a river, aren't they? If I plot y vs x, assuming that samples were taken at regular time intervals, I see a peak of concentration, therefore time dependence seems to be some kind of physical and/or chemical phenomenon that could be modeled as y = f(b, x) + e, where f is a function of parameters b of chemical/physical phenomena and x represents time. Term e is the random error, in chemistry usually samples are independently measured thus e ~ N(0, s^2). Then you fit f(b, x) by nls.

user_89218110

Can you use generalized additive model (GAM)? GAM is powerful and flexible, but it is difficult to interpret the model coefficients. So the decision would be depends on your purpose. If the purpose is to evaluate trend, or the purpose is to predict concentration (within the known time range), then GAM could be a good choice.

library(mgcv)
library(ggplot2)

dat <- data.frame(x = 1:length(y), y = y)

fit_gam <- gam(y ~ s(x, k = 20), data = dat) 

ggplot(dat, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = data.frame(x = x, y = fit_gam$fitted.values),
            color = "red") +
  ggtitle("Data") +
  xlab("Cocentration") +
  ylab("Time") +
  theme_bw() +
  theme(panel.grid = element_blank())

enter image description here

The following is another option to apply stat_smooth to fit the same GAM model.

ggplot(dat, aes(x = x, y = y)) +
  geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "tp", k = 20)) +
  ggtitle("Data") +
  xlab("Cocentration") +
  ylab("Time") +
  theme_bw() +
  theme(panel.grid = element_blank())

enter image description here

user_76698090

I spoke to a buddy who's a wiz in python, and he helped me construct the right skewed normal distribution equation. I've posted the R script bellow.

What I was wanting to do was replace the normal distribution model with the right-skewed distribution model. It wasn't so much the script writing that caught me up, but my ability to write a general equation for a right skewed distribution (my buddy is also a math wiz).

I gave the thumbs up to www because for all intensive purposes they answered my question. I like that they took a different approach using a GAM as well, although I am interested in the coefficients produced by the model.

My next plan is to integrate the area under the model curve, as well as the area under the confidence interval curves.

First experience with stackoverflow was good. Thank you all. 

f <- function(x, theta)  { 
  m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4]; k <- theta[5]
  a*exp(k*((x - m))/s - sqrt(((x - m))/s*((x - m))/s+1)) + b
}

# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y)); k.0 <- -0.5

# Do the fit.  (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b, k)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0, k=k.0))

# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]

par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
     xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)
user_132829050
Related