Beta Renan
pkg <- c("ggplot2", "knitr","minpack.lm")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## ggplot2 knitr minpack.lm
## TRUE TRUE TRUEsetwd("/home/epi/Dropbox/MyR/Análise otros/Renan")# lab
#setwd("/media/DATA/Dropbox/MyR/Análise otros/Renan") # Dell
load("~/Dropbox/MyR/Análise otros/Renan/mono renan.RData")Dataset
# Dados planilha laboratorio
dat = read.csv("dat.csv", dec=",", header=T, sep="\t", check.names=FALSE)
head(dat)## exp hm temp rep germ
## 1 1 6 5 1 0.0000000
## 2 1 6 5 2 0.0000000
## 3 1 6 5 3 0.0000000
## 4 1 6 5 4 0.0000000
## 5 1 6 10 1 0.2222222
## 6 1 6 10 2 0.2608696names(dat) = c("exp", "hm", "x", "rep", "y")
dat$symb <- c(1,19)[match(dat$exp, c(1,2))]
head(dat)## exp hm x rep y symb
## 1 1 6 5 1 0.0000000 1
## 2 1 6 5 2 0.0000000 1
## 3 1 6 5 3 0.0000000 1
## 4 1 6 5 4 0.0000000 1
## 5 1 6 10 1 0.2222222 1
## 6 1 6 10 2 0.2608696 1explot = ggplot(dat, aes(x, y, group =exp)) +
theme_bw() +
geom_point(shape=exp, size=3) +
geom_smooth(color="red")+
facet_wrap(~hm) +
coord_cartesian(xlim = c(0, 38))
explot## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
Criamos a função de beta generalizada
beta.reg<-function(x, Yopt,Topt,Tmax, b1) {
ifelse(x>Tmax,0,
#ifelse(x<Tmin,0,
Yopt*((x-5)/(Topt-5))^(b1*(Topt-5)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
)
}6 horas de molhamento
dat6 = subset(dat, hm=="6")
head(dat6)## exp hm x rep y symb
## 1 1 6 5 1 0.0000000 1
## 2 1 6 5 2 0.0000000 1
## 3 1 6 5 3 0.0000000 1
## 4 1 6 5 4 0.0000000 1
## 5 1 6 10 1 0.2222222 1
## 6 1 6 10 2 0.2608696 1grid6 <- expand.grid(Yopt=0.85,
Topt=seq(18,22,len=10),
#Tmin= 5,
Tmax= seq(33,35,len=10),
b1 = seq(1,10,len=10))
mod.lst6 <- apply(grid6,1,function(gr){
nlsLM(y ~ beta.reg(x, Yopt,Topt,Tmax, b1), data=dat6,
start=c(gr),control=nls.control(maxiter=800)) })
rss6 <- sapply(mod.lst6,function(m)sum(residuals(m)^2))
mod6 <- mod.lst6[[which.min(rss6)]] # fit with lowest RSS
coef(summary(mod6)) ## Estimate Std. Error t value Pr(>|t|)
## Yopt 0.8482705 0.01857049 45.678411 1.202682e-43
## Topt 21.7681865 0.33240794 65.486361 1.214838e-51
## Tmax 30.4609039 0.22445312 135.711650 5.377749e-68
## b1 0.8845041 0.13706307 6.453263 3.679673e-08sum(residuals(mod6)^2)## [1] 0.2200961with(dat6, c(plot(y~x, pch=symb)))## NULLpoints(fitted(mod6)~I(dat6$x), pch=4)
with(as.list(coef(mod6)),curve(beta.reg(x, Yopt,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod6))
qqline(residuals(mod6))
12 horas de molhamento
dat12 = subset(dat, hm=="12")
dat12$y = dat12$y + 0.001
head(dat12)## exp hm x rep y symb
## 57 1 12 5 1 0.0010000 1
## 58 1 12 5 2 0.0010000 1
## 59 1 12 5 3 0.0010000 1
## 60 1 12 5 4 0.0010000 1
## 61 1 12 10 1 0.4642035 1
## 62 1 12 10 2 0.3213463 1grid12 <- expand.grid(Yopt = 0.91,
#Tmin = seq(-12, 0, len = 10),
Topt = 20,
Tmax = seq(34, 37, len = 10),
b1 = seq(3, 4, len = 10))
mod.lst12 <- apply(grid12,1,function(gr){nlsLM(y ~ beta.reg(x, Yopt,Topt,Tmax, b1), data=dat12,
start=c(gr),control=nls.control(maxiter=800)) })
rss12 <- sapply(mod.lst12,function(m)sum(residuals(m)^2))
mod12 <- mod.lst12[[which.min(rss12)]] # fit with lowest RSS
coef(summary(mod12)) ## Estimate Std. Error t value Pr(>|t|)
## Yopt 0.8543953 0.02299932 37.148725 4.082527e-39
## Topt 20.9190313 0.39993688 52.305832 1.215251e-46
## Tmax 30.2110944 0.16347300 184.807855 5.910371e-75
## b1 0.7997566 0.13388465 5.973475 2.120313e-07sum(residuals(mod12)^2)## [1] 0.3432969with(dat12, c(plot(y~x, pch=symb)))## NULLpoints(fitted(mod12)~I(dat12$x), pch=4)
with(as.list(coef(mod12)),curve(beta.reg(x, Yopt,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod12))
qqline(residuals(mod12))
24 horas de molhamento
dat24 = subset(dat, hm=="24")
head(dat24)## exp hm x rep y symb
## 113 1 24 5 1 0.0000000 1
## 114 1 24 5 2 0.0000000 1
## 115 1 24 5 3 0.0000000 1
## 116 1 24 5 4 0.0000000 1
## 117 1 24 10 1 0.6228070 1
## 118 1 24 10 2 0.4692982 1grid24 <- expand.grid(Yopt=0.9,
Topt=19,
Tmax= seq(30,35,len=10),
b1 = 0.9)
mod.lst24 <- apply(grid24,1,function(gr){
nlsLM(y ~ beta.reg(x, Yopt,Topt,Tmax, b1), data=dat24,
start=c(gr),control=nls.control(maxiter=800)) })
rss24 <- sapply(mod.lst24,function(m)sum(residuals(m)^2))
mod24 <- mod.lst24[[which.min(rss24)]] # fit with lowest RSS
coef(summary(mod24)) ## Estimate Std. Error t value Pr(>|t|)
## Yopt 0.9290584 0.03235545 28.714124 1.442615e-33
## Topt 19.8192428 0.45874153 43.203507 2.024978e-42
## Tmax 30.6587900 0.47262414 64.869285 1.976022e-51
## b1 0.9845703 0.21420385 4.596417 2.776770e-05sum(residuals(mod24)^2)## [1] 0.6694874with(dat24, c(plot(y~x, pch=symb)))## NULLpoints(fitted(mod24)~I(dat24$x), pch=4)
with(as.list(coef(mod24)),curve(beta.reg(x, Yopt,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod24))
qqline(residuals(mod24))
36 horas de molhamento
dat36 = subset(dat, hm=="36")
head(dat36)## exp hm x rep y symb
## 169 1 36 5 1 0.0000000 1
## 170 1 36 5 2 0.0000000 1
## 171 1 36 5 3 0.0000000 1
## 172 1 36 5 4 0.0000000 1
## 173 1 36 10 1 0.7955556 1
## 174 1 36 10 2 0.7066667 1grid36 <- expand.grid(Yopt=0.8,
Topt=seq(18,22,len=4),
Tmax= seq(30,35,len=10),
b1 = seq(1,10,len=10))
mod.lst36 <- apply(grid36,1,function(gr){
nlsLM(y ~ beta.reg(x, Yopt,Topt,Tmax, b1), data=dat36,
start=c(gr),control=nls.control(maxiter=800)) })
rss36 <- sapply(mod.lst36,function(m)sum(residuals(m)^2))
mod36 <- mod.lst36[[which.min(rss36)]] # fit with lowest RSS
coef(summary(mod36)) ## Estimate Std. Error t value Pr(>|t|)
## Yopt 0.8736399 0.02692609 32.445852 3.493950e-36
## Topt 18.5995730 0.66188867 28.100757 4.147665e-33
## Tmax 30.3486431 0.33333568 91.045287 5.096160e-59
## b1 0.4701389 0.13807139 3.405042 1.282108e-03sum(residuals(mod36)^2)## [1] 0.4834685with(dat36, c(plot(y~x, pch=symb)))## NULLpoints(fitted(mod36)~I(dat36$x), pch=4)
with(as.list(coef(mod36)),curve(beta.reg(x, Yopt,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod36))
qqline(residuals(mod36))