Weibull Distribution – How to Fit a Weibull Function in R

rweibull distribution

I am trying to model some data that follows a sigmoid curve relationship. In my field of work (psychophysics), a Weibull function is usually used to model such relationships, rather than probit.

I am trying to create a model using R. My model will have a 'shape' and 'scale' value. In order to generate this model, I would have thought that R would need to have input vectors from both my X and Y axes (I am plotting size of a visual stimulus against the probability of it being seen). However, the R help page gives details of functions that only require either a vector of quantiles (size of stimulus) or a vector of probabilities (presumably, this means my 'probability of seeing' data).

I don't see how R could possibly fit a model on the grounds of only having half of my data (i.e. using my X axis (quantile) data, but none of my Y axis (probability of seeing) data).

Am I missing something? I'm not a statistician so I fear there's something fundamental that I've overlooked. Any help would be appreciated.

Best Answer

What you are trying to do (if I understand you correctly), is to model the probability of the stimulus being seen by a Weibull distribution... but with the Weibull parameters depending on a covariate, i.e., the size of the stimulus. So you are really doing a so-called "Weibull regression". The weibreg() function in the eha package looks like it should do what you want.

Related Solutions

Solved – Multiple regression with categorical and numeric predictors

Try this:

fit <- glm(wealth_indicator ~ 
           factor(ranking) + age_in_years + factor(ranking) * age_in_years)

The factor() command will make sure that R knows that your variable is categorical. This is especially useful if your categories are indicated by integers, otherwise glm will interpret the variable as continuous.

The factor(ranking) * age_in_years term lets R know that you want to include the interaction term.

Solved – Fitting Weibull distribution in R

Here's the fitted pdf and cdf (Weibull) for each of locations 1 to 3:

Let's break down what we need to do here, keeping in mind that the end goal is to estimate the cumulative proportion of area planted with a certain crop at some value for the random variable time $X$:

The first step is to fit a distribution (e.g. $Weibull\left(a:= \text{shape},b := \text{scale}\right)$ or $Beta\left(\alpha,\beta\right)$).

You made an error in fitting the data on a Weibull distribution because the function fitdist (which uses MLE by default) expects a vector of observations, which would be 'time $X_i$ when some area $i$ was planted with a certain crop'.

Here is some R that turns the cumulative data into a vector of observations, which can then be used to fit the distribution with fitdist:

dat <- dat[order(dat$loc.id,dat$year.id),]
dat$loc.id <- factor(dat$loc.id)
dat$year.id <- factor(dat$year.id)

library(plyr)
#hacking the original data from the cumulative
to.density <- function(cdf) {  #given input (a cdf vector), returns respective densities
  y <- c(); cdf <- cdf[,'cum.per.plant']
  for (cum in seq(1,NROW(cdf))) {y <- c(y, ifelse(cum==1, cdf[cum], cdf[cum]-cdf[cum-1]))}
  return(y)
}

#apply to each location X year in dat, then kbind density data to dat
densities <- dlply(dat, .(loc.id, year.id), to.density) 
dat <- cbind(dat,dens=unlist(densities)) 

#forming a dataset of the variable planting time
#--because the proportions go to the thousandths place, a pseudo-sample of time.id's with
#--n=1000 should allow usage of fitdist
rep.vec <- function(df) {
  y <- c()
  for (row in rownames(df)) {y <- c(y, rep(df[row,'time.id'], df[row,'dens']*1000))}
  return(y)
}

#apply to each location in dat - may want to include year.id in the split (excluded)
time.sample <- dlply(dat, .(loc.id), rep.vec)

Be sure to assess the 'fit' of your estimated parameters (perhaps by calculating the mean square error) on the data.

Your code does not indicate that each location and year will be fitted with its own distribution, although you suggest you want to calculate the cumulative area planted for a given location or year.

Here is some R for fitting each location:

#fitting separate Weibull distributions for each loc.id (may want to include year.id in 
#the split)
library(fitdistrplus)
fit.weibull <- function(loc) {
  y <- summary(fitdist(loc,'weibull'))[[1]]
  return(y)
}
params <- lapply(time.sample, fit.weibull) #apply to each element in time sample

Finally, consider the inclusion of a location parameter, which shifts the graph of the pdf in a negative or positive direction along the x-axis; this should be appropriate because in many locations, no area gets plotted until the $X=2^{nd}$ week. The inclusion of a location parameter would likely improve the fit of your distributions.

Use the fitted cdf (with the parameters informed by the previous step) to predict the cumulative proportion of area planted on a certain day for a given location.

It is important to understand what you are doing when you want to rescale the random variable time $X$ to represent days rather than weeks, which simply involves dividing the data (a vector of time observations) by 7. This can be seen in the equation for the Weibull pdf itself.

I've demonstrated it through a few lines of R:

#creating predictive model - on a daily rather than weekly domain
predict1.cum.plant <- function(day, loc.id, params) {
  #'scaling' 7x the scale parameter
  pweibull(day, shape=params[[loc.id]][[1]], scale=params[[loc.id]][[2]]*7)  
}

predict2.cum.plant <- function(day, loc.id, params) {
  #'scaling' 1/7x the random variable time
  pweibull(day/7, shape=params[[loc.id]][[1]], scale=params[[loc.id]][[2]])  
}

#example - the models are equal, which is obvious from the cdf equation!
predict1.cum.plant(35, 2, params)
predict2.cum.plant(35, 2, params) 
pweibull(5, shape=params[[2]][[1]], scale=params[[2]][[2]]) #by week

Some supplemental code of mine can be found here.

Very similar methods can be used to fit a Beta distribution.

Best Answer

Related Solutions

Solved – Multiple regression with categorical and numeric predictors

Solved – Fitting Weibull distribution in R

Related Question