r - How to plot raw data and add predict + confidence intervals using maximum likelihood function in R gpplot2?
问题描述
I'm trying to plot the raw data along with the curve and confidence intervals for maximum likelihood model function in R. My data is stored in df(x.number, y.size)
. I'm trying to run MLL using skewfunc
.
I want to use estimates from the model to plot the curve + 95% confidence intervals along with the raw data, but I am not sure how to proceed:
- It's not a regular polynomial so I can't use the
geom_smooth
function. - The function also has a
ifelse
that complicates matters. - I'm also unclear on how to extract the 95% confidence intervals from the model.
library(stats4)
library(ggplot2)
y.size <- c(2.69,4.1,8.04,3.1,5.27,5.033333333,3.2,7.25,6.29,4.55,6.1,2.65,3.145,3.775,3.46,5.73,5.31,4.425,3.725,4.32,5,3.09,5.25,5.65,3.48,6.1,10,9.666666667,6.06,5.9,2.665,4.32,3.816666667,3.69,5.8,5,3.72,3.045,4.485,3.642857143,5.5,6.333333333,4.75,6,7.466666667,5.03,5.23,4.85,5.59,5.96,5.33,4.92,4.255555556,6.346666667,4.13,6.33,4,7.35,6.35,4.63,5.13,7.4,4.28,4.233333333,4.3125,6.18,4.3,4.47,4.88,4.5,2.96,2.1,3.7,3.62,5.42,3.8,5.5,3.27,3.36,3.266666667,2.265,3.1,2.51,2.51,4.4,2.64,4.38,4.53,2.29,2.87,3.395,3.26,2.77,3.22,4.31,4.73,4.05,3.48,4.8,4.7,3.05,4.21,5.95,4.39,4.55,4.27,4.955,4.65,3.32,3.48,3.828571429,4.69,4.68,3.76,3.91,4,4.41,4.19,4.733333333,4.32,2.83,3.41,4.42,3.47,3.84,4.39)
x.number <- c(69,62,8,80,13,12,2,22,19,49,840,44,31,56,33,58,91,8,15,86,11,69,12,24,32,27,1,4,26,4,28,33,1516,41,20,58,44,29,58,14,3,3,6,3,26,52,26,29,92,30,18,11,27,19,38,78,57,52,17,45,56,7,37,7,14,13,164,76,82,14,273,122,662,434,126,374,1017,522,374,602,164,5,191,243,134,70,23,130,306,516,414,236,172,164,92,53,50,17,22,27,92,48,30,55,28,296,35,12,350,17,22,53,97,62,92,272,242,170,37,220,452,270,392,314,150,232)
df <- data.frame(x.number, y.size)
df <- df[df$x.number < 750,] #data
aa = -1.02 ; K = 12; Ka = 15; q = -0.9; sigma = 0.1
skewfun <- function(aa, K, Ka, q, sigma){
p <- df$x.number / K
lnqp <- if (q == 0) log(p) else (p^q - 1) / q
y.pred <- (aa * (p * K / Ka - 1) - 1) * lnqp
ll <- -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
ll
}
mle2.model <- mle(skewfun, start = list(aa = -1.02, K = 12, Ka = 15, q = -0.9, sigma = 0.1))
summary(mle2.model)
-logLik(mle2.model)
AIC(mle2.model)
ggplot(data=df,aes(x=x.number,y=y.size))+
geom_point() #Raw Data
解决方案
推荐阅读
- python - Py.test:如何在某些值应返回错误时进行参数化
- python-3.x - 如何修复诸如 T&C 之类的标记短语被拆分为 'T' '&' 'C'
- java - JBoss - 如何在多数据库环境中设置默认数据源?
- excel - 有没有办法在UTC格式的excel“一般”单元格中添加一天?
- spring - Spring数据@Transaction没有按顺序执行
- python - 将一个函数的返回值传递给另一个函数
- asp.net - cshtml aspfor getter 不工作,但 setter
- postgresql - 在带有 SELECT 语句的 INSERT INTO 语句中返回
- javascript - 如何进行超过 6 个并发的 REST API 调用
- javascript - 自己的时间选择器不取默认值