Commit 6cc93d0b authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Aplicações em rpanel para distribuições e reg. não linear.

parent 6ab0dff4
##----------------------------------------------------------------------
## Função que a partir de uma expressão do modelo cria deslizadores para
## estudo do comportamento da função com relação a alterações nos
## valores dos parâmetros.
##----------------------------------------------------------------------
## Definições da sessão.
require(rpanel)
##----------------------------------------------------------------------
## Função.
eyefun <- function(model,
start,
xlim=c(0,1),
ylim=c(0,1),
length.out=101,
...){
## PORÇÃO DE PROTEÇÕES DA FUNÇÃO.
## Nome da variável dependente.
## all.vars(model[[2]])
vardep <- "y"
## Nome da variável independente, testa se é apenas uma.
varindep <- "x"
if (length(varindep)!=1){
stop("just one independent variable is expected!")
}
## Nome dos parâmetros.
parnames <- intersect(all.vars(model[[3]]), names(start))
## Testa se os nomes dos parâmetros seguem o padrão th1, th2, ...
test.start.names <- length(
grep("^th[1-5]$", parnames))==length(parnames)
if (!test.start.names){
stop(paste("in start and model formula parameter names",
"should follow the pattern: th1, th2, ..., th5!"))
}
## Função que converte uma formula em uma função que retorna valores
## da função testa e adequa os nomes dos elementos dos vetores
## dentro da lista start.
startnames <- sapply(start,
function(x){
!(!is.null(names(x)) &
all(names(x)%in%c("init","from","to")))
})
if (any(startnames)){
message(paste("at least one start element is not named.",
"Using the current order to name it as init,",
"from and to."))
for(j in which(startnames)){
names(start[[j]]) <- c("init","from","to")
}
}
## PORÇÃO DE FUNÇÕES INTERNAS.
## Função que converte uma fórmula da nls() em uma função para
## predição.
form2func <- function(formu){
arg1 <- all.vars(formu)
arg2 <- vector("list", length(arg1))
names(arg2) <- arg1
Args <- do.call("alist", arg2)
fmodele <- as.function(c(Args, formu))
return(fmodele)
}
fmodele <- form2func(model[[3]])
## Função que é controlada e passar a curva por meio dos pontos.
nlr.draw <- function(panel){
vindepseq <- seq(xlim[1], xlim[2], length.out=length.out)
listparvar <- c(list(vindepseq), panel[parnames])
names(listparvar)[1] <- varindep
fx <- do.call("fmodele", listparvar)
plot(vindepseq, fx, col=1, lty=1,
xlim=xlim, ylim=ylim, type="l",
...)
panel
}
## PORÇÃO COM OS CONTROLADORES.
action <- nlr.draw
## Abre o painel e as caixas de seleção.
nlr.panel <- rp.control(title="Ajuste",
size=c(200, 200), model=model)
## Abre os deslizadores para controlar o valor dos parâmetros.
if (any(names(start)=="th1")){
rp.slider(panel=nlr.panel, var=th1,
from=start[["th1"]]["from"],
to=start[["th1"]]["to"],
initval=start[["th1"]]["init"],
showvalue=TRUE, action=action)
}
if (any(names(start)=="th2")){
rp.slider(panel=nlr.panel, var=th2,
from=start[["th2"]]["from"],
to=start[["th2"]]["to"],
initval=start[["th2"]]["init"],
showvalue=TRUE, action=action)
}
if (any(names(start)=="th3")){
rp.slider(panel=nlr.panel, var=th3,
from=start[["th3"]]["from"],
to=start[["th3"]]["to"],
initval=start[["th3"]]["init"],
showvalue=TRUE, action=action)
}
if (any(names(start)=="th4")){
rp.slider(panel=nlr.panel, var=th4,
from=start[["th4"]]["from"],
to=start[["th4"]]["to"],
initval=start[["th4"]]["init"],
showvalue=TRUE, action=action)
}
if (any(names(start)=="th5")){
rp.slider(panel=nlr.panel, var=th5,
from=start[["th5"]]["from"],
to=start[["th5"]]["to"],
initval=start[["th5"]]["init"],
showvalue=TRUE, action=action)
}
invisible()
}
##----------------------------------------------------------------------
## Usos da função.
## Essa função foi feita a muito e para simplificar você pode só deve
## usar y e x para indicar as variáveis dependente e independente. Os
## parâmentros devem ser indicados por th seguido de um número de 1 à
## 5. Modelos com mais de 5 parâmentros precisam que a função seja
## modificada.
## Equação da reta: A+B*x.
model <- y~th1+th2*x
start <- list(th1=c(init=0.2, from=0, to=0.5),
th2=c(init=0.6, from=0.4, to=0.8))
eyefun(model, start)
## Modelo beta generalizado.
model <- y~th1*(x-th2)^th3*(th4-x)^th5
start <- list(th1=c(1,0,3),
th2=c(0,0,3),
th3=c(1,0,3),
th4=c(1,0,1),
th5=c(1,0,3))
eyefun(model, start)
##----------------------------------------------------------------------
##----------------------------------------------------------------------
## Definições da sessão.
sourceUrl <-
"https://raw.githubusercontent.com/walmes/wzRfun/master/R/rp.nls.R"
## source(sourceUrl)
download.file(url=sourceUrl,
dest="rp.nls.R")
source("rp.nls.R")
library(rpanel)
library(latticeExtra)
##----------------------------------------------------------------------
## Regressão linear simples.
rp.nls(model=dist~b0+b1*speed,
data=cars,
start=list(
b0=c(init=0, from=-20, to=20),
b1=c(init=2, from=0, to=10)),
assignTo="cars.fit")
cars.fit
##----------------------------------------------------------------------
## Um exemplo mais interessante.
xyplot(rate~conc, groups=state, data=Puromycin)
rp.nls(model=rate~Int+(Top-Int)*conc/(Half+conc),
data=Puromycin,
start=list(
Int=c(init=50, from=20, to=70),
Top=c(init=150, from=100, to=200),
Half=c(init=0.1, from=0, to=0.6)),
subset="state",
assignTo="Puro.fit",
startCurve=list(col=3, lty=3, lwd=1),
fittedCurve=list(col=4, lty=1, lwd=1.5),
extraplot=function(Int, Top, Half){
abline(h=c(Int, Top), v=Half, col=2, lty=2)
},
finalplot=function(Int, Top, Half){
abline(h=c(Int, Top), v=Half, col=3, lty=1)
},
xlab="Concentration",
ylab="Rate",
xlim=c(0, 1.2),
ylim=c(40, 220))
length(Puro.fit)
sapply(Puro.fit, coef)
sapply(Puro.fit, logLik)
sapply(Puro.fit, deviance)
##======================================================================
## MAIS CASOS.
##----------------------------------------------------------------------
## 1. Ajuste de curvas de retenção de água do solo.
cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt",
header=TRUE, sep="\t")
cra$tens[cra$tens==0] <- 0.1
cras <- subset(cra, condi=="LVA3,5")
cras <- aggregate(umid~posi+prof+tens, data=cras, FUN=mean)
cras$caso <- with(cras, interaction(posi, prof))
cras$ltens <- log(cras$tens)
xyplot(umid~ltens|posi, groups=prof, data=cras, type=c("p","a"))
## modelo: van Genuchten com retrição de Mualem.
## x: representado por ltens (log da tensão matricial, psi).
## y: representado por umid, conteúdo de água do solo (%).
## th1: assíntota inferior, mínimo da função, quando x -> +infinito.
## th2: assíntota superior, máximo da função, quando x -> -infinito.
## th3: locação, translada o ponto de inflexão.
## th4: forma, altera a taxa ao redor da inflexão.
model <- umid~Ur+(Us-Ur)/(1+exp(n*(al+ltens)))^(1-1/n)
start <- list(Ur=c(init=0.2, from=0, to=0.5),
Us=c(init=0.6, from=0.4, to=0.8),
al=c(1, -2, 4),
n=c(1.5, 1, 4))
rp.nls(model=model, data=cras,
start=start, subset="caso",
assignTo="cra.fit")
sapply(cra.fit, coef)
lapply(cra.fit, summary)
##----------------------------------------------------------------------
## 2. Curva de produção em função da desfolha do algodão.
cap <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt",
header=TRUE, sep="\t", encoding="latin1")
cap$desf <- cap$desf/100
cap <- subset(cap, select=c(estag, desf, pcapu))
cap$estag <- factor(cap$estag, labels=c("vegetativo","botão floral",
"florescimento","maçã","capulho"))
str(cap)
xyplot(pcapu~desf|estag, data=cap, layout=c(5,1),
xlab="Nível de desfolha artificial", ylab="Peso de capulhos")
## modelo: potência.
## x: representado por desf (nível de desfolha artifical).
## y: representado por pcapu (peso de capulhos), produto do algodão.
## th1: intercepto, valor da função quando x=0 (situação ótima).
## th2: diferença no valor da função para x=0 e x=1 (situação extrema).
## th3: forma, indica como a função decresce, se th3=0 então função
## linear.
model <- pcapu~f0-delta*desf^exp(curv)
start <- list(f0=c(30,25,35), delta=c(8,0,16), curv=c(0,-2,4))
rp.nls(model=model, data=cap,
start=start, subset="estag",
assignTo="cap.fit")
model <- pcapu~f0-f1*desf^((log(5)-log(f1))/log(xde))
start <- list(f0=c(30,25,35), f1=c(8,0,16), xde=c(0.5,0,1))
x11()
rp.nls(model=model, data=cap,
start=start, subset="estag",
assignTo="cap.fit",
extraplot=function(f0,f1,xde){
abline(v=xde, h=c(f0, f0-f1), lty=2, col=2)
})
length(cap.fit)
sapply(cap.fit, coef)
lapply(cap.fit, summary)
##----------------------------------------------------------------------
## 3. Curva de produção em função do nível de potássio no solo.
soja <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
header=TRUE, sep="\t", encoding="latin1", dec=",")
soja$agua <- factor(soja$agua)
str(soja)
xyplot(rengrao~potassio|agua, data=soja)
## modelo: linear-plato.
## x: representado por potássio, conteúdo de potássio do solo.
## y: representado por rengrao, redimento de grãos por parcela.
## f0: intercepto, valor da função quando x=0.
## tx: taxa de incremento no rendimento por unidade de x.
## brk: valor acima do qual a função é constante.
model <- rengrao~f0+tx*potassio*(potassio<brk)+tx*brk*(potassio>=brk)
start <- list(f0=c(15,5,25), tx=c(0.2,0,1), brk=c(50,0,180))
rp.nls(model=model, data=soja,
start=start, subset="agua",
assignTo="pot.fit")
sapply(pot.fit, coef)
##----------------------------------------------------------------------
## 4. Curva de lactação.
lac <- read.table(
"http://www.leg.ufpr.br/~walmes/data/saxton_lactacao1.txt",
header=TRUE, sep="\t", encoding="latin1")
lac$vaca <- factor(lac$vaca)
str(lac)
xyplot(leite~dia|vaca, data=lac)
## modelo: de Wood (nucleo da densidade gama).
## x: representado por dia, dia após parto.
## y: representado por leite, quantidade produzida.
## th1: escala, desprovido de interpretação direta.
## th2: forma, desprovido de interpretação direta.
## th3: forma, desprovido de interpretação direta.
model <- leite~th1*dia^th2*exp(-th3*dia)
start <- list(th1=c(15,10,20), th2=c(0.2,0.05,0.5),
th3=c(0.0025,0.0010,0.0080))
rp.nls(model=model, data=lac,
start=start, subset="vaca",
assignTo="lac.fit", xlim=c(0,310))
sapply(lac.fit, coef)
##----------------------------------------------------------------------
## 5. Curvas de crescimento em placa de petri.
cre <- read.table(
"http://www.leg.ufpr.br/~walmes/data/cresmicelial.txt",
header=TRUE, sep="\t", encoding="latin1")
cre$isolado <- factor(cre$isolado)
cre$mmdia <- sqrt(cre$mmdia)
str(cre)
xyplot(mmdia~temp|isolado, data=cre)
## modelo: quadrático na forma canônica.
## x: representado por temp, temperatura da câmara de crescimento.
## y: representado por mmdia, taxa média de crescimento.
## thy: valor da função no ponto de máximo.
## thc: curvatura ou grau de especificidade à condição ótima.
## thx: ponto de máximo, temperatura de crescimento mais rápido.
model <- mmdia~thy+thc*(temp-thx)^2
start <- list(thy=c(4,0,7), thc=c(-0.05,0,-0.5), thx=c(23,18,30))
rp.nls(model=model, data=cre,
start=start, subset="isolado",
assignTo="mic.fit",
extraplot=function(thy, thx, thc){
abline(v=thx, h=thy, lty=2, col=2)
},
xlim=c(17,31), ylim=c(0,6))
t(sapply(mic.fit, coef))
##----------------------------------------------------------------------
## 6. Curva de secagem do solo em microondas.
sec <- read.table("http://www.leg.ufpr.br/~walmes/data/emr11.txt",
header=TRUE, sep="\t", encoding="latin1")
str(sec)
xyplot(umrel~tempo|nome, data=sec)
## modelo: logístico.
## x: representado por tempo, período da amostra dentro do microondas.
## y: representado por umrel, umidade relativa o conteúdo total de água.
## th1: assíntota superior.
## th2: tempo para evaporar metade do conteúdo total de água.
## th3: proporcional à taxa máxima do processo.
model <- umrel~th1/(1+exp(-(tempo-th2)/th3))
start <- list(th1=c(1,0.8,1.2), th2=c(15,0,40), th3=c(8,2,14))
rp.nls(model=model, data=sec,
start=start, subset="nome",
assignTo="sec.fit",
extraplot=function(th1, th2, th3){
abline(v=th2, h=th1/(1:2), lty=2, col=2)
})
sapply(sec.fit, coef)
lapply(sec.fit, summary)
##----------------------------------------------------------------------
## Informações da sessão.
Sys.time()
sessionInfo()
##----------------------------------------------------------------------
##----------------------------------------------------------------------
## Definições da sessão.
## Pacote para janelas interativas.
require(rpanel)
##----------------------------------------------------------------------
## Distribuição binomial.
pb <- function(panel){
with(panel,
{
x <- 0:size
px <- dbinom(x, size=size, prob=prob)
plot(px~x, type="h", ylim=c(0,max(c(px), 0.5)))
if(showEX){
abline(v=size*prob, col=2)
}
})
panel
}
panel <- rp.control(title="Binomial", size=c(300,100))
rp.slider(panel, size, from=2, to=80, initval=10, resolution=1,
action=pb, showvalue=TRUE, title="size")
rp.slider(panel, prob, from=0.01, to=0.99, initval=0.5, resolution=0.01,
action=pb, showvalue=TRUE, title="prob")
rp.checkbox(panel, showEX, action=pb, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição de Poisson.
pp <- function(panel){
with(panel,
{
x <- 0:100
px <- dpois(x, lambda=lambda)
plot(px~x, type="h", ylim=c(0,max(c(px),0.5)))
if(showEX){
abline(v=lambda, col=2)
}
})
panel
}
panel <- rp.control(title="Poisson", size=c(300,100))
rp.slider(panel, lambda, from=0.5, to=90, initval=10, resolution=0.25,
action=pp, showvalue=TRUE, title="lambda")
rp.checkbox(panel, showEX, action=pp, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição binomial negativa.
pbn <- function(panel){
with(panel,
{
x <- 0:size
px <- dnbinom(x, size=size, prob=prob)
plot(px~x, type="h", ylim=c(0,max(c(px), 0.5)))
if(showEX){
abline(v=size*(1-prob)/prob, col=2)
}
})
panel
}
panel <- rp.control(title="Binomial negativa", size=c(300,100))
rp.slider(panel, size, from=2, to=80, initval=10, resolution=1,
action=pbn, showvalue=TRUE, title="size")
rp.slider(panel, prob, from=0.01, to=0.99, initval=0.5, resolution=0.01,
action=pbn, showvalue=TRUE, title="prob")
rp.checkbox(panel, showEX, action=pbn, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição hipergeométrica.
## m the number of white balls in the urn.
## n the number of black balls in the urn.
## k the number of balls drawn from the urn.
## x the number of white balls drawn without replacement.
ph <- function(panel){
with(panel,
{
x <- max(c(0, k-m)):min(c(k, m))
## p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)
px <- dhyper(x, m=m, n=n, k=k)
plot(px~x, type="h", ylim=c(0, max(c(px), 0.5)))
if(showEX){
abline(v=k*m/(m+n), col=2)
}
})
panel
}
panel <- rp.control(title="Hipergeométrica", size=c(300,100))
rp.slider(panel, m, from=5, to=30, initval=10, resolution=1,
action=ph, showvalue=TRUE, title="Brancas")
rp.slider(panel, n, from=2, to=15, initval=5, resolution=1,
action=ph, showvalue=TRUE, title="Pretas")
rp.slider(panel, k, from=2, to=15, initval=5, resolution=1,
action=ph, showvalue=TRUE, title="Retiradas")
rp.checkbox(panel, showEX, action=ph, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição normal.
pn <- function(panel){
with(panel,
{
curve(dnorm(x, mean=mean, sd=sd), -5, 5,
ylim=c(0,1))
if(showEX){
abline(v=mean, col=2)
}
if(showVX){
d <- dnorm(mean+sd, mean, sd)
segments(mean-sd, d, mean+sd, d, col=2)
}
})
panel
}
panel <- rp.control(title="Normal", size=c(300,100))
rp.slider(panel, mean, from=-4, to=4, initval=0, resolution=0.1,
action=pn, showvalue=TRUE, title="mean")
rp.slider(panel, sd, from=0, to=3, initval=1, resolution=0.1,
action=pn, showvalue=TRUE, title="sd")
rp.checkbox(panel, showEX, action=pn, title="E(X)",
labels="Mostrar o valor esperado?")
rp.checkbox(panel, showVX, action=pn, title="sd(X)",
labels="Mostrar o desvio-padrão?")
##----------------------------------------------------------------------
## Distribuição beta.
pbt <- function(panel){
with(panel,
{
curve(dbeta(x, shape1=exp(sh1), shape2=exp(sh2)), 0, 1,
ylim=c(0,7))
if(showEX){
abline(v=exp(sh1)/(exp(sh1)+exp(sh2)), col=2)
}
})
panel
}
panel <- rp.control(title="Beta", size=c(300,100))
rp.slider(panel, sh1, from=-5, to=5, initval=0, resolution=0.1,
action=pbt, showvalue=TRUE, title="log(shape1)")
rp.slider(panel, sh2, from=-5, to=5, initval=0, resolution=0.1,
action=pbt, showvalue=TRUE, title="log(shape2)")
rp.checkbox(panel, showEX, action=pbt, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição gamma.
pg <- function(panel){
with(panel,
{
curve(dgamma(x, shape=shape, scale=scale), 0, 50,
ylim=c(0, 0.25))
if(showEX){
abline(v=shape*scale, col=2)
}
})
panel
}
panel <- rp.control(title="Gamma", size=c(300,100))
rp.slider(panel, shape, from=0.1, to=20, initval=5, resolution=0.1,
action=pg, showvalue=TRUE, title="shape")
rp.slider(panel, scale, from=0.1, to=10, initval=3, resolution=0.1,
action=pg, showvalue=TRUE, title="scale")
rp.checkbox(panel, showEX, action=pg, title="E(X)",
labels="Mostrar o valor esperado?")
##----------------------------------------------------------------------
## Distribuição Weibull.
pw <- function(panel){
with(panel,