4
$\begingroup$

I have the following data (prop is like empirical CDF):

td <- data.frame(a = 3:14, prop=c(0, 0, 0.026, 0.143, 0.21, 0.361, 0.535, 
                                  0.719, 0.814, 0.874, 0.950, 0.964))

I want to fit a normal CDF using an appropriate mean and standard deviation. How can I find such a mean and a standard deviation value that fits the data above?

mean <- 8.8 # <-- How can I find a best fitting number?
sd <- 2.3 # <-- How can I find a best fitting number?
x <- seq(from = 2, to = 15, by = .1)
cdf <- data.frame(x = x, y = pnorm(q = x, mean = mean,  sd = sd))
library(ggplot2)
ggplot(data = td, aes(x = a, y = prop)) + geom_point() +
  geom_line(data = cdf, aes(x = x, y = y))

enter image description here

EDIT:

The motivation behind this question is to replicate a graph I saw on a book. The book used the same proportions and fitted the normal ogive to the data. It looks like normal ogive fitted so well and I couldn't replicate it. There is no way the author used the raw data because the data is from a 100 years old book and the author's book published about 17 years ago. Here is the graph with its caption:

enter image description here

Caption: "Proportions correct on item 46, plotted against age, with a fitted normal ogive."

$\endgroup$
10
  • $\begingroup$ Could you describe the raw data from which you computed the proportions in the td data frame? $\endgroup$ Commented Sep 8, 2019 at 17:53
  • $\begingroup$ Unfortunately, the raw data is not available. a represents the age of a child and prop represents the percentage of children who correctly answered a question in that age group. I believe the total sample size is about 3,000. $\endgroup$
    – HBat
    Commented Sep 8, 2019 at 18:02
  • 1
    $\begingroup$ You need the raw data. If you have the number of children n in each group you can fit a binomial regression model by doing glm(prop ~ a, binomial(link="probit"), weight=n). Using the default logit link function instead of the probit would probably make the results easier to interpret (in terms of oddsratios). $\endgroup$ Commented Sep 8, 2019 at 18:11
  • $\begingroup$ Obviously, having raw data would be great. But, I prefer an answer that uses the given information only, i.e. proportions. I specifically want probit, I'm not interested in any other distribution even if they fit better. $\endgroup$
    – HBat
    Commented Sep 8, 2019 at 18:21
  • 1
    $\begingroup$ One potential problem I can see is that the data is not an empirical CDF: The first two proportions are 0. $\endgroup$ Commented Sep 9, 2019 at 5:45

2 Answers 2

6
$\begingroup$

Edit

@JarleTufto is right to suggest a binomial glm with a probit link.

library(ggplot2)

dat <- data.frame(a = 3:14, prop=c(0, 0, 0.026, 0.143, 0.21, 0.361, 0.535, 
                                   0.719, 0.814, 0.874, 0.950, 0.964))

mod <- glm(prop~a, family = binomial(link = "probit"), data = dat)

xx <- seq(min(dat$a), max(dat$a), length.out = 1000)
pred_frame <- data.frame(a = xx)

pred_frame$fitted <- predict(mod, newdata = pred_frame, type = "response")    

theme_set(theme_bw())
ggplot(data = dat, aes(x = a, y = prop)) + geom_point(size = 3) +
  geom_line(data = pred_frame, aes(x = a, y = fitted), size = 1, colour = "steelblue")

glmprobit

Original answer

I mean you could just minimze the sum of squared residuals (least square fit). The best fitting normal distribution has a mean of $8.826$ and a standard deviation of $2.396$. The R code:

library(ggplot2)

dat <- data.frame(a = 3:14, prop=c(0, 0, 0.026, 0.143, 0.21, 0.361, 0.535, 
                                  0.719, 0.814, 0.874, 0.950, 0.964))

foo <- function(parms, x, y) {
  sum((pnorm(x, mean = parms[1], sd = parms[2]) - y)^2)
}

fit <- optim(c(8, 2), fn = foo, x = dat$a, y = dat$prop)

$par
[1] 8.826028 2.396412

xx <- seq(min(dat$a), max(dat$a), length.out = 1000)
fitted <- pnorm(xx, fit$par[1], fit$par[2])

dat2 <- data.frame(a = xx, prop = fitted)

theme_set(theme_bw())
ggplot(data = dat, aes(x = a, y = prop)) + geom_point(size = 3) +
  geom_line(data = dat2, size = 1, colour = "steelblue")

CDFfit

I'm unsure if that makes much sense, though.

$\endgroup$
2
  • $\begingroup$ Visually, it better fits better compared the solution I proposed. Much closer to the graph I saw on the book. Still wonder why fitdist I proposed did not work as well as this one (visually I mean). $\endgroup$
    – HBat
    Commented Sep 9, 2019 at 6:28
  • 2
    $\begingroup$ -1: Surely, even though the underlying count data is unavailable, you would get a more efficient estimate by running a glm taking into account the smaller variance of the proportions that are close to 0 or 1. $\endgroup$ Commented Sep 9, 2019 at 7:36
0
$\begingroup$

I found a possible solution. I feel that there might be a better solution/fit and open to other suggestions.

# Create data from proportions, since proportions 
# are accurate to 3 digits, I use 1000:
freq <- td$prop * 1000
freq <- c(freq[1], freq[-1] - freq[-length(freq)])
freq <- rep(x = 3:14, freq)

# Fit normal distribution to the data
library(fitdistrplus)
fitdist(freq, 'norm')

Here are the fitted values:

Fitting of the distribution ' norm ' by maximum likelihood 
Parameters:
     estimate Std. Error
mean     9.20     0.0698
sd       2.17     0.0493

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.