Índice
O Curso
Material de Apoio
Área dos Alunos
Visitantes
Forum
notaR
Área Restrita
Cursos Anteriores
IBUSP
Outras Insitutições
IBUSP
Outras Insitutições
Doutorando em Biologia Vegetal, Instituto de Biologia, Unicamp.
Interesses em ecologia de comunidades vegetais e de ecossistemas.
Pensei numa função “curve finder” onde seria utilizada o AIC para que indicasse a curva que melhor se ajusta em uma relação bivariada (dois vetores). Os principais modelos de curvas a serem analisados nessas relações bivariadas seriam: no modelo linear (simples, quadrático e - talvez - polinomial), exponencial (e o modificado), o logarítmico, geométrico, (e talvez outros). A inspiração foi no curve expert, que faz bem isso, mas usa estatística F para definir a curva com melhor ajuste e possui uns graficos bem feinhos. Aí o usuário do R receberia como resposta os graficos X por Y acrescentado as linhas para cada um dos modelos e a ordem dos modelos com os melhores ajustes e suas fórmulas.
Boa proposta. Falta só definir se todos os modelos presumirão distribuição normal dos resíduos. Isto vai afetar a função que vc usará para ajuste. Se sim, será a função de regressão gaussiana não linear (nls
). Se não, terás que partir para glm's, mas aí é complicar demais. O importante é deixar claro que os modelos são Gaussianos.
fit.model=function(x,y, res.normal=TRUE) { if(res.normal==TRUE) { lm(y~x)->ls lm(y~x+I(x^2))->lq lm(y~x+I(x^2)+I(x^3))->lp nls(y~a+exp(b*x), start=c("a"=1, "b"=2)) ->expo nls(y~a+b*log(x), start=c("a"=1e-5, "b"=1e-5))->loga nls(y~a*x/(b+x), start=c("a"=1, "b"=1))->sat nls(y~a*sin(b*x), start=c("a"=1, "b"=2))->sinus nls(y~a*I(x^b), start=c("a"=1, "b"=1)) -> pm AIC(ls, lq, lp, expo, loga, sat, sinus, pm)->AT AT[,2]->AT AT[AT<0]*(-1)->Akaike layout(matrix(c(1:9), 3, 3, byrow=TRUE)) plot(x,y, main="Linha de tendencia"); panel.smooth(x, y) plot(x, y, main="Modelo Linear Simples"); abline(x, y, col="red") plot(x, y, main="Modelo Linear Quadratico"); lines(x, lq$coef[1]+lq$coef[2]*x+lq$coef[3]*x^2, col="red") plot(x, y, main="Modelo linear Polinomial"); lines(x, lp$coef[1]+lp$coef[2]*x+lp$coef[3]*x^2+lp$coef[4]*x^4, col="red") plot(x, y, main="Modelo Exponencial"); coef(expo)->coexp; lines(coexp[1]+exp(coexp[2]*x), col="red") plot(y~x, main="Modelo Logarítmico"); coef(loga)->cologa; lines(cologa[1]+cologa[2]*log(x), col="red") plot(y~x, main="Modelo de Saturação"); coef(sat)->cosat; lines(cosat[1]*x/(cosat[2]+x), col="red")->l1 plot(x,y, main="Modelo Sinusoidal"); coef(sinus)->cosin; lines(cosin[1]*sin(cosin[2]*x), col="red") plot(x,y, main="Modelo de Potência"); coef(pm)->copm; lines(copm[1]*x^copm[2], col="red") x11() AIC(ls, lq, lp, expo, loga, sat, sinus, pm)->AT AT[,2]->AT AT[AT<0]*(-1)->Akaike layout(matrix(c(1:9), 3, 3, byrow=TRUE)) plot(y~x, main="Dispersão dos pontos") qqnorm(residuals(ls), main="Modelo Linear simples", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo Linear Quadrático", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo Linear Polinomial", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo Exponencial", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Normalidade dos Resíduos da Regressão \n no modelo Logarítmico", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo de Saturação", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo Sinusoidal", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) qqnorm(residuals(ls), main="Modelo de Potência", xlab="Quantis Teóricos", ylab="Quantis amostrados"); qqline(residuals(ls)) c("a+b*x", "a+b*x+c*x^2", "a+b*x+c*x^2+d*x^3", "a+exp(b*x)", "a+b*log(x)", "a*x/(b+x)", "a*sin(b*x)", "a*(x^b)")->formulas Modelo=c("Linear Simples", "Linear Quadrático", "Linear Polinomial", "Exponencial", "Logarítmico", "de Saturação", "Sinusoidal", "de Potência") as.vector(ls$coef)-> lsc; as.vector(lq$coef)-> lqc; as.vector(lq$coef)->lpc; as.vector(coef(expo))->coexp; as.vector(coef(loga))->cologa; as.vector(coef(sat))->cosat; as.vector(coef(sinus))->cosin; as.vector(coef(pm))->copm a=c(lsc[1], lqc[1], lpc[1], coexp[1], cologa[1], cosat[1], cosin[1], copm[1]) b=c(lsc[2], lqc[2], lpc[2], coexp[2], cologa[2], cosat[2], cosin[2], copm[2]) c=c("NA", lqc[3], lpc[3], "NA", "NA", "NA", "NA", "NA") d=c("NA", "NA", lpc[4], "NA", "NA", "NA", "NA", "NA") AIC(ls, lq, lp, expo, loga, sat, sinus, pm)->AT AkaikeIndexCriterion<-AT[,2] Ordem=order(AkaikeIndexCriterion) resultado=data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) print("Alerta: Verifique se a relação entre os quantis teóricos e os quantis amostrados dos resíduos dos modelos seguem a normalidade. \n Caso os resíduos dos modelos não sigam a normalidade: Use 'res.normal=FALSE' na função.") return(resultado) } if(res.normal==FALSE) { formulas=c("y~a+b*x","y~a+b*x+c*x^2", "y~a+b*x+c*x^2+d*x^3") glm(y~x)->simp glm(y~x+I(x^2))->quad glm(y~x+I(x^2)+I(x^3))->poli AIC(simp, quad, poli)->ATG ATG[,2]->AkaikeIndexCriterion Ordem=order(AkaikeIndexCriterion) a=c(simp$coef[1], quad$coef[1], poli$coef[1]) b=c(simp$coef[2], quad$coef[2], poli$coef[2]) c=c("NA", quad$coef[3], poli$coef[3]) d=c("NA", "NA", poli$coef[4]) Gaussiano<-data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) Y=y/max(y) X=x/max(x) glm(Y~X, family='binomial')->binos glm(Y~X+I(X^2), family=binomial())->binoq glm(Y~X+I(X^2)+I(X^3), family=binomial()) ->binop AIC(binos, binoq, binop)->ATG ATG[,2]->AkaikeIndexCriterion Ordem<-order(AkaikeIndexCriterion) a=c(binos$coef[1], binoq$coef[1], binop$coef[1]) b=c(binos$coef[2], binoq$coef[2], binop$coef[2]) c=c("NA", binoq$coef[3], binop$coef[3]) d=c("NA", "NA", binop$coef[4]) Binomial<-data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) glm(y~x, family=gamma()) -> gamms glm(y~x+I(x^2)+I(x^3), family=inverse.gaussian()) ->gammq glm(y~x+I(x^2)+I(x^3), family=inverse.gaussian()) ->gammp AIC(gamms, gammq, gammp)->ATG ATG[,2]->AkaikeIndexCriterion Ordem<-order(AkaikeIndexCriterion) a=c(gamms$coef[1], gammq$coef[1], gammp$coef[1]) b=c(gamms$coef[2], gammq$coef[2], gammp$coef[2]) c=c("NA", gammq$coef[3], gammp$coef[3]) d=c("NA", "NA", gammp$coef[4]) Gamma<-data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) glm(y~x, family=inverse.gaussian()) ->igauss glm(y~x+I(x^2)+I(x^3), family=inverse.gaussian()) ->igausq glm(y~x+I(x^2)+I(x^3), family=inverse.gaussian()) ->igausp AIC(iguass, igausq, igausp)->ATG ATG[,2]->AkaikeIndexCriterion Ordem<-order(AkaikeIndexCriterion) a=c(iguass$coef[1], iguasq$coef[1], iguasp$coef[1]) b=c(iguass$coef[2], iguasq$coef[2], iguasp$coef[2]) c=c("NA", iguasq$coef[3], iguasp$coef[3]) d=c("NA", "NA", iguasp$coef[4]) GaussianoInverso<-data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) glm(y~x, family=poisson())->poiss glm(y~x+I(x^2), family=poisson()) -> poisq glm(y~x+I(x^2)+I(x^3), family=poisson()) ->poisp AIC(poiss, poisq, poisp)->ATG ATG[,2]->AkaikeIndexCriterion Ordem<-order(AkaikeIndexCriterion) a=c(poiss$coef[1], poisq$coef[1], poisp$coef[1]) b=c(poiss$coef[2], poisq$coef[2], poisp$coef[2]) c=c("NA", poisq$coef[3], poisp$coef[3]) d=c("NA", "NA", poisp$coef[4]) Poisson<-data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) glm(y~x, family=quasi())->quasis glm(y~x+I(x^2), family=quasi()) -> quasiq glm(y~x+I(x^2)+I(x^3), family=quasi()) ->quasip AIC(quasis, quasiq, quasip)->ATG ATG[,2]->AkaikeIndexCriterion Ordem<-order(AkaikeIndexCriterion) a=c(quasis$coef[1], quasiq$coef[1], quasip$coef[1]) b=c(quasis$coef[2], quasiq$coef[2], quasip$coef[2]) c=c("NA", quasiq$coef[3], quasip$coef[3]) d=c("NA", "NA", quasip$coef[4]) QuasiGaussiano=data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) Y=y/max(y) X=x/max(x) glm(Y~X, family=quasibinomial())->quabis glm(Y~X+I(X^2), family=quasibinomial()) -> quabiq glm(Y~X+I(X^2)+I(X^3), family=quasibinomial()) ->quabip AIC(quabis, quabiq, quabip)->ATG ATG[,2]->AkaikeIndexCriterion Ordem=order(AkaikeIndexCriterion) a=c(quabis$coef[1], quabiq$coef[1], quabip$coef[1]) b=c(quabis$coef[2], quabiq$coef[2], quabip$coef[2]) c=c("NA", quabiq$coef[3], quabip$coef[3]) d=c("NA", "NA", quabip$coef[4]) QuasiBinomial=data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) glm(y~x, family=quasipoisson())->quapos glm(y~x+I(x^2), family=quasipoisson()) -> quapoq glm(y~x+I(x^2)+I(x^3), family=quasipoisson()) ->quapop AIC(quapos, quapoq, quapop)->ATG ATG[,2]->AkaikeIndexCriterion Ordem=order(AkaikeIndexCriterion) a=c(quapos$coef[1], quapoq$coef[1], quapop$coef[1]) b=c(quapos$coef[2], quapoq$coef[2], quapop$coef[2]) c=c("NA", quapoq$coef[3], quapop$coef[3]) d=c("NA", "NA", quapop$coef[4]) QuasiPoisson=data.frame(Ordem, Modelo, formulas, a, b, c, d, AkaikeIndexCriterion) resultglm=array(Gaussiano, Binomial, Gamma, GaussianoInverso, Poisson, QuasiGaussiano, QuasiBinomial, QuasiPoisson) return(resultglm) print("\n\n Os valores dos vetores deve estar entre 0 e 1 para que ocorra o modelo generalizado binomial \n\n Cuidado: Para o modelo generalizado Gamma, ao menos 1 argumento deve passar por 'gamma'.") } }
fit.model package:nenhum R Documentation Ajuste de modelos entre dois vetores numéricos, usando o Índice de Critério de Akaike. Description Ajusta modelos lineares e não-lineares na relação entre dois vetores, classificando o melhor ajuste e no caso de distribuição normal dos resíduos no modelo (res.norm=TRUE) Usage: fit.model(x, y, res.normal = TRUE) Arguments: x um objeto da classe "vector", sendo com valores numéricos. Variável preditora y: um objeto da classe "vector", sendo com valores numéricos. Varipavel resposta res.normal: para resíduos com distribuição não normal use res.normal=FALSE para resíduos com distribuição gaussiana use o argumento res.normal=TRUE(default) Details: A função promove o ajuste do melhor modelo entre dois vetores, usando o Índice de Critério de Akaike (Akaike Index Criterion - AIC). Para a distribuição normal dos resíduos (res.normal=TRUE), também são dados os gráficosda variável resposta vs a variável preditora com o as curvas dos modelos. Além disso, é fornecido a distribuição teórica gaussiana vs observada dos quantis dos resíduos dos modelos. Dessa forma, o usuário pode ver se os resíduos do melhor modelo ajustado possui distribuição gaussiana. Para modelos generalizados (res.normal=FALSE), são ajustados três tipos de modelos: simples, quadrático e polinomial. Nesse caso, a resposta fornecida será um "array" com os dados de AIC e as fórmulas para cada um dos modelos generalizados testados (Gaussiano, Binomial, Gamma, GaussianoInverso, Poisson, QuasiGaussiano, QuasiBinomial, QuasiPoisson). Esses modelos são chamados de "family" pela função "glm". Warning: Os vetores não podem possuir NA's ou NAN's. Author(s): Versão original por Juliano van Melis jvmelis@yahoo.com.br References: Kelley, C. T.(1999). "Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics", no 18, ISBN: 0-89871-433-8. Nelder, J. & Wedderburn, R. (1972). "Generalized Linear Models". Journal of the Royal Statistical Society. Series A (General) 135 (3): 370–384. See Also: glm nls Examples:
Envolveria estatística circular. A minha idéia seria comparar a riqueza ou abundância de determinada(s) espécie(s) em diferentes ecounidades florestais (Oldeman 1990: clareira-construção-bioestasia-degradação). Transformaria cada ecounidade em valores angulares e os valores a serem comparados nos eixos que saem do centro do gráfico. O teste V (Jammalamadaka & SenGupta 2001) seria utilizado para verificar a preferência da(s) espécies e/ou da sinúsia vegetal por determinado ângulo (ecounidade). Como o número de parcelas em cada ecounidade não é uniforme (p.ex: eco-unidades maduras 28-37%, Cassola 2008), teria que ser utilizado valores proporcionais. Talvez fazer uma rarefação, reduzindo o número de parcelas para aquela que tiver menor número.
Novamente, a inspiração vem de outro software (Oriana), que eu já tive que usar e que tem saídas gráficas sofríveis…Além de não ser um freeware.
Cassola, H. (2008)Aspectos da estrutura fitossociológica e silvigenética em fragmentos de floresta estacional semidecídua com diferentes histórias de perturbação em Botucatu, SP. Dissertação Mestrado - ESALQ/USP, Piracicaba
Jammalamadaka, S.R., SenGupta, A. (2001)Topics in Circular Statistics. Series on multivariate analysis, vol. 5. World Scientific Publishing Co.. Singapore
Oldeman, R.A.A. (1990) Forests: Elements of Silvology. Springer-Verlag. Berlin
Também muito boa. A minha crítica não é à função em si, mas á subsjetivdiade de converter cada econunidade em um ângulo. Se vc tivesse tempo de regeneração a coisa ficaria mais precisa, concordas?
Mas isto não invalida o exercício de construir a função, que é adequado.
Usei no relatório FAPESP do ano passado essa “subjetividade” pois o meu problema era exatamente como descrever a floresta amostrada. Eu tinha classificado essas ecounidades por Oldeman (1990), e como em campo eu não achava que cada parcela de 100m² era considerada como uma única ecounidade s.s., as vezes eu achava que uma parcela estava entre clareira e construção. Por isso não teríamos somente quatro angulos, mas 8 ângulos.
A “subjetividade” ficou na análise visual em campo, pois não tinha experiência nessa classificação (Oldeman 1990). Para refinar a caracterização das ecounidades terei que usar análise multivariada para ver se as ecounidades formam grupos usando algumas variáveis características (altura do dossel médio, n de árvores mortas, IIC médio, n de árvores do presente/passado/futuro (Oldeman 1990), etc…). Não vou colocar o n de lianas pq talvez enviese a relação que eu quero (ecounidades vs sinúsia de trepadeiras).
Concordo que se eu tivesse o tempo de regenaração (que talvez eu consiga usando os dados do Roque (2001 e 2007), o uso da estatística circular seria mais robusta.
Se as outras duas propostas forem muito complicadas, também seria útil para um dos meus projetos a análise de covariância (ANCOVA), onde as retas de dois modelos lineares seriam comparadas, verificando se as retas são iguais (mesma inclinação) e paralelas (diferentes interceptos), ou se são distintas (valores de intercepto e inclinação significativamente diferentes).
Este é tão fácil de fazer diretamente no R que não vejo muito sentido em criar uma função para isto. Mas as duas propostas anteriores são viáveis.
Viro hippie.
Não é original, mas é sempre uma opção.
Não é original, mas seria desafiadora. O único artesanato que sou capaz de fazer são aviões de papel, no formato "estrela da morte" (by Star Wars).