====== Lígia Haselmann Apostólico ======
\\
{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:10400875_10202394260078202_7026659431929350814_n.jpg?150|}}
Mestranda em Zoologia, Instituto de Biociências, USP.
O título da minha dissertação é: "Caracterização do dimorfismo intrassexual masculino de Doryteuthis plei (Mollusca: Cephalopoda), como base para a compreensão dos mecanismos de competição espermática em lulas", orientado pelo prof. José Eduardo A.R. Marian.
\\
===== Meus Exercícios =====
\\
Link para meus exercícios resolvidos - [[.:exec]]
\\
===== Propostas - Trabalho Final =====
\\
Link para minhas propostas de trabalho final - [[.:trabfinal]]
\\
\\
===== Trabalho Final =====
\\
==== Help da Função ====
isd.analysis package:base R Documentation
Função para testar a existência de dimorfismo intrassexual a partir de medidas morfométricas.
Description:
Investigação inicial sobre a possível existência de dimorfismo intrassexual dentro de uma mesma população ou espécie.
A função baseia-se na aplicação de três modelos. O Modelo 1 testa se a relação entre as variáveis morfométricas de interesse é linear. Caso isto ocorra, o modelo conclui que não há dimorfismo intressexual e a investigação é interrompida. No entanto, caso a relação não seja linear, serão testados outros dois modelos para averiguar a existência de um switchpoint, ou seja, um valor específico a partir do qual a relação torna-se não-linear. O Modelo 2 testa se a relação entre as variáveis torna-se descontínua a partir do switchpoint, enquanto que o Modelo 3 testa se a relação entre as variáveis é alterada, mas sem descontinuidade.
A função retorna ao usuário o sumário dos modelos aplicados e o valor de switchpoint, caso este exista, que melhor se adequa aos modelos. A função retorna também, em sua janela gráfica, gráficos de dispersão e os gráficos padrão resultantes dos modelos lineares.
Usage:
isd.analysis (dados, switchpoint = 100, eixo.x = "Variável x", eixo.y = "Variável y", plot = "todos")
Arguments:
dados Conjunto de dados. Deve ser matriz ou data.frame com duas colunas. Cada linha deve representar um indivíduo. A primeira coluna deve representar a variável preditora (e.g., tamanho corporal dos espécimes); a segunda coluna deve representar a variável resposta (característica de interesse, e.g., peso gonadal, tamanho de espinhos, tamanho de armamentos).
switchpoint Valor numérico. Número de switchpoints que serão simulados no modelo. No modo default, equivale a 100.
eixo.x Título para o eixo x. Nome que será inserido no argumento xlab da função plot. No modo default, equivale a "Variável x".
eixo.y Título para o eixo y. Nome que será inserido no argumento xlab da função plot. No modo default, equivale a "Variável y".
plot Tipos de gráficos gerados pela função. Para plot="todos" (default), serão plotados gráfico de dispersão da variável y em função da variável x e os quatro gráficos padrão gerados em modelos lineares pela inserção do comando plot(lm()).
... Argumentos opicionais para plot: ver seção "Note".
Details:
Nos processos de seleção sexual, a forte competição entre machos pode resultar em dimorfismo intrassexual, caracterizado pela descontinuidade de traços morfológicos, fisiológicos e de ciclo de vida entre indivíduos do mesmo sexo. O objetivo da função isd.analysis é fazer uma investigação inicial sobre a possível existência de dimorfismo intrassexual dentro da população ou espécie de interesse, partindo-se da premissa de que uma das formas de se detectar este tipo de dimorfismo ocorre através da identificação de descontinuidade nos traços morfológicos.
A análise realizada pela função é composta pela aplicação de três modelos lineares.
A primeira parte da análise tem como objetivo testar se a relação entre as variáveis morfométricas de interesse apresenta-se de forma linear. Para isso, testa-se um modelo linear (Modelo 1), cuja equação é representada por: ln(Y)= α0 + α1*ln(X) + α2*ln(X)^2 + ε, na qual ln=log natural, α=coeficientes de regressão e ε=erro associado. Caso a relação entre as variáveis seja linear, o modelo conclui que não há dimorfismo intressexual e a investigação é interrompida, sem a aplicação dos demais modelos.
Caso o coeficiente α2 seja significantemente diferente de zero (i.e., a relação entre as variáveis não seja linear), serão testados dois novos modelos para averiguar a existência de um switchpoint, ou seja, um valor a partir do qual a relação entre as variáveis torna-se não-linear.
No Modelo 2, a equação utilizada é: Y= β0 + β1*X + β2*(X-X0)*D + β3*D + ε, na qual β=coeficientes de regressão, X0=switchpoint testado, ε=erro associado e D=constante condicional. Neste modelo, testa-se se, a partir do switchpoint (X0) a relação entre as variáveis torna-se descontínua (ver figura 1a para detalhes). Caso o coeficiente β3 seja significantemente diferente de zero, conclui-se que o dimorfismo ocorre – de forma descontínua – a partir do switchpoint selecionado. Neste caso, a investigação é interrompida, sem a aplicação do último modelo.
Caso o coeficiente β3 não seja significantemente diferente de zero, será aplicado o Modelo 3, cuja equação é: Y= β0 + β1*X + β2*(X-X0)*D + ε. Neste modelo, testa-se se, a partir do switchpoint (X0), a relação linear entre as variáveis é alterada (evidenciada por alteração na inclinação da reta), mas sem descontinuidade (ver figura 1b para detalhes).
Em ambos os Modelos 2 e 3, os valores de X0 testados dependem do número inserido pelo usuário no argumento switchpoint. A partir da amplitude de valores da variável x e do número de switchpoints que o usuário deseja testar, serão calculados os intervalos que a variável X0 irá assumir. Por exemplo, se os valores de x variam entre 50mm (mínimo) e 150mm (máximo), e o usuário deseja testar intervalos de 5 a 5mm, deve ser inserido switchpoint = 20.
Para determinar qual valor de switchpoint mais se adequa ao modelo, valores de X0 são simulados na equação do Modelo 2 e seleciona-se aquele que apresenta maior valor de R^2 ajustado.
Para mais detalhes a respeito do modelo proposto por Eberhard & Gutiérrez (1991), ver referências e arquivo anexo.
Value:
A função isd.analysis retorna ao usuário, na forma de mensagem no console, um resumo da função, contendo quais modelos foram aplicados e as principais conclusões baseadas nos valores de p associados aos coeficientes dos modelos.
Caso os Modelos 2 e 3 sejam aplicados, serão retornados no console também o valor do switchpoint selecionado e o valor do R^2 ajustado para o melhor switchpoint.
A função também retorna no console o summary dos modelos aplicados, na forma de lista.
Se somente o Modelo 1 for aplicado, a lista retorna como:
[[1]] : summary do Modelo 1
[[2]] : summary do Modelo 0 (ln(Y)= α0 + α1*ln(X) + ε)
Se os Modelos 1 e 2 forem aplicados, a lista retorna como:
[[1]] : summary do Modelo 1
[[2]] : summary do Modelo 2
Se os Modelos 1, 2 e 3 forem aplicados, a lista retorna como:
[[1]] : summary do Modelo 1
[[2]] : summary do Modelo 2
[[3]] : summary do Modelo 3
Na janela gráfica, a função retorna 5 gráficos: (I) um gráfico de dispersão de x e y com retas ajustadas de acordo com o melhor modelo que descrever o conjunto de dados e (II) os quatro gráficos padrão gerados a partir de modelos lineares.
Para saber quais gráficos são retornados ao usuário, dependendo do valor inserido no argumento plot, ver seção "Note".
Warning:
A função é interrompida e mensagens de erro são retornadas ao usuário caso seja inserido, no argumento dados: (1) algum objeto que seja diferente de uma matriz ou um dataframe ou (2) a matriz ou data.frame não contenha exatamente duas colunas.
A função é executada, mas retorna mensagem de atenção, caso o conjunto de dados (inserido no argumento dados) apresente valores faltantes, valores iguais a zero ou NAs. Neste caso, a função elimina as linhas que contenham estes dados e é executada normalmente.
Note:
Detalhes sobre argumento plot:
plot == "todos" Saída gráfica com cinco gráficos: gráfico de dispersão de x e y com retas ajustadas de acordo com o melhor modelo que descrever o conjunto de dados e quatro gráficos padrão gerados a partir de modelos lineares
plot == "resultado" Saída gráfica com apenas o gráfico de dispersão de x e y com retas ajustadas de acordo com o melhor modelo que descrever o conjunto de dados
plot == "modelo" Saída gráfica com quatro gráficos padrão gerados a partir de modelos lineares (gráfico de dispersão resíduo versus valor ajustado, gráfico quantil-quantil normal dos resíduos, gráfico de dispersão da raiz quadrada do valor absoluto do resíduo padronizado versus valor ajustado e gráfico de dispersão do resíduo padronizado versus leverage, com a distância de Cook)
plot == "nenhum" Nenhum gráfico será gerado na função
Author(s):
Lígia Haselmann Apostólico
ligia.haselmann@ib.usp.br
References:
Eberhard, W.G. & Gutiérrez, E.E. Male dimorphisms in beetles and earwings and the question of developmental constrains. Evolution, v.45, n.1, p.18-28, 1991.
Iwata, Y. & Sakurai, Y. Threshold dimorphism in ejaculate characteristics in the squid Loligo bleekeri. Mar Ecol Prog Ser, v.345, p.141-146, 2007.
See Also:
Arquivo "Modelos de Eberhard & Gutiérrez.pdf", anexado abaixo
Figura 1, retirada de Eberhard & Gutiérrez, 1991
Examples:
Os exemplos utilizados para testar a função são referentes a medidas reais obtidas em exemplares adultos de machos de lulas da espécie Doryteuthis plei coletadas em diferentes períodos ao longo dos anos de 1995 a 1998 no litoral norte do estado de São Paulo (SP,Brasil)
#################
### Exemplo 1 ###
#################
#Dados empíricos de comprimento do manto (mm) e peso do sistema reprodutor (g) de 494 lulas
#sistema reprodutor total equivale a gônada (testículo) e órgãos acessórios (órgão espermatofórico + saco espermatofórico)
#vetor numérico com 494 medidas de comprimento do manto (mm)
exemplo1.x <-c(240,165,148,193,98,257,134,117,269,203,114,100,123,171,273,212,335,282,182,275,232,207,200,296,296,298,235,318,340,339,236,278,248,227,268,194,217,281,213,156,119,169,103,164,174,164,205,117,207,212,194,197,223,218,273,153,126,204,94,227,208,112,173,260,235,199,224,248,155,220,113,247,131,201,174,192,152,273,239,168,218,237,287,222,256,258,231,295,154,267,211,200,106,217,188,174,211,255,83,217,118,233,238,287,243,198,169,192,243,222,217,154,195,117,134,75,253,98,144,216,121,229,271,226,186,229,154,213,198,127,227,188,228,179,215,117,191,248,203,248,210,118,210,153,208,177,144,101,188,152,219,177,142,156,59,201,216,125,143,222,230,150,214,293,256,131,164,237,194,185,224,228,204,225,249,194,191,227,239,170,237,241,139,182,217,162,223,224,182,143,235,209,127,171,267,227,122,150,168,155,134,165,192,199,164,89,179,259,274,218,164,234,250,222,244,174,273,208,132,237,187,219,257,162,187,62,161,222,295,227,272,228,177,125,203,226,176,145,128,177,141,151,127,128,232,106,229,149,130,102,182,144,117,176,121,223,184,283,178,263,254,253,192,264,153,229,288,279,278,270,192,278,163,96,266,229,226,166,236,80,111,164,106,132,225,328,208,263,222,129,228,176,191,111,233,201,170,160,189,226,201,147,223,161,206,138,111,121,198,164,177,194,158,158,144,188,198,126,67,233,150,110,222,235,156,187,158,222,114,164,131,156,118,120,256,240,186,156,119,145,148,182,119,156,100,193,211,322,244,125,112,143,158,210,178,170,321,153,260,155,124,135,211,139,174,194,235,220,121,209,138,184,151,131,152,77,128,144,70,173,163,90,115,203,103,102,85,192,196,250,220,116,111,196,236,161,114,94,171,198,118,132,185,141,231,114,150,130,143,171,139,205,181,244,238,191,229,191,166,100,225,91,125,152,133,231,190,121,196,176,150,280,82,174,250,169,130,252,211,186,230,222,112,232,105,153,251,211,222,203,121,185,143,197,148,210,184,198,233,140,133,193,92,169,212,235,215,173,185,149,126,111,115,246,106,163,206,136,132,120,126,146,176,241,219,203,149,166,196,139,153,145,187,236)
#vetor numérico com 494 medidas de peso do sistema reprodutor total (g)
exemplo1.y <-c(1.93,1.31,0.63,0.83,0.33,2.56,0.84,0.54,1.34,1.14,0.55,0.37,0.54,0.82,3.04,0.89,3.64,2.75,0.35,1.73,1.98,0.98,0.70,2.15,0.92,2.13,1.66,1.85,2.82,0.76,1.47,1.60,2.32,1.30,1.78,1.50,1.37,2.18,2.12,0.97,0.53,1.26,0.51,0.95,1.39,1.29,2.00,0.48,0.74,2.09,0.91,2.47,1.33,1.74,1.96,0.86,0.73,1.02,0.73,0.75,0.70,0.53,0.72,1.14,0.79,1.09,0.69,1.94,1.13,0.36,0.55,1.76,0.64,1.46,0.42,1.28,0.51,2.65,2.32,1.18,1.49,2.33,0.93,1.38,1.62,0.97,1.49,0.93,0.66,3.35,2.01,1.77,0.45,3.03,1.59,1.12,0.88,2.55,0.26,2.30,0.63,1.60,2.04,3.24,2.03,1.74,1.41,1.09,1.97,1.87,1.78,0.93,1.42,0.45,0.88,0.23,1.73,0.24,0.53,1.63,0.24,2.04,1.82,2.23,1.62,1.38,0.84,2.30,1.21,1.76,2.63,1.98,2.38,1.03,2.41,1.58,1.48,2.78,1.72,2.40,1.12,0.27,2.40,1.45,1.95,1.51,0.90,0.59,1.04,1.19,1.88,0.94,0.53,0.92,0.51,2.01,1.52,0.45,0.65,1.14,1.59,1.21,2.26,3.07,1.28,0.57,1.51,1.68,1.26,1.68,2.01,1.30,2.21,1.23,1.69,1.31,1.88,1.62,1.59,1.59,1.44,2.86,0.34,1.06,1.37,0.88,1.99,1.43,1.18,0.66,1.86,1.05,0.61,1.99,3.13,2.65,0.45,1.06,1.76,1.50,1.03,0.99,1.57,1.65,0.64,0.40,1.25,2.02,1.75,2.37,1.60,2.70,1.43,1.33,1.13,1.49,2.51,2.64,1.41,1.24,2.30,2.86,2.55,1.11,2.08,2.95,1.36,2.58,2.08,2.22,2.58,2.81,1.45,2.46,1.80,1.21,1.31,0.94,0.60,1.07,0.51,0.65,0.39,0.26,0.47,0.14,0.63,0.62,0.46,0.41,1.07,0.49,0.22,0.92,0.47,1.49,1.47,0.87,0.46,1.03,1.34,2.40,0.70,2.12,0.54,0.92,1.13,0.75,1.36,0.83,1.10,1.08,0.35,0.33,1.56,2.35,1.89,1.40,1.19,0.17,0.43,0.93,0.32,0.36,1.61,3.69,1.34,3.66,0.99,0.61,1.71,0.95,0.69,0.40,2.36,0.60,0.67,1.28,0.66,0.74,1.29,0.61,1.25,0.61,0.87,0.59,0.37,0.35,2.09,0.63,0.64,0.80,0.37,0.57,0.41,0.96,0.84,0.39,0.37,0.82,0.73,0.33,2.33,0.69,0.85,0.37,1.28,0.81,0.42,0.82,0.54,0.29,0.36,0.15,1.10,2.10,1.18,0.40,0.42,0.59,0.55,0.98,0.74,0.64,0.17,0.89,1.14,1.74,2.11,0.65,0.28,0.73,0.75,0.89,0.76,0.54,2.28,0.44,1.21,0.49,0.50,0.45,0.89,0.78,1.30,1.35,2.08,2.31,0.93,1.13,0.61,1.47,1.02,0.69,0.49,0.43,0.67,0.81,0.30,0.94,1.42,0.24,0.58,0.87,0.46,0.46,0.18,0.95,1.13,2.05,1.05,0.40,0.41,1.50,2.45,0.62,0.44,0.13,0.62,1.15,0.27,0.43,1.48,0.61,1.44,0.86,0.84,0.66,0.95,0.65,1.13,0.89,1.09,1.66,0.69,1.17,1.87,1.23,1.35,0.58,2.26,0.40,0.42,0.29,0.52,0.84,1.26,0.37,1.04,1.00,0.46,0.87,0.50,0.99,1.25,0.90,0.60,3.46,2.69,1.49,2.71,2.71,0.70,3.33,0.44,1.25,2.48,1.20,1.02,2.27,0.98,0.68,0.78,0.72,1.12,1.73,2.00,1.40,2.33,0.78,1.17,1.63,0.64,1.51,1.86,2.37,1.75,1.47,1.00,1.03,1.30,1.12,0.88,3.72,0.60,1.74,1.84,1.32,0.85,0.48,0.77,0.79,1.53,1.63,1.12,0.84,1.23,1.23,1.89,0.81,0.69,0.89,1.46,1.55)
exemplo1 <- data.frame (exemplo1.x,exemplo1.y) #criação de um data.frame com duas colunas para análise dos dados
isd.analysis (dados=exemplo1, switchpoint=100, eixo.x="Comprimento do manto (mm)", eixo.y="Peso do sistema reprodutor (g)", plot="todos") #aplicação da função
#################
### Exemplo 2 ###
#################
#Dados empíricos de comprimento do manto (mm) e peso úmido total (g) de 293 lulas
#vetor numérico com 293 medidas de comprimento do manto (mm)
exemplo2.x <- c(172,206,257,258,260,235,190,246,156,176,182,224,213,285,200,188,216,231,168,201,267,222,214,240,165,148,193,98,257,134,117,269,203,114,100,123,171,273,212,335,323,282,182,199,275,232,207,200,296,296,298,235,318,340,339,117,134,222,75,253,98,144,216,121,127,128,232,106,229,149,130,102,182,144,117,176,121,223,184,283,178,263,254,253,192,264,153,229,288,279,278,270,192,278,163,253,266,229,226,166,236,80,111,164,106,132,225,328,208,263,222,129,228,176,247,191,111,213,233,201,170,218,160,189,226,201,147,200,223,161,206,138,111,121,198,164,177,194,158,158,144,188,198,126,67,233,150,110,222,235,156,187,158,222,114,164,131,156,118,120,256,240,186,156,119,145,148,182,142,119,156,100,193,211,220,243,302,270,226,232,195,207,141,219,199,132,308,189,136,186,168,148,203,150,272,151,135,350,260,293,280,322,316,221,288,203,265,116,241,174,181,193,209,214,224,188,208,197,135,201,213,204,193,210,187,176,204,224,151,265,286,268,295,205,255,283,276,279,335,228,281,267,281,271,287,250,254,274,258,273,137,253,162,287,249,328,287,282,205,279,266,238,260,242,298,340,187,229,151,215,271,293,271,301,136,316,252,255,253,231,303,222,294,249,272,218,259,269,236,261,294,305,295)
#vetor numérico com 293 medidas de peso úmido total (g)
exemplo2.y <- c(83.22,137.83,198.50,193.74,183.60,179.96,113.55,179.22,63.05,100.10,89.70,128.60,121.60,176.90,104.70,114.90,143.20,161.10,72.20,113.10,200.30,164.60,147.40,103.40,51.80,50.60,65.40,19.10,135.70,38.00,33.00,112.90,70.40,31.80,22.70,33.10,46.70,167.10,140.90,308.00,220.40,258.70,91.20,84.00,202.20,160.10,113.60,95.10,213.50,228.40,237.30,144.10,161.20,225.40,263.40,37.50,42.02,144.30,14.00,133.90,17.00,39.00,92.90,25.50,31.50,37.10,95.50,23.00,93.30,38.80,37.60,26.60,78.90,37.10,30.00,70.60,29.40,93.60,53.70,151.50,68.00,137.20,125.00,121.70,63.70,132.70,57.30,78.80,133.30,140.60,119.10,108.00,62.50,104.10,50.50,140.10,134.50,116.60,92.20,65.80,110.40,10.30,19.30,72.10,23.80,51.40,101.20,281.90,110.50,189.60,114.30,41.90,135.70,66.50,123.80,78.90,25.80,86.60,112.20,71.40,53.30,74.10,50.50,58.16,83.10,59.80,37.30,65.40,89.40,45.60,70.70,44.70,29.60,36.80,108.70,61.80,68.30,73.80,65.00,44.60,50.60,89.00,63.40,38.30,11.70,86.00,57.00,34.40,135.40,130.30,66.60,104.90,61.30,88.60,32.90,66.60,35.80,42.00,36.90,34.20,163.60,168.50,107.80,51.80,36.50,62.60,48.70,98.60,64.20,46.50,73.00,24.00,90.90,130.00,155.53,168.96,223.55,209.42,109.50,159.40,104.90,114.90,58.30,141.20,114.60,76.40,238.80,120.20,45.80,97.40,86.90,52.10,72.60,59.50,179.50,65.30,52.00,210.50,165.70,240.40,181.50,246.70,252.70,141.80,198.50,98.00,191.10,38.20,177.70,75.40,96.60,86.20,117.40,135.00,157.80,102.50,108.60,129.20,47.20,102.40,131.20,109.00,93.00,113.40,107.50,86.40,128.40,111.20,66.80,174.50,149.20,122.50,220.70,104.70,182.60,156.20,210.30,185.10,262.90,173.00,202.20,157.60,183.00,176.60,204.20,183.20,167.40,176.60,170.20,155.50,42.90,125.30,77.00,187.50,156.50,247.10,221.60,154.00,109.30,132.60,149.10,148.10,160.00,127.80,192.70,226.80,72.00,118.90,53.40,90.30,164.40,171.90,174.00,203.00,41.40,242.70,221.70,113.70,187.60,147.90,187.90,142.10,202.20,168.50,171.00,108.00,163.00,148.00,123.00,165.00,159.00,242.00,193.00)
exemplo2 <- data.frame (exemplo2.x,exemplo2.y) #criação de um data.frame com duas colunas para análise dos dados
isd.analysis(dados=exemplo2, switchpoint=100, eixo.x="Comprimento do manto (mm)", eixo.y="Peso úmido total (g)", plot="resultado") #aplicação da função
#################
### Exemplo 3 ###
#################
#Dados empíricos de comprimento do manto (mm) e peso úmido total (g) de 486 lulas
#vetor numérico com 486 medidas de comprimento do manto (mm)
exemplo3.x <- c(130,186,125,192,236,278,248,227,268,194,217,281,213,156,119,169,103,164,174,164,205,117,207,212,194,197,223,218,273,153,126,204,94,227,208,112,173,260,235,199,224,248,155,220,113,247,131,201,174,192,152,273,239,168,218,237,287,222,256,258,231,295,154,267,211,200,106,217,188,174,211,255,83,217,118,233,238,287,243,288,198,169,192,243,222,217,154,195,229,271,226,186,229,154,213,198,127,227,188,228,179,215,117,191,248,203,248,210,118,315,210,153,208,177,144,101,188,152,219,177,142,156,59,201,216,125,143,222,230,150,214,293,256,131,164,237,194,185,224,228,204,225,249,194,223,191,227,239,228,170,237,241,139,182,217,162,223,224,182,143,235,209,127,171,267,227,122,150,168,155,134,165,192,199,164,89,179,142,259,274,218,164,234,250,222,244,174,273,208,132,237,187,219,257,162,187,62,161,222,252,295,227,272,228,177,125,203,226,176,145,128,177,141,151,120,133,113,126,158,155,96,271,249,264,216,271,245,240,263,231,295,144,70,173,163,111,282,294,283,276,265,237,231,256,282,266,270,161,266,235,192,324,250,200,236,249,296,222,253,119,174,183,254,218,276,137,111,164,322,182,203,136,159,178,86,213,213,157,260,233,243,120,214,203,195,210,272,269,222,212,223,161,261,290,246,286,171,222,196,195,219,237,250,270,213,260,252,230,228,131,252,224,139,192,122,186,177,162,164,169,120,212,230,119,132,228,97,194,250,158,247,215,235,269,232,226,186,207,263,291,159,240,207,178,284,123,238,118,250,234,201,204,277,177,273,278,186,215,254,312,292,210,259,250,302,285,254,91,125,272,261,121,242,117,205,106,257,176,234,137,145,174,130,268,193,139,266,120,161,108,131,104,101,87,242,219,185,227,208,122,93,248,257,266,248,268,187,162,95,217,216,194,252,298,186,297,249,246,263,234,174,215,242,215,250,143,166,232,264,252,276,280,177,106,232,202,215,227,133,232,213,295,262,211,250,220,183,195,252,313,226,285,217,257,208,220,176,210,299,238,271,221,176,241,260,240,225,190,203,202,197,216,255,224,193,170,203,233,147,285,212,241,204,185,226,264)
#vetor numérico com 486 medidas de peso úmido total (g)
exemplo3.y <- c(38.30,68.50,35.80,99.10,116.40,132.30,125.60,108.80,131.10,76.30,86.70,129.70,97.00,45.50,30.10,56.60,23.40,46.80,54.40,59.50,74.80,25.40,70.70,78.00,68.50,113.50,130.40,95.90,157.80,65.11,35.60,67.70,33.30,116.60,90.60,27.50,53.10,154.30,105.40,72.10,102.90,139.20,66.30,67.70,30.30,97.30,36.90,67.10,59.30,64.30,56.90,186.70,101.90,49.40,107.10,123.20,43.70,78.60,116.60,156.60,115.60,133.70,26.40,158.90,78.30,69.90,19.60,75.80,75.70,61.80,67.80,143.40,13.40,109.40,43.10,106.40,95.40,194.40,160.70,239.60,79.50,61.20,77.10,114.40,108.40,94.40,55.00,74.30,92.30,169.00,92.50,73.70,125.00,41.60,116.20,105.90,97.80,129.30,86.60,154.60,69.10,108.10,26.80,65.40,148.50,100.30,116.40,85.50,25.50,217.70,107.70,58.00,122.30,66.30,43.30,23.90,54.10,46.70,50.00,65.50,34.30,52.60,9.70,78.10,104.20,28.60,40.00,92.00,76.00,37.20,104.80,172.50,115.70,40.10,58.30,123.20,77.40,66.90,94.20,98.90,89.80,99.40,124.70,80.20,99.80,74.90,107.00,107.40,105.50,54.30,125.10,153.00,34.20,48.40,102.10,44.20,91.20,114.50,64.30,45.20,148.00,88.50,38.90,88.30,197.60,124.20,33.30,47.60,59.30,59.90,44.40,72.50,66.30,78.40,51.10,17.00,63.10,49.50,164.40,205.60,88.90,51.50,140.40,138.30,144.60,70.10,63.10,153.60,106.80,41.90,105.30,60.50,95.60,146.40,47.40,73.00,163.00,93.10,121.10,121.80,151.60,84.70,132.10,83.70,58.70,20.80,109.40,107.50,66.90,48.80,36.50,64.40,38.80,41.30,32.00,45.90,26.20,30.80,51.10,42.30,18.50,174.30,138.20,134.60,110.40,184.70,160.40,137.70,176.20,116.90,202.80,47.26,11.96,65.24,55.40,32.50,175.70,177.70,174.20,153.10,139.00,121.80,113.00,172.20,188.20,140.40,139.90,47.20,184.50,135.20,65.40,226.40,168.50,91.70,129.30,147.30,212.90,123.00,153.10,35.60,68.90,73.80,186.00,131.10,240.10,37.70,32.00,61.10,286.50,94.90,100.80,55.70,78.10,86.20,21.70,127.20,95.70,50.50,115.50,95.00,146.20,23.00,89.50,132.30,112.00,113.90,176.30,179.00,152.00,131.90,135.60,67.80,172.50,220.70,144.80,204.50,77.60,120.10,103.90,93.40,124.00,161.60,167.20,151.30,121.50,175.70,156.30,169.00,126.70,52.90,119.40,110.50,42.70,60.50,32.20,78.50,83.40,48.50,68.40,46.40,27.90,95.70,117.70,31.20,41.00,95.10,18.40,75.30,123.40,50.30,129.90,151.50,133.00,154.90,136.10,123.40,74.00,102.40,166.30,162.90,61.80,98.00,62.30,67.80,175.10,33.20,128.40,34.60,122.10,112.20,67.30,62.60,166.90,52.90,130.50,171.90,65.60,85.10,123.20,184.40,160.60,80.20,130.90,161.30,188.00,165.70,152.10,20.20,37.40,175.00,163.20,39.30,94.30,33.20,98.60,24.00,111.50,79.70,139.50,40.50,48.70,68.70,38.70,152.40,55.90,44.90,117.60,26.70,54.60,23.60,35.30,23.60,21.00,14.90,131.60,108.70,65.90,130.50,88.70,29.20,16.50,108.50,161.60,165.40,145.10,172.50,70.80,58.00,20.00,108.30,94.90,81.70,126.70,175.20,67.50,204.50,146.20,133.00,145.40,111.10,51.70,111.20,129.10,89.90,43.50,38.00,64.60,86.30,126.50,142.90,176.90,174.30,83.90,26.70,125.60,82.10,116.90,123.10,32.30,124.40,87.20,242.80,160.30,86.80,95.40,125.10,64.80,52.80,127.60,242.70,135.40,174.90,87.20,142.50,103.30,82.50,65.40,107.70,243.40,95.60,164.90,106.30,77.20,113.30,163.70,144.20,168.70,90.20,94.30,102.40,93.50,106.80,129.30,112.60,96.10,72.10,118.50,110.20,40.80,173.90,105.50,112.90,86.40,58.80,82.30,148.30)
exemplo3 <- data.frame (exemplo3.x,exemplo3.y) #criação de um data.frame com duas colunas para análise dos dados
isd.analysis(dados=exemplo3, switchpoint=100, eixo.x="Comprimento do manto (mm)", eixo.y="Peso úmido total (g)", plot="resultado") #aplicação da função
## FIM ##
==== Arquivo da Função ====
Link para arquivo da função (.r):{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:funcao_isd.analysis.r|}}
==== Código da Função ====
#################################################################################
########### Função para detecção de dimorfismo intrassexual em machos ###########
################################ (isd.analysis) #################################
#################################################################################
#
######################################
##### Descrição básica da função #####
######################################
#
# Modelo baseado na análise alométrica proposta por Eberhard & Gutiérrez (1991) - ver referência
# O modelo parte da premissa que o dimorfismo intrassexual entre machos caracteriza-se pela descontinuidade de traços morfológicos. Em outras palavras, o modelo prevê que é possivel detectar dimorfismo se a relação entre as variáveis morfométricas não for linear
#
# Essa análise é feita em trê etapas:
# (1) Verifica-se se a relação entre as variáveis é linear (Equação 1 do modelo)
# (2) Se a relação não é linear, assume-se que o dimorfismo pode se apresentar de duas formas:
# (2A) De forma descontínua, i.e., com uma interrupção na linearidade entre as variáveis (Equação 2 do modelo) - ver figura 1a para detalhes
# (2B) De forma contínua, i.e., com uma mudança na inclinação da linearidade entre as variáveis (Equação 3 do modelo) - ver figura 1b para detalhes
#
######################################
####### Função isd.analysis() ########
######################################
#
isd.analysis <- function (dados, switchpoint=100, plot="todos", eixo.x="Variável x", eixo.y="Variável y") { #início da função; para descrição dos argumentos, ver arquivo de Ajuda
#
#######################################
### conferência da entrada de dados ###
#######################################
#
if (is.matrix(dados) != TRUE & is.data.frame(dados) != TRUE) #teste lógico; se dados não forem matriz ou data.frame, parar a função
{
stop ("Erro na função:","\n","\t","Argumento 'dados' deve ser matriz ou data.frame","\n","\t","Corrija esta condição antes de continuar","\n","\t","Leia arquivo de Ajuda da função para mais detalhes") #se argumento 'dados' não for matriz ou data.frame, a função é interrompida e o usuário recebe uma mensagem de erro na tela
}
if (is.matrix(dados) == TRUE) #teste lógico para saber se os dados estão no formato de matriz
{
dados <- as.data.frame(dados) #se os dados estiverem no formato de matriz, eles serão transformados em data.frame
}
if (ncol(dados) != 2) #teste lógico para saber se o conjunto de dados tem somente 2 colunas
{
stop ("Erro na função:","\n","\t", "Argumento 'dados' deve conter apenas 2 colunas","\n","\t","Corrija esta condição antes de continuar","\n","\t","Leia arquivo de Ajuda da função para mais detalhes")
#acima: caso o argumento 'dados' não tenha exatamente 2 colunas, a função é interrompida e o usuário recebe uma mensagem de erro na tela
}
#
colnames(dados) <- c("x","y") #altera o nome das colunas do argumento 'dados';primeira coluna corresponde à variável x e segunda coluna à variável y
#
##################################
### remoção de dados faltantes ###
##################################
#
#como a análise trata de dados morfométricos, valores iguais a "zero" são considerados com dados faltantes
dados[dados == 0] <- NA #transforma valores iguais a zero em NAs; caso haja zeros no conjunto de dados, este comando irá gerar um alerta ao final da função, avisando sobre a produção de NAs
#
#se houver células com valores faltantes (com símbolos como "-" na planilha original de dados,por exemplo)
dados$x <-as.numeric(dados$x) #transforma valores em números; células com "-" serão transformados em NA
dados$y <-as.numeric(dados$y) #transforma valores em números; células com "-" serão transformados em NA
#
#se houver NAs no dataframe (presentes no conjunto de dados original, ou criados pelos comandos acima), eles serão removidos para a análise
if (sum(is.na(dados$x)) != 0 | sum(is.na(dados$y)) != 0) #teste lógico; se houver NA em alguma das colunas, o valor da soma será diferente de zero
{
dados <- na.omit(dados) #linhas que contêm NAs são omitidas do conjunto de dados
cat("\n","Aviso:","\n","\t","Há NAs ou valores iguais a zero no conjunto de dados","\n","\t","\t","Linhas contendo NAs e/ou zero foram removidas para a análise","\n","\n")
#acima: usuário receberá mensagem de aviso na tela (ao final da função) caso haja NAs ou valores iguais a zero em seu conjunto de dados
}
#
#após os testes lógicos para saber se a entrada de dados está correta, começar a função propriamente dita
#
#################################
##### Aplicação dos Modelos #####
#################################
#
#####################
###### Modelo 1 #####
#####################
#
# A proposta da função é fazer uma investigação inicial sobre a possível existência de dimorfismo intrassexual entre machos
# O primeiro modelo a ser aplicado é um modelo de regressão quadrática. Esse modelo tem como objetivo determinar se a relação entre tamanho do corpo e a característica de interesse é não linear
# A equação do modelo linear que será aplicada é: ln(Y) = a0 + a1*ln(X) + a2*ln(X^2) + E, na qual: ln = log natural, a = coeficientes de regressao ("alfa"), E = erro associado
#
dados$ln.x <- log(dados$x) #cria uma coluna no data.frame com os valores do log natural da variável x
dados$ln.y <- log(dados$y) #cria uma coluna no data.frame com os valores do log natural da variável y
modelo1 <- lm (ln.y ~ ln.x + I(ln.x^2), data=dados) #aplicação do modelo de regressão quadrática baseado na equação descrita acima
pvalue.a2 <- round(summary(modelo1)$coef[3,4],digits=4) #extrai o valor (arredondado) do p associado ao coeficiente a2 da equação
#valor de p associado ao coeficiente a2 está na posição coef[3,4] do sumário do modelo
#valor de p corresponde à probabilidade do coeficiente a2 ser igual zero (hipótese nula do modelo)
#
#################################
### análise do coeficiente a2 ###
#################################
#
#se o coeficiente a2 não for significantemente diferente de 0 (isto é, pvalue.a2 >= 0.05, hipótese nula aceita)
#conclui-se que a relação entre as variáveis não apresenta desvios significativos da linearidade
#portanto, não há descontinuidade nos traços morfológicos ou dimorfismo intrassexual
#obs: valor crítico de significância para este modelo = 0.05
#
if (pvalue.a2 >= 0.05) #teste lógico; se a probabilidade de a2 ser igual a 0 é maior ou igual a 0.05, conclui-se que a2 não é significantemente diferente de zero
{
#se a2 não é diferente de 0, testar novo modelo em que a relação entre as variáveis é linear
modelo0 <- lm (ln.y ~ ln.x, data=dados) #modelo para testar a relação entre os logs naturais das variáveis é linear
modelo0.r2 <- round(summary(modelo0)$adj.r.squared,digits=3) #guarda o valor de r2 ajustado do modelo (valor arredondado)
resulta <- list (summary(modelo1),summary(modelo0)) #salva um objeto com o sumário dos modelos 1 e do modelo nulo na forma de lista
cat("\n","Resumo da função:","\n","\n","\t","Modelo 1","\n","\t","\t","Equação do modelo: ln(Y) = a0 + a1*ln(X) + a2*ln(X^2) + E")
cat("\n","\t","\t","\t","Coeficiente a2 não é significantemente diferente de zero","\n","\t","\t","\t","Relação entre as variáveis é linear","\n","\t","\t","\t","Não foi detectado dimorfismo intrassexual")
cat("\n","\n","\t","\t","\t","P-value associado ao coeficiente a2=",pvalue.a2,"\n","\t","\t","\t","Valor de p é maior ou igual ao valor crítico (0.05)","\n","\t","\t","\t","Coeficiente não é significantemente diferente de zero")
cat("\n","\n","\n","\t","Modelo 0","\n","\t","\t","Equação do modelo: ln(Y) = a0 + a1*ln(X) + E")
cat("\n","\t","\t","\t","Este modelo linear descreve melhor os dados do que o modelo quadrático")
cat("\n","\n","\n","Seguem abaixo os sumários do Modelo 1 e do Modelo 0, respectivamente","\n","\n")
#acima:mensagens exibidas na tela para o usuário com os resultados da função
#
if (plot == "todos") #se usuário inserir "todos", ou não inserir nada (modo default), no argumento plot
{
#serão plotados gráficos de dispersão (y em função de x) E do modelo de regressão quadrática
par(mfrow=c(3,2), bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(ln.y~ln.x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 0: ln(Y) = a0 + a1*ln(X) + E",lab=c(7,5,5),xlim= range(dados$ln.x),ylim=range(dados$ln.y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
abline(modelo0, col="red",lty=1,lwd=1.5) #insere linha de tendência a partir do Modelo 0
r2 <- vector('expression',1) #cria um vetor de valor único, com o objeto "expressão"
r2[1] <- substitute(expression(italic(R)^2 == valor.r),list(valor.r = modelo0.r2))[2] #substitui o objeto do vetor anterior pelo r2 do modelo 0
legend("topleft",legend = r2,bty ='n',cex=0.7) #insere legenda no gráfico com o valor de r2
plot (modelo0) #plota os 4 gráficos padrão gerados pelo modelo 0
par(mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "resultado") #se usuário inserir "resultado" no argumento plot
{
#só será plotado o gráfico de dispersão (y em funcao de x)
par (mfrow=c(1,1),bty="l",mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(ln.y~ln.x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 0: ln(Y)= a0 + a1*ln(X) + E",lab=c(7,5,5),xlim= range(dados$ln.x),ylim=range(dados$ln.y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
abline(modelo0, col="red",lty=1,lwd=1.5) #insere linha de tendência a partir do Modelo 0
r2 <- vector('expression',1) #cria um vetor de valor único, com o objeto "expressão"
r2[1] <- substitute(expression(italic(R)^2 == valor.r),list(valor.r = modelo0.r2))[2] #substitui o objeto do vetor anterior pelo r2 do modelo 0
legend("topleft",legend = r2, bty = 'n',cex=0.7) #insere legenda no gráfico com o valor de r2
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "modelo") #se usuário inserir "modelo" no argumento plot
{
#só serão plotados os gráficos padrão resultantes do modelo linear
par (mfrow=c(2,2), bty="l",mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot (modelo0) #plota os 4 gráficos padrão gerados no modelo linear (modelo 0)
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "nenhum") #se usuário inserir "nenhum" no argumento plot
{
#não será plotado nenhum gráfico
}
}
#
#se o coeficiente a2 for significantemente diferente de 0 (isto é, pvalue.a2 < 0.05, hipótese nula rejeitada)
#
if (pvalue.a2 < 0.05) #teste lógico; se a probabilidade de a2 ser igual a 0 é menor que 0.05, conclui-se que a2 é significantemente diferente de zero
{
#
#####################
###### Modelo 2 #####
#####################
#
#se coeficiente a2 for significantemente diferente de zero, ou seja, a relação entre as variáveis não for linear,
#será feita nova análise (modelo 2) para determinar se existe um "switchpoint" a partir do qual a linearidade da relação entre as variáveis x e y torna-se descontínua
#novo modelo linear que será aplicado segue a seguinte equação: Y = b0 + b1*X + b2*(X-X0)*D + b3*D + E, na qual: b = coeficientes de regressao ("beta"), X0 = switchpoints a serem testados, E = erro associado, D = constante condicional do modelo (pré-definida como: D=0 para X < X0, D=1 para X > X0)
#
### Definição dos valores de switchpoint que serão testados ###
#os switchpoints devem ser valores contidos dentro da amplitude da variável x, ou seja, devem estar entre os valores mínimo e o máximo da variável
#o usuário escolhe o número de switchpoints que deseja testar (argumento "switchpoint" da função; default =100) e a função então calcula quais serão os valores testados
x0 <- round(seq(from=min(dados$x),to=max(dados$x),length.out=switchpoint),digits=3) #cria uma sequência de valores de switchpoint para teste
#para testar cada valor de x0 criado, serão feitas simulações da equação 2 com cada valor obtido
modelo2.r2 <- rep (x=NA, times=length(x0)) #cria vetor (de tamanho igual ao número de switchpoints testados) com NAs para guardar conjunto de valores da simulação
#
for (i in 1:length(x0)) #cria um loop para testar cada valor de switchpoint na equação do modelo 2 descrita acima
{
dados$x.x0 <- dados$x - x0[i] #cria coluna no data.frame com o valor da subtração de cada valor de X (cada linha) pelo valor de X0 (switchpoint simulado)
#para cálculo do D (D=1 para X > X0 e D=0 para X < X0),será usada uma operação lógica, na qual, para X > XO, retorna TRUE (valor 1) ou FALSE (valor 0)
dados$D <- as.numeric (dados$x > x0[i]) #cria coluna com resultado da operacao lógica entre parênteses (com coerção para número, i.e., valores 0 ou 1)
dados$D.x.x0 <- dados$D * dados$x.x0 #cria coluna com a multiplicação: (X-X0)*D
modelo2 <- lm (y ~ x + D.x.x0 + D, data=dados) #aplicação do modelo 2, baseado na equação descrita acima
modelo2.r2[i] <- summary(modelo2)$adj.r.squared #guarda o valor do r2 ajustado no objeto criado antes do loop
}
#
#para determinar qual switchpoint mais se adequa ao modelo, seleciona-se a simulação que apresenta maior valor de r2 ajustado
modelo2.max.r2 <- which(modelo2.r2 == max(abs(modelo2.r2))) #guarda a posição do elemento com maior valor (em módulo) de r2 ajustado, dentre todos os switchpoints simulados
#se o comando anterior selecionar mais do que um valor máximo, i.e., caso => selecionar apenas o primeiro valor com maior r2
modelo2.max.r2 <- modelo2.max.r2[1] #seleciona somente o primeiro valor do vetor anterior (caso existam duas simulações com valor de r2 máximo igual)
#comando anterior garante que só será selecionado um valor ideal de switchpoint para o modelo
modelo2.r2 <- round(modelo2.r2[modelo2.max.r2],digits=3) #seleciona o valor do maior r2 ajustado dentre os valores simulados
modelo2.x0 <- round(x0[modelo2.max.r2],digits=2) #extrai o valor de switchpoint que corresponde ao maior valor de r2 simulado
#
#################################
### analise do coeficiente b3 ###
#################################
#
#para o switchpoint com maior valor de r2 ajustado, será testado se o coeficiente b3 difere significativamente de zero
#para isso, aplica-se a equação novamente - somente para o melhor valor de switchpoint
dados$m2.x.x0 <- dados$x - modelo2.x0 #cria coluna no data.frame com o valor da subtração de cada valor de x pelo switchpoint (x0) selecionado
dados$m2.D <- as.numeric(dados$x > modelo2.x0) #cria coluna com resultado da operacao lógica entre parênteses (com coerção para número, i.e., valores 0 ou 1)
dados$m2.D.x.x0 <- dados$m2.D * dados$m2.x.x0 #cria coluna com a multiplicação: (X-X0)*D
modelo2.switchp <- lm (y ~ x + m2.D.x.x0 + m2.D, data=dados) #aplicação do modelo 2, baseado na equação descrita acima, para o melhor valor de switchpoint
pvalue.b3 <- round(summary(modelo2.switchp)$coef[4,4],digits=4) #extrai o valor (arredondado) do p associado ao coeficiente b3 da equação
#valor de p associado ao coeficiente b3 está na posição coef[4,4] do sumário do modelo
#valor de p corresponde à probabilidade do coeficiente b3 ser igual zero (hipótese nula do modelo)
#
#se o coeficiente b3 for significantemente diferente de 0 (pvalue < 0.05, hipótese nula rejeitada)
#conclui-se que o dimorfismo ocorre e que ele é descontínuo a partir do switchpoint encontrado (conforme figura 1a)
#obs: valor crítico de significância para este modelo = 0.05
#
if (pvalue.b3 < 0.05) #teste lógico; se a probabilidade de b3 ser igual a 0 é menor que 0.05, conclui-se que b3 é significantemente diferente de zero
{
b0 <- coef(modelo2.switchp)[1] #guarda o coeficiente b0 do modelo
b1 <- coef(modelo2.switchp)[2] #guarda o coeficiente b1 do modelo
b2 <- coef(modelo2.switchp)[3] #guarda o coeficiente b2 do modelo
b3 <- coef(modelo2.switchp)[4] #guarda o coeficiente b3 do modelo
resulta <- list (summary(modelo1),summary(modelo2.switchp)) #salva lista com sumário do modelo 1 e do modelo 2 para o switchpoint selecionado
cat("\n","Resumo da função:","\n","\n","\t","Modelo 1","\n","\t","\t","Equação do modelo: ln(Y) = a0 + a1*ln(X) + a2*ln(X^2) + E")
cat("\n","\t","\t","\t","Coeficiente a2 é significantemente diferente de zero","\n","\t","\t","\t","Relação entre as variáveis não é linear")
cat("\n","\n","\n","\t","Modelo 2","\n","\t","\t","Equação do modelo: Y = b0 + b1*X + b2*(X-X0)*D + b3*D + E")
cat("\n","\t","\t","\t","Coeficiente b3 é significantemente diferente de zero","\n","\t","\t","\t","Dimorfismo é descontínuo a partir do switchpoint")
cat("\n","\t","\t","\t","Switchpoint=",modelo2.x0,"\n","\t","\t","\t","R2 ajustado para melhor switchpoint=",modelo2.r2)
cat("\n","\n","\t","\t","\t","P-value associado ao coeficiente b3=",pvalue.b3,"\n","\t","\t","\t","Valor de p é menor que valor crítico (0.05)","\n","\t","\t","\t","Coeficiente é significantemente diferente de zero")
cat("\n","\n","\n","Seguem abaixo os sumários do Modelo 1 e do Modelo 2, respectivamente","\n","\n")
#acima:mensagens exibidas na tela para o usuário com os resultados da função
#
if (plot == "todos") #se usuário inserir "todos", ou não inserir nada (modo default), no argumento plot
{
#serão plotados gráficos de dispersão (y em função de x) E do modelo de regressão quadrática
par(mfrow=c(2,3),bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 2: Y= b0 + b1*X +b2*(X-X0)*D + b3*D + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
segments(x0=0,x1=modelo2.x0,y0=b0,y1=b0+(b1*modelo2.x0),col="red",lwd=1.5) #primeira reta (até o switchpoint) => relação linear entre as variáveis
segments(x0=modelo2.x0,x1=max(dados$x),y0=b0+(b1*modelo2.x0)+b3,y1=b0+b1*max(dados$x)+b2*(max(dados$x)-modelo2.x0)+b3,col="red",lwd=1.5) #segunda reta (a partir do switchpoint) => mostra o local de "quebra" da linearidade; há descontinuidade na relação entre as variáveis
abline(v=modelo2.x0,col="red",lty=3,lwd=1.5) #reta no valor do switchpoint
plot (modelo2.switchp) #plota os 4 gráficos padrão gerados pelo modelo 2
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "resultado") #se usuário inserir "resultado" no argumento plot
{
#será plotado somente o gráfico de dispersão (y em função de x)
par (mfrow=c(1,1),bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 2: Y= b0 + b1*X +b2*(X-X0)*D + b3*D + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
segments(x0=0,x1=modelo2.x0,y0=b0,y1=b0+(b1*modelo2.x0),col="red",lwd=1.5) #primeira reta (até o switchpoint) => relação linear entre as variáveis
segments(x0=modelo2.x0,x1=max(dados$x),y0=b0+(b1*modelo2.x0)+b3,y1=b0+b1*max(dados$x)+b2*(max(dados$x)-modelo2.x0)+b3,col="red",lwd=1.5) #segunda reta (a partir do switchpoint) => mostra o local de "quebra" da linearidade; há descontinuidade na relação entre as variáveis
abline(v=modelo2.x0,col="red",lty=3,lwd=1.5) #reta no valor do switchpoint
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "modelo") #se usuário inserir "modelo" no argumento plot
{
#so serão plotados os gráficos resultantes do modelo linear
par (mfrow=c(2,2), bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot (modelo2.switchp) #plota os 4 gráficos padrão gerados no modelo linear (modelo 2)
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "nenhum") #se usuário inserir "nenhum" no argumento plot
{
#não será plotado nenhum gráfico
}
}
#
#se o coeficiente b3 não for significantemente diferente de 0 (p >= 0.05, hipótese nula aceita)
#será feita uma última análise para testar se há mudança na linearidade da relação entre as variáveis a partir do switchpoint (ver figura 1b)
#
if (pvalue.b3 >= 0.05) #teste lógico; se a probabilidade de b3 ser igual a 0 é maior ou igual a 0.05, conclui-se que b3 não é significantemente diferente de zero
{
#
#####################
###### Modelo 3 #####
#####################
#
#aplicação do modelo 3 (testar se há mudança na linearidade da relação entre as variáveis a partir do switchpoint)
#equação do modelo: Y = b0 + b1*X + b2*(X-X0)*D + E
#para testar cada valor de x0, serão feitas simulações da equação 3 com cada valor de x0 já criado no modelo 2
modelo3.r2 <- rep (x=NA, times=length(x0)) #cria vetor (de tamanho igual ao número de switchpoints testados) com NAs para guardar conjunto de valores da simulação
#
for (i in 1:length(x0)) #cria um loop para testar cada valor de switchpoint na equação do modelo 3 descrita acima
{
dados$x.x0 <- dados$x - x0[i] #cria coluna no data.frame com o valor da subtração de cada valor de X (cada linha) pelo valor de X0 (switchpoint simulado)
#para cálculo do D (D=1 para X > X0 e D=0 para X < X0),será usada uma operação lógica, na qual, para X > XO, retorna TRUE (valor 1) ou FALSE (valor 0)
dados$D <- as.numeric(dados$x > x0[i]) #cria coluna com resultado da operacao lógica entre parênteses (com coerção para número, i.e., valores 0 ou 1)
dados$D.x.x0 <- dados$D * dados$x.x0 #cria coluna com a multiplicação: (X-X0)*D
modelo3 <- lm (y ~ x + I(D.x.x0), data=dados) #aplicação do modelo 3, baseado na equação descrita acima
modelo3.r2[i] <- summary(modelo3)$adj.r.squared #guarda o valor do r2 ajustado no objeto criado antes do loop
}
#
#para determinar qual switchpoint mais se adequa ao modelo, seleciona-se a simulação que apresenta maior valor de r2 ajustado
modelo3.max.r2 <- which(modelo3.r2 == max(abs(modelo3.r2))) #guarda a posição do elemento com maior valor (em módulo) de r2 ajustado, dentre todos os switchpoints simulados
#se o comando anterior selecionar mais do que um valor máximo, i.e., caso => selecionar apenas o primeiro valor com maior r2
modelo3.max.r2 <- modelo3.max.r2[1] #seleciona somente o primeiro valor do vetor anterior (caso existam duas simulações com valor de r2 máximo igual)
#comando anterior garante que só será selecionado um valor ideal de switchpoint para o modelo
#comando necessário em casos (raros, mas possíveis) de duas simulações apresentarem mesmo r2
modelo3.r2 <- round(modelo3.r2 [modelo3.max.r2],digits=3) #seleciona o valor do maior r2 ajustado dentre os valores simulados
modelo3.x0 <- round(x0[modelo3.max.r2],digits=2) #extrai o valor de switchpoint que corresponde ao maior valor de r2 simulado
#
#################################
### análise do coeficiente b2 ###
#################################
#
#para o switchpoint com maior valor de r2 ajustado, será testado se o coeficiente b3 difere significativamente de zero
#para isso, aplica-se a equação novamente - somente para o melhor valor de switchpoint
dados$m3.x.x0 <- dados$x - modelo3.x0 #cria coluna no data.frame com o valor da subtração de cada valor de x pelo switchpoint (x0) selecionado
dados$m3.D <- as.numeric(dados$x > modelo3.x0) #cria coluna com resultado da operacao lógica entre parênteses (com coerção para número, i.e., valores 0 ou 1)
dados$m3.D.x.x0 <- dados$m3.D * dados$m3.x.x0 #cria coluna com a multiplicação: (X-X0)*D
modelo3.switchp <- lm (y ~ x + m3.D.x.x0 + m3.D, data=dados) #aplicação do modelo 3, baseado na equação descrita acima, para o melhor valor de switchpoint
pvalue.b2 <- round(summary(modelo3.switchp)$coef[3,4],digits=4) #extrai o valor (arredondado) do p associado ao coeficiente b2 da equação
#valor de p associado ao coeficiente b2 está na posição coef[3,4] do sumário do modelo
#valor de p corresponde à probabilidade do coeficiente b2 ser igual zero (hipótese nula do modelo)
#
#se o coeficiente b2 for significantemente diferente de 0 (pvalue < 0.05, hipótese nula rejeitada)
#conclui-se que o dimorfismo ocorre e que há diferença na relação linear entre x e y a partir do swtichpoint (conforme figura 1b)
#obs: valor crítico de significância para este modelo = 0.05
#
#
if (pvalue.b2 < 0.05) #teste lógico; se a probabilidade de b2 ser igual a 0 é menor que 0.05, conclui-se que b2 é significantemente diferente de zero
{
b0 <- coef(modelo2.switchp)[1] #guarda o coeficiente b0 do modelo
b1 <- coef(modelo2.switchp)[2] #guarda o coeficiente b1 do modelo
b2 <- coef(modelo2.switchp)[3] #guarda o coeficiente b2 do modelo
resulta <- list (summary(modelo1),summary(modelo2.switchp),summary(modelo3.switchp)) #salva lista com sumário do modelo 1 e dos modelos 2 e 3 para o melhor switchpoint selecionado
cat("\n","Resumo da função:","\n","\n","\t","Modelo 1","\n","\t","\t","Equação do modelo: ln(Y) = a0 + a1*ln(X) + a2*ln(X^2) + E")
cat("\n","\t","\t","\t","Coeficiente a2 é significantemente diferente de zero","\n","\t","\t","\t","Relação entre as variáveis não é linear")
cat("\n","\n","\n","\t","Modelo 2","\n","\t","\t","Equação do modelo: Y = b0 + b1*X + b2*(X-X0)*D + b3*D + E")
cat("\n","\t","\t","\t","Coeficiente b3 não é significantemente diferente de zero","\n","\t","\t","\t","Dimorfismo não é descontínuo a partir do switchpoint")
cat("\n","\n","\n","\t","Modelo 3","\n","\t","\t","Equação do modelo: Y = b0 + b1*X + b2*(X-X0)*D + E")
cat("\n","\t","\t","\t","Coeficiente b2 é significantemente diferente de zero","\n","\t","\t","\t","Há diferença na relação linear entre x e y a partir do switchpoint")
cat("\n","\t","\t","\t","Switchpoint=",modelo3.x0,"\n","\t","\t","\t","R2 ajustado para melhor switchpoint=",modelo3.r2)
cat("\n","\n","\t","\t","\t","P-value associado ao coeficiente b2=",pvalue.b2,"\n","\t","\t","\t","Valor de p é menor que valor crítico (0.05)","\n","\t","\t","\t","Coeficiente é significantemente diferente de zero")
cat("\n","\n","\n","Seguem abaixo os sumários do Modelo 1, Modelo 2 e Modelo 3, respectivamente","\n","\n")
#acima:mensagens exibidas na tela para o usuário com os resultados da função
#
if (plot == "todos") #se usuário inserir "todos", ou não inserir nada (modo default), no argumento plot
{
#serão plotados gráficos de dispersão (y em função de x) E do modelo de regressão quadrática
par(mfrow=c(2,3),bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 3: Y= b0 + b1*X +b2*(X-X0)*D + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
segments(x0=0,x1=modelo3.x0,y0=b0,y1=b0+(b1*modelo3.x0),col="red",lwd=1.5) #primeira reta (até o switchpoint) => relação linear entre as variáveis
segments(x0=modelo3.x0,x1=max(dados$x),y0=b0+(b1*modelo3.x0),y1=b0+b1*max(dados$x)+b2*(max(dados$x)-modelo3.x0),col="red",lwd=1.5) #segunda reta (a partir do switchpoint) => mostra o local de "quebra" da linearidade; há descontinuidade na relação entre as variáveis
abline(v=modelo3.x0,col="red",lty=3,lwd=1.5) #reta no valor do switchpoint
plot (modelo3.switchp) #plota os 4 gráficos padrão gerados pelo modelo 3
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "resultado") #se usuário inserir "resultado" no argumento plot
{
#será plotado somente o gráfico de dispersão (y em função de x)
par (mfrow=c(1,1),bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 3: Y= b0 + b1*X +b2*(X-X0)*D + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
segments(x0=0,x1=modelo3.x0,y0=b0,y1=b0+(b1*modelo3.x0),col="red",lwd=1.5) #primeira reta (até o switchpoint) => relação linear entre as variáveis
segments(x0=modelo3.x0,x1=max(dados$x),y0=b0+(b1*modelo3.x0),y1=b0+b1*max(dados$x)+b2*(max(dados$x)-modelo3.x0),col="red",lwd=1.5) #segunda reta (a partir do switchpoint) => mostra o local de "quebra" da linearidade; há descontinuidade na relação entre as variáveis
abline(v=modelo3.x0,col="red",lty=3,lwd=1.5) #reta no valor do switchpoint
par(mfrow=c(1,1)) #volta ao padrao original do dispositivo grafico
}
if (plot == "modelo") #se usuário inserir "modelo" no argumento plot
{
#so serão plotados os gráficos resultantes do modelo linear
par (mfrow=c(2,2), bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot (modelo3.switchp) #plota os 4 gráficos padrão gerados no modelo linear (modelo 3)
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "nenhum") #se usuário inserir "nenhum" no argumento plot
{
#não será plotado nenhum gráfico
}
}
#
#se a probabilidade de b2 ser igual a zero for maior ou igual a 0.05 (p >= 0.05, hipótese nula aceita)
#conclui-se que o coeficiente b2 não é significantemente diferente de 0
#e portanto a relação entre as variáveis é linear (Y = b0 + b1*X)
#
if (pvalue.b2 >= 0.05) #teste lógico; se a probabilidade de b2 ser igual a 0 é maior ou igual a 0.05, conclui-se que b2 não é significantemente diferente de zero
{
modelo4 <- lm (y ~ x, data=dados) #se b2 não é diferente de 0, aplicar novo modelo em que a relação entre as variáveis é linear
modelo4.r2 <- round(summary(modelo4)$adj.r.squared,digits=3) #guarda o valor arredondado de r2 ajustado em um novo objeto
resulta <- list (summary(modelo3),summary(modelo4)) #salva um objeto com o sumário dos modelos 3 e 4 na forma de lista
cat("\n","Resumo da função:","\n","\n","\t","Modelo 4","\n","\t","\t","Equação do modelo: Y = b0 + b1*X + E")
cat("\n","\t","\t","\t","Este é o modelo que descreve melhor os dados")
cat("\n","\n","\n","Segue abaixo o sumário do Modelo 4","\n","\n")
#acima:mensagens exibidas na tela para o usuário com os resultados da função
#
if (plot == "todos") #se usuário inserir "todos", ou não inserir nada (modo default), no argumento plot
{
#serão plotados gráficos de dispersão (y em função de x) E do modelo de regressão quadrática
par(mfrow=c(2,3),bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 4: Y= b0 + b1*X + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
abline(modelo4, col="red",lty=1,lwd=1.5) #insere linha de tendência a partir do Modelo 4
r2 <- vector('expression',1) #cria um vetor de valor único, com o objeto "expressão"
r2[1] <- substitute(expression(italic(R)^2 == valor.r),list(valor.r = modelo4.r2))[2] #substitui o objeto do vetor anterior pelo r2 do modelo 0
legend("topleft",legend = r2,bty ='n',cex=0.7) #insere legenda no gráfico com o valor de r2
plot (modelo4) #plota os 4 gráficos padrão gerados pelo modelo 4
par(mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "resultado") #se usuário inserir "resultado" no argumento plot
{
#será plotado somente o gráfico de dispersão (y em função de x)
par (mfrow=c(1,1), bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot(y~x,data=dados,xlab="",ylab="",main="Relação entre variáveis x e y \n Modelo 4: Y= b0 + b1*X + E",lab=c(7,5,5),xlim= range(dados$x),ylim=range(dados$y)) #gráfico de x em função de y com vários parâmetros ajustados
mtext(c(eixo.x,eixo.y),side=c(1,2),line=1.5,family="serif",cex=0.8) #insere nome aos eixos x e y de acordo com os argumentos inseridos pelo usuário
abline(modelo4, col="red",lty=1,lwd=1.5) #insere linha de tendência a partir do Modelo 4
r2 <- vector('expression',1) #cria um vetor de valor único, com o objeto "expressão"
r2[1] <- substitute(expression(italic(R)^2 == valor.r),list(valor.r = modelo4.r2))[2] #substitui o objeto do vetor anterior pelo r2 do modelo 0
legend("topleft",legend = r2,bty ='n',cex=0.7) #insere legenda no gráfico com o valor de r2
plot (modelo4) #plota os 4 gráficos padrão gerados pelo modelo 4
par(mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "modelo") #se usuário inserir "modelo" no argumento plot
{
#so serão plotados os gráficos resultantes do modelo linear
par (mfrow=c(2,2), bty="l", mar=c(3,3,3,3),family="serif",mgp=c(1.5,0.4,0),tck=0.02,cex.axis=0.8,cex=0.8,cex.main=1,pch=20,col="black") #mudanças feitas no dispositivo gráfico e nos parâmetros dos gráficos
plot (modelo4) #plota os 4 gráficos padrão gerados pelo modelo 4)
par (mfrow=c(1,1)) #volta ao padrão original do dispositivo gráfico
}
if (plot == "nenhum") #se usuário inserir "nenhum" no argumento plot
{
#não será plotado nenhum gráfico
}
}
}
}
return (resulta)
}
### FIM ###
\\
==== Exemplos ====
Os exemplos abaixo são referentes a medidas((CManto ou ML = Comprimento do manto, em milímetros; Ptotal = Peso total, em gramas; Prep = Peso total do sistema reprodutor, em gramas )) obtidas em exemplares adultos de machos de lulas da espécie //Doryteuthis plei//.
**Link (.r)**::{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:exemplos_isd.analysis.r|}}
**Exemplo 1**:{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:exemplo1.cmanto.prep.txt|}}
**Exemplo 2**:{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:exemplo2.cmanto.ptotal.txt|}}
**Exemplo 3**:{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:exemplo3.cmanto.ptotal.txt|}}
\\
==== Anexos ====
**Arquivo I** Arquivo de Ajuda (.pdf)
{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:help_isd.analysis.pdf|}}
**Arquivo II** Modelos de Eberhard & Gutiérrez
{{:bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:modelos_de_eberhard_gutierrez.pdf|}}
**Figura 1** Gráficos dos Modelos 2 e 3, retirados de Eberhard & Gutiérrez, 1991
{{ :bie5782:01_curso_atual:alunos:trabalho_final:ligia.apostolico:figura.1_modelos_.png |}}
\\
==== Comentários finais ====
//Alguns comentários sobre a função final...//
A função foi modificada em relação à proposta inicial. Em primeiro lugar, resolvi remover os gráficos relacionados à variável x - tanto histogramas quanto qqplots ou gráficos de densidade cumulativa. Achei que seria desnecessário colocar estes gráficos como resultantes da função, uma vez que o usuário pode aplicá-los durante a análise exploratória dos dados (antes de executar a função //isd.analysis//) e assim ajustar todos os parâmetros que desejar. Acho que o objetivo da função está mais relacionado a aplicação dos modelos e aos gráficos resultantes desta aplicação.
Em segundo lugar, com relação ao critério utilizado para seleção do melhor modelo, mantive a proposta inicial (de utilizar o R2 ajustado), ao invés da sugestão do monitor Diogo (de utilizar algum critério baseado em informação, como o de Akaike). O modelo ((modelo de Eberhard & Gutiérrez,1991)) no qual me baseei para a função utiliza a informação do R2, assim como todos os artigos e análises que replicaram o modelo. Eu entendo que o R2 não seja o melhor parâmetro para comparar equações, mas decidi mantê-lo. Além do mais, não estou familiarizada com o uso de AIC, então achei melhor optar por um caminho mais seguro.
----//[[ligia.haselmann@ib.usp.br|Lígia H.A.]]//