微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

如何估计指数 MLL 函数 + 在 R gpplot2 中添加置信区间?

如何解决如何估计指数 MLL 函数 + 在 R gpplot2 中添加置信区间?

我正在尝试绘制原始数据以及指数曲线拟合,并为 R 中的最大似然模型函数添加置信区间。我的数据存储在 df(x.number,y.size) 中。我正在尝试使用 R 中的 skewfun 运行 MLL

我想使用模型的估计值来绘制曲线 + 95% 置信区间以及原始数据,但我一直遇到错误 (Error in solve.default(oout$hessian) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0)

library(stats4)
library(ggplot2)

#Data
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,10,9.666666667,6.06,5.9,2.665,3.816666667,3.69,5.8,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,3.27,3.36,3.266666667,2.265,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,4.8,4.7,3.05,4.21,5.95,4.39,4.27,4.955,4.65,3.32,3.828571429,4.69,4.68,3.76,3.91,4.41,4.19,4.733333333,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,15,86,11,69,24,32,27,1,26,28,1516,41,20,29,14,3,52,92,30,18,38,78,57,17,45,7,37,164,76,82,273,122,662,434,126,374,1017,522,602,191,243,134,70,23,130,306,516,414,236,172,53,50,48,55,296,35,350,97,272,242,170,220,452,270,392,314,150,232)
df <- data.frame(x.number,y.size)
df <- df[df$x.number < 750,] #data

#Plotting using stat_smooth function
ggplot(data=df,aes(x=x.number,y=y.size))+ geom_point(shape=21,fill="white",color="red",size=3)+ stat_smooth(method="nls",formula =  y~(a*exp(-x*b) + c),method.args=list(start=c(a=10,b=0.05,c=3)),se=F,color="red")+ stat_smooth(method="lm",formula =  y~a*exp(-x*b)+c,color="blue") + theme_classic()

enter image description here

 #MLL function
a = 10; b = 0.05; c = 3; sigma = 1
expreg <- function(a,b,c,sigma){
  y.pred <- a * exp(-x.number*b) + c
  ll <-   -sum(dnorm(y.size,mean = y.pred,sd = sigma,log=TRUE ))
  # cat(a,sigma,ll,"\n")
  ll
}
#Here I keep running into an error that reads
Error in solve.default(oout$hessian) :    Lapack routine dgesv: system is exactly singular: U[1,1] = 0

#Computing MLL

mle2.expreg.model <- mle(expreg,start = list(a = 9,b = 0.03,c = 3,sigma = 1))
summary(mle2.expreg.model)
-logLik(mle2.expreg.model)
AIC(mle2.expreg.model)

#Plot with MLL estimates
a <- 8.66
b <- 12.67
c <- 4.48
ggplot(data=df,size=3)+ stat_function(fun=function(x)a*exp(-x*b) + c,color = "blue")  + My_Theme + theme(panel.background = element_blank()) + ylab("Female Size (mm)") + xlab ("Number of Females/ Host")

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。