3
$\begingroup$

I'm trying to build a model to forecast direct mail marketing campaign responses. In the "response" vector are the average number of responses from a marketing campaign from day 1 to day 63 (8 weeks). Before I took the average, I first normalized the individual campaign responses to adjust for seasonality (day-of-the-week, holidays, promotions, etc.) and lined up the campaigns from the first day to the last. My goal is to model this "response time series" with a smooth curve (i.e. continuous probability) via a Weibull curve. Then I want to differentiate between days to get the % of total responses that I'll get in a certain day. In the past I've had success using the Weibull curve to model data that takes on this "response curve shape" but this time it's not working well. Any suggestions?

Here's my code:

#The "responses" vector contains an average number of responses to a Direct Mail marketing campaign from day 1 to day 63 (8 weeks)
#these repsonses have been normalized to account for day-of-the-week variability, holidays, promotions, etc., then lined up by the first starting day

responses <- c(
24.16093706,
41.59607507,
68.20083052,
85.19109064,
100.0704403,
58.6600221,
86.08475816,
88.97439581,
65.58341418,
49.25588053,
53.63602085,
47.03620672,
29.71552264,
32.85862747,
31.29118096,
23.67961069,
19.81261675,
18.69300933,
17.25738435,
12.01161679,
12.36734071,
14.32360673,
11.02390849,
9.108021409,
9.647965622,
8.815576548,
5.67225654,
5.739220185,
6.233999138,
5.527376627,
5.024065761,
5.565266355,
4.626749364,
3.480761716,
4.621902301,
4.518554271,
4.075985188,
3.204946787,
3.174020873,
2.966915873,
2.129178828,
2.673009031,
2.410429043,
2.331287075,
2.509300578,
2.13820695,
2.53433787,
1.603934405,
1.555813592,
1.834605068,
1.842905685,
1.454045577,
2.08684322,
1.318276487,
0.807666643,
1.333167088,
1.004526525,
1.180110123,
1.078079735,
1.151394678,
1.426747942,
0.699119833,
0.583347236)




set.seed(2)
install.packages("MASS")
library("MASS")


shape_and_scale <- fitdistr(responses,'weibull')

#check the shape and scale
shape_and_scale


#now use the curve () function along with dweibull () and the shape/scale factors to fit a smooth curve. 
#to look at the curve from the interval from 0 to 63, set 0 and 63 as parameters 
?dweibull 
curve_results <- curve(dweibull(x,shape_and_scale$estimate[1],shape_and_scale$estimate[2]),from=0, to=63, n=500)
#the above only works if the first argument of dweibull is x (because it's nested in the curve function and it's passing 0 to 63 through). 
#...in order to run dweibull() as a stand alone function you'd need to assign a vector of numeric values as its first argument, or explicitly state a value
#the more "n's" you choose, the tighter the fit. 

curve_results
#by default the curve() function evaluates at 101 points. This is why it increases in increments of .63 instead of 1. 

#notice that when you integrate the dweibull function from 0 to 63 the area (i.e. density) under this part of the curve is only .94.
integrate(dweibull, 0, 63, shape_and_scale$estimate[1],shape_and_scale$estimate[2])
#in order to have the density equal 1, you'd need to integrate from 0 to Inf
#integrate(dweibull, 0, Inf, shape = 0.70730466,scale = 13.79467490)


#These are the densities from 0 to 63
diff(pweibull(0:63,shape_and_scale$estimate[1],shape_and_scale$estimate[2]))
#notice that they sum to .94
sum(diff(pweibull(0:63,shape_and_scale$estimate[1],shape_and_scale$estimate[2])))


#now we want to scale the densities so that represent a % of the total area under the weibull distribution from 0 to 63
diff(pweibull(0:63,shape_and_scale$estimate[1],shape_and_scale$estimate[2]))

scaled_densities <- diff(pweibull(0:63,shape_and_scale$estimate[1],shape_and_scale$estimate[2]))/sum(diff(pweibull(0:63,shape_and_scale$estimate[1],shape_and_scale$estimate[2])))
sum(scaled_densities)

Compare the actual response data with the smooth curve created using the weibull () function.

Actual Responses

Weibull output

Comparison:

enter image description here


After @tristen suggested that I need to make my "responses" vector the actual observations (as opposed to a histogram of the total by day), this is the new "responses" vector I should be using. However, as you'll see, the curve is still not capturing the height of the actual data:

responses <- c(
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  1,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  2,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  4,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  5,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  6,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  7,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  8,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  9,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  10,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  11,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  12,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  13,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  14,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  15,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  16,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  17,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  18,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  19,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  20,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  21,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  22,
  23,
  23,
  23,
  23,
  23,
  23,
  23,
  23,
  23,
  23,
  23,
  24,
  24,
  24,
  24,
  24,
  24,
  24,
  24,
  24,
  25,
  25,
  25,
  25,
  25,
  25,
  25,
  25,
  25,
  25,
  26,
  26,
  26,
  26,
  26,
  26,
  26,
  26,
  26,
  27,
  27,
  27,
  27,
  27,
  27,
  28,
  28,
  28,
  28,
  28,
  28,
  29,
  29,
  29,
  29,
  29,
  29,
  30,
  30,
  30,
  30,
  30,
  30,
  31,
  31,
  31,
  31,
  31,
  32,
  32,
  32,
  32,
  32,
  32,
  33,
  33,
  33,
  33,
  33,
  34,
  34,
  34,
  35,
  35,
  35,
  35,
  35,
  36,
  36,
  36,
  36,
  36,
  37,
  37,
  37,
  37,
  38,
  38,
  38,
  39,
  39,
  39,
  40,
  40,
  40,
  41,
  41,
  42,
  42,
  42,
  43,
  43,
  44,
  44,
  45,
  45,
  45,
  46,
  46,
  47,
  47,
  47,
  48,
  48,
  49,
  49,
  50,
  50,
  51,
  51,
  52,
  53,
  53,
  54,
  55,
  56,
  57,
  58,
  59,
  60,
  61,
  62,
  63
)

However, the curve is still not capturing the height of my original data (see picture below). Any thoughts?

enter image description here


@tristen suggested that I use a log-normal distribution as opposed to a Weibull distribution which essentially solved my problem. BIG THANK YOU! See below for the finalized product:

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Your post is very unclear. What are the numbers in "response" and how were they obtained? Even your first plot doesn't make sense to me (at least not yet). It seems to be a plot of values against the observation number. What values is it showing? Why are they plotted against the numbers 1 to 63? Are you fitting a density to a single variable, or a functional curve to a y-x relationship? $\endgroup$ – Glen_b Apr 14 '15 at 0:40
  • $\begingroup$ Sorry, let me clear up my post- the numbers contained in the "response" vector are an average number of responses to a marketing campaign from day 1 to day 63 (8 weeks). Before I took the average, I normalized the individual campaign responses to adjust for seasonality (day-of-the-week, holidays, promotions, etc.). I'll update my original post. $\endgroup$ – Ryan Chase Apr 14 '15 at 16:14
  • $\begingroup$ So these are averages across different campaigns? $\endgroup$ – Glen_b Apr 14 '15 at 22:50
  • $\begingroup$ yes. First I "de-seasonalized" the responses so that normal day-of-the-week variability doesn't get mistaken for what a normal "response lag" looks like. The reason I averaged was to get a generalized response lag curve so that it wasn't defined based on only one campaign. $\endgroup$ – Ryan Chase Apr 14 '15 at 23:15
  • 1
    $\begingroup$ @tristan: I wouldn't ignore your comment. A gamma curve fits well, including the height and a longer tail. And a plot of the density of the log looks normalish. (I'm not sure where curve-fitting and modeling part ways, though.) $\endgroup$ – Wayne Apr 21 '15 at 14:01
3
$\begingroup$

Your current usage of MASS fitdistr is incorrect, it expects to be supplied with actual observations rather than essentially a histogram which you have supplied. Your best bet I think is to simulate a dataset with 24 values 0.5, 42 values 1.5, 68 values 2.5 etc and then run fitdistr on that. You can get lognormal or log-logistic or gamma by replacing 'weibull' in your call as appropriate.

EDIT

I assumed that people would not all respond exactly at the end of day 1 but would be spread out across the day with average response roughly half way through the day.

You will presumably then want to use the model to predict how many responses will arrive at time $T > t$ given you have received $m(t)$ responses up to time $t$. You would calculate this as

$$m(t) \frac{\Pr(T > t)}{\Pr(T \le t)}= \frac{m(t)\exp\left\{-\left(\frac{t}{b}\right)^a\right\}}{1-\exp\left\{-\left(\frac{t}{b}\right)^a\right\}}$$

Where $a$ is the shape parameter and $b$ the scale parameter.

$\endgroup$
  • $\begingroup$ ah, that makes a lot of sense. Can you explain why I'd use 0.5 for the first 24 values and not just 1 to represent the responses on day 1? $\endgroup$ – Ryan Chase Apr 16 '15 at 1:40
  • $\begingroup$ @RyanChase I have given my reasoning in an edit to my question $\endgroup$ – tristan Apr 16 '15 at 8:48
  • $\begingroup$ I have done what you insisted and changed my responses vector to the actual observations and added it to the bottom of my initial question (is there a better way on StackExchange or Crossvalidated to show the updates in the code you've made since the original post? ...the comments section doesn't give me enough space to add the entire "resposes" vector. I'm still getting an issue thought where the curve does not quite get to the "height" of the actual data. I included a chart showing the new results. Any thoughts? Thank you! $\endgroup$ – Ryan Chase Apr 21 '15 at 13:14
  • $\begingroup$ @RyanChase I took your data and also fit log-normal and gamma distributions. I then used AIC(fit.lognormal), etc. to compare model fits (note that this assumes that the number of data points you have supplied is actually representative of the amount of information you have). I found that log-normal was superior to weibull and gamma and verified this graphically. I suggest you need to pick a different distribution to Weibull and ideally do further model fit checks. $\endgroup$ – tristan Apr 21 '15 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.