Re: [R-br] Regressão segmentada!
Descobri o erro na função. O correto seria: segment_lpl <- function(x,b0,b1,b2,X0,X1){ (b0+b1*x)*(x<=X0)+ (b0+b1*X0)*(X0<x & x<X1)+ (b0+b1*X0-b2*(X1-x))*(x>=X1) } Só não entendo o porque que não estou conseguindo convergência, uma vez que estou fornecendo exatamente os valores que o gráfico interativo está fornecendo. Se alguém descobrir, fico grato. (S,f,P) Allaman \begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5076 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature}
Ivan,
O modelo que você passou não é contínuo nos "cotovelos". Pelos dados que
passou, penso que você queira um modelo contínuo nas trocas de taxa. Eu fiz
uns ajustes quanto à isso
# modelo não contínuo
curve(segment_lpl(x, b0=7.8, b1=-0.02, b2=-0.01, K=7.5, X0=10, X1=20), -1,
30)
# modelo tri-linear segmentado contínuo
tri.lin <- function(x, x0, b0, x1, b1, x2, b2){
## x0: intercepto
## b0: inclinação/taxa do primeiro segmento
## x1: primeiro "cotovelo", mudança de taxa
## b1: inclinação/taxa do segundo segmento
## x2: segundo "cotovelo", mudança de taxa
## b2: inclinação/taxa do terceiro segmento
(x0+b0*x)*(x<x1)+
(x0+b0*x1+b1*(x-x1))*(x>x1)*(x<x2)+
(x0+b0*x1+b1*(x2-x1)+b2*(x-x2))*(x>x2)
}
curve(tri.lin(x, x0=0, b0=1, x1=2, b1=-1, x2=4, b2=1), -1, 6)
abline(v=c(0,2,4), lty=2)
x <- c(0,0,0,10,10,10,20,20,20,35,35,35)
y <- c(7.77, 7.81, 7.81, 7.60, 7.67, 7.65, 7.60, NA, 7.66, 7.60, 7.42, 7.50)
plot(y~x)
curve(tri.lin(x, x0=7.8, b0=-0.017, x1=10, b1=0, x2=20, b2=-0.01), add=TRUE)
da <- data.frame(x=x, y=y)
da <- da[complete.cases(da),]
> n0 <- nls(y~tri.lin(x, x0, b0, x1, b1, x2, b2), data=da, trace=TRUE,
+ start=list(x0=7.8, b0=-0.017, x1=10, b1=0.001, x2=20,
b2=-0.01))
Erro em nlsModel(formula, mf, start, wts) :
matriz gradiente singular com estimativas de parâmetros iniciais
>
> tri.lin.model <- deriv3(~(x0+b0*x)*(x<x1)+
+ (x0+b0*x1+b1*(x-x1))*(x>x1)*(x<x2)+
+ (x0+b0*x1+b1*(x2-x1)+b2*(x-x2))*(x>x2),
+ c("x0","b0","x1","b1","x2","b2"),
+ function(x, x0, b0, x1, b1, x2, b2){NULL})
Erro em deriv3.formula(~(x0 + b0 * x) * (x < x1) + (x0 + b0 * x1 + b1 * :
Função '`>`' não está na tabela de derivadas
>
Porém, o problema da não convergência não é qualidade do chute mas a função
do modelo em si. Ao usar o o operador lógico "<" (e outros) a nls() não
consegue interpretar para criar a matriz gradiente. O curioso que o modelo
bi-segmentado ela faz, não sei como, mas. O mais apropriado, dado que se os
valores de x0, x1 e x2 fossem conhecidos isso seria um problema de modelo
linear, seria escrever o modelo não linear tirando vantagem disso, dessa
linearidade parcial, e usar o algorithm="plinear" na nls(). No livro MASS
tem usa seção curta, embora suficiente, para entender como fazer isso. Não
sei se os dados que você me mandou são reais ou não, mas entre ficar com
uma trilinear (6 parâmetros) passando por 4 pontos no domínio da
covariável, eu prefiro ajustar uma cúbica (pelo menos usa menos parâmetros,
é contínua). Esse trilinear é ligar com segmentos as médias, e as
inclinações são como se fossem contrastes entre essas médias
(delta_y/delta_x). Note que trilinear requer mais parâmetros do o temos em
grau de liberdade para oferecer!! Outra maneira seria usar a optim() para
fazer o ajuste.
À disposição.
Walmes.
==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
Estava vendo a discussão, e não consegui rodar isto...
# modelo não contínuo
curve(segment_lpl(x, b0=7.8, b1=-0.02, b2=-0.01, K=7.5, X0=10, X1=20), -1, 30)
Erro em eval(expr, envir, enclos) :
não foi possível encontrar a função "segment_lpl"
>
Ivan,
O modelo que você passou não é contínuo nos "cotovelos". Pelos dados que passou, penso que você queira um modelo contínuo nas trocas de taxa. Eu fiz uns ajustes quanto à isso
# modelo não contínuo
curve(segment_lpl(x, b0=7.8, b1=-0.02, b2=-0.01, K=7.5, X0=10, X1=20), -1, 30)
# modelo tri-linear segmentado contínuo
tri.lin <- function(x, x0, b0, x1, b1, x2, b2){
## x0: intercepto
## b0: inclinação/taxa do primeiro segmento
## x1: primeiro "cotovelo", mudança de taxa
## b1: inclinação/taxa do segundo segmento
## x2: segundo "cotovelo", mudança de taxa
## b2: inclinação/taxa do terceiro segmento
(x0+b0*x)*(x<x1)+
(x0+b0*x1+b1*(x-x1))*(x>x1)*(x<x2)+
(x0+b0*x1+b1*(x2-x1)+b2*(x-x2))*(x>x2)
}
curve(tri.lin(x, x0=0, b0=1, x1=2, b1=-1, x2=4, b2=1), -1, 6)
abline(v=c(0,2,4), lty=2)
x <- c(0,0,0,10,10,10,20,20,20,35,35,35)
y <- c(7.77, 7.81, 7.81, 7.60, 7.67, 7.65, 7.60, NA, 7.66, 7.60, 7.42, 7.50)
plot(y~x)
curve(tri.lin(x, x0=7.8, b0=-0.017, x1=10, b1=0, x2=20, b2=-0.01), add=TRUE)
da <- data.frame(x=x, y=y)
da <- da[complete.cases(da),]
> n0 <- nls(y~tri.lin(x, x0, b0, x1, b1, x2, b2), data=da, trace=TRUE,
+ start=list(x0=7.8, b0=-0.017, x1=10, b1=0.001, x2=20, b2=-0.01))
Erro em nlsModel(formula, mf, start, wts) :
matriz gradiente singular com estimativas de parâmetros iniciais
>
> tri.lin.model <- deriv3(~(x0+b0*x)*(x<x1)+
+ (x0+b0*x1+b1*(x-x1))*(x>x1)*(x<x2)+
+ (x0+b0*x1+b1*(x2-x1)+b2*(x-x2))*(x>x2),
+ c("x0","b0","x1","b1","x2","b2"),
+ function(x, x0, b0, x1, b1, x2, b2){NULL})
Erro em deriv3.formula(~(x0 + b0 * x) * (x < x1) + (x0 + b0 * x1 + b1 * :
Função '`>`' não está na tabela de derivadas
>
Porém, o problema da não convergência não é qualidade do chute mas a função do modelo em si. Ao usar o o operador lógico "<" (e outros) a nls() não consegue interpretar para criar a matriz gradiente. O curioso que o modelo bi-segmentado ela faz, não sei como, mas. O mais apropriado, dado que se os valores de x0, x1 e x2 fossem conhecidos isso seria um problema de modelo linear, seria escrever o modelo não linear tirando vantagem disso, dessa linearidade parcial, e usar o algorithm="plinear" na nls(). No livro MASS tem usa seção curta, embora suficiente, para entender como fazer isso. Não sei se os dados que você me mandou são reais ou não, mas entre ficar com uma trilinear (6 parâmetros) passando por 4 pontos no domínio da covariável, eu prefiro ajustar uma cúbica (pelo menos usa menos parâmetros, é contínua). Esse trilinear é ligar com segmentos as médias, e as inclinações são como se fossem contrastes entre essas médias (delta_y/delta_x). Note que trilinear requer mais parâmetros do o temos em grau de liberdade para oferecer!! Outra maneira seria usar a optim() para fazer o ajuste.
À disposição.
Walmes.
==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
Ivan conhecendo o spontos de corte o modelo se torna linear portanto a nls pode ser usada com a opçãp "plinear" e a iteracao numérica é feirta em uma dimensao menos (apenas dos pontos de conte) Isto deve ajudar muito na convergencia Em Qui, 2011-11-17 às 05:20 -0800, Ivan Bezerra Allaman escreveu:
Descobri o erro na função. O correto seria:
segment_lpl <- function(x,b0,b1,b2,X0,X1){ (b0+b1*x)*(x<=X0)+ (b0+b1*X0)*(X0<x & x<X1)+ (b0+b1*X0-b2*(X1-x))*(x>=X1) }
Só não entendo o porque que não estou conseguindo convergência, uma vez que estou fornecendo exatamente os valores que o gráfico interativo está fornecendo. Se alguém descobrir, fico grato.
(S,f,P) Allaman
\begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5076 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature} _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
participantes (4)
-
Ivan Bezerra Allaman -
Mauro Sznelwar -
Paulo J Ribeiro Jr -
Walmes Zeviani