======Riguel F Contente======
{{:bie5782:01_curso2009:alunos:trabalho_final:Riguel.jpg}}
Doutorando em Oceanografia Biológica - IOUSP

====Proposta principal=====

Mudanças alométricas relacionadas ao tamanho corpóreo podem afetar comparações de 
determinadas grandezas morfométricas entre populações. Lleonart et al. 2000 
desenvolveram uma metodologia de remover esse efeito alométrico, a qual consiste na 
seguinte normalização: Z=Y(Tlo/Tl)^b, onde Y=valor da variável (e.g. tamanho da 
mandíbula), Tl= comprimento máximo do exemplar (em peixes, por exemplo, comprimento 
total), Tlo=representa um valor de referência (e.g. média ou mediana do comprimento 
total) através do qual todos as indivíduos são reduzidos. De um modo geral, o método 
remove do modelo o efeito do coeficiente 'a' ao qual se associa as informações de 
alometria. Pretendo realizar uma função que efetue essa normalização, produza gráficos 
de comparações de grandezas entre duas amostras (ex. boxplot, médias) e realize um teste 
de hipótese (teste-t ou anova) para avaliar possíveis diferenças significativas.

Referência

LLEONART J.; SALAT J.; TORRES G.J. Removing Allometric Effects of Body Size in Morphological Analysis Journal of Theoretical Biology, Volume 205, Number 1, 7 July 2000, pp. 85-93(9)

=== Comentários ===
== Paulo ==
  - Faça como o Fulano: formate a página no padrão que indicamos, po favor!
  - Legal a proposta, e bem dimensionada. Apenas lembre que um teste t ou ANOVA tm várias premissas, e seria legal pelo indicar ao usuário se estas premissas parecem se cumprir com os dados.

======Página de Ajuda======

<code>

RemAllom                package:nenhum                R Documentation



Teste de alometria e comparação de grandezas morfométricas com remoção do efeito da alometria 

Description:

Realiza teste-t para avaliar possíveis desvios do crescimento isométrico, realiza ANOVA para testar se há relação entre variável resposta e preditora, remove potenciais efeitos de alometria através de uma normalização e, utilizando ANOVA, compara simultaneamente médias das estruturas morfométricas entre duas amostras/populações/espécies. Produz gráficos dos modelos alométricos e box-plots de comparação entre as amostras/populações/espécies. Exporta as matrizes com dados normalizados para cada população/espécie para uso em posteriores análises multivariadas.

Usage:
	
	RemAllom<-function(A,B,coef.A,coef.B,nomes,diag.allom=TRUE,TL0media=TRUE,diag.norm=TRUE,WT=TRUE)


Arguments:

	A: data.frame ou matriz da pop/espécie A. Estruturas morfométricas nas colunas e indivíduos nas linhas 
	B: data.frame ou matriz da pop/espécie B. Mesmo formato 
	coef.A: vetor numérico contendo os coeficientes b esperados para cada variável em A
	coef.B: vetor numérico contendo os coeficientes b esperados para cada variável em B
	nomes: objeto contendo a denominação para cada amostra/população/espécie
	diag.allom: lógico. Plotar gráficos diagnósticos para avaliação das premissas de regressão linear?
	TL0media: lógico. Reduzir os valores pela média? Se FALSE, será pela mediana
	diag.norm: lógico. Plotar gráficos exploratórios para avaliação de normalidade e homogeneidade das variâncias?
	WT: lógico. Exportar em .txt matrizes com valores normalizados?
	
Details:

Através do test-t, testa-se a hipótese nula de que o coeficiente b do modelo alométrico log-transformado (ln) é igual ao valor esperado de acordo com a dimensão da variável morfométrica em questão. Para comparação das variáveis entre duas populações/espécies, os dados são normalizados pela média ou mediana do comprimento do tamanho corpóreo, segundo a expressão:Z=Y(Tlo/Tl)^b, onde Y=valor da variável, Tl= comprimento corpóreo máximo do exemplar (ou outra medida de referência), Tlo= valor de referência, média ou mediana do comprimento total. O efeito da espécie/população sobre o tamanho médio de cada variável é posteriormente testado pela análise de variância.

Value:

Gráficos das regressões e box-plots de comparações serão gerados. Um data frame é gerado apresentando o resultado do teste-t para alometria e da análise de variância para avaliação de associação entre variável resposta e independente. Outro data frame exibe a tabela ANOVA de comparação entre pop/espécies. Dois arquivos .txt serão gerados no diretório de trabalho. 

Warning:

Atente se os dados exibem distribuição normal, e se as variâncias das variáveis em comparação são homogêneas. Efetue os testes diagnósticos para explorar os dados.
O designe da ANOVA é para amostras balanceadas, logo NAs (que serão automaticamente removidos) podem tornar a análise inválida a subsequentes interpretações. 
Note que não há dimensão nos box-plot, pois Z é adimensional

Author(s):

Riguel Feltrin Contente

riguel@io.usp.br


References:

Slope of a regression line t-test - http://en.wikipedia.org/wiki/T-test
 
LLEONART J.; SALAT J.; TORRES G.J. Removing Allometric Effects of Body Size in Morphological Analysis. Journal of Theoretical Biology, Volume 205, Number 1, 7 July 2000, pp. 

85-93(9)

Examples:

data(iris) ## Abrir "iris"
setosa1=iris[iris$Species=="setosa",] ## Selecionar apenas os vetores com as variáveis morfométricas da sépala e pétala da espécie 'setosa' 
setosa<-setosa[,-5]
head(setosa)
virginica1=iris[iris$Species=="virginica",] ## Selecionar apenas os vetores com as variáveis morfométricas da sépala e pétala da espécie 'virginica'
virginica<-virginica[,-5]
head(virginica)
setosa.b<-c(1,1,1,1) ## Criar um vetor numérico com os coeficientes esperados de acordo com a dimensão da variável 
virginica.b<-c(1,1,1,1)
nomes<-c("Setosa","Virginica") ## Criar objeto com o nome das espécies
RemAllom(setosa,virginica,setosa.b,virginica.b,nomes,diag.allom=F,TL0media=T,diag.norm=F,WT=T) ## Rodar a função na qual a normalização será efetuada pela média da sépala

</code>

======Código da Função======

<code>


RemAllom<-function(A,B,coef.A,coef.B,nomes,diag.allom=T,TL0media=T,diag.norm=T,WT=T)
{
##Matriz A
##Remoção das linhas com NA, transformação ln e calculo da media e mediana
	x.1<-na.omit(A)
	x<-log(x.1)
	MM<- matrix(rep(NA,dim(x)[2]*2),nrow=2)
	for(j in 1:dim(x)[2])
	{
		MM[1,j]=mean(x[,j])
		MM[2,j]=median(x[,j])
	}
##Remoção das linhas com NA e calculo da media e mediana
	x.3<-na.omit(A)
	MM2<- matrix(rep(NA,dim(x.3)[2]*2),nrow=2)
	for(j in 1:dim(x.3)[2])
	{
		MM2[1,j]=mean(x.3[,j])
		MM2[2,j]=median(x.3[,j])
	}
	mm=MM2
	med=MM
##Teste de hipótese (test-t) para crescimento alométrico 
	Stat<- matrix(rep(NA,(dim(x)[2]*8)-8),nrow=8)
	for(j in 2:dim(x)[2])
	{
		som.q<-sum((x[,1]-med[1,j])^2)
		mod.l<-lm(x[,j]~x[,1],data=x)
		gl<-length(x[,1])-2
		test.coef<-((coef(mod.l)[2]-coef.A[j])*sqrt(gl))/(sqrt(sum(residuals(mod.l)^2)/som.q))
		p.1<-2*pt(test.coef,df=gl,lower.tail=F)
		Stat[1,c(j-1)]=round(test.coef,3)
		Stat[2,c(j-1)]=round(p.1,5)
		Stat[3,c(j-1)]=round(coef(mod.l)[2],3)
		Stat[4,c(j-1)]=coef.A[j]
		Stat[6,c(j-1)]<-round(summary(mod.l)$r.squared,3)
		Stat[7,c(j-1)]<-round(anova(mod.l)$"F value"[1],3)
		Stat[8,c(j-1)]<-round(anova(mod.l)$"Pr(>F)"[1],7)
		StatA<-data.frame(Stat)
	}
	for(i in 1:dim(StatA)[2])
	{
		if(StatA[2,i]>=0.05)
			{StatA[5,i]=c("isom")}
		else{
			if(StatA[2,i]<0.05&StatA[3,i]>StatA[4,i])
				{StatA[5,i]=c("alom+")}
			else{
				if(StatA[2,i]<0.05&StatA[3,i]<StatA[4,i])
					{StatA[5,i]=c("alom-")
		}}}
	}
	for(k in 1:dim(StatA)[2])
	{
		if(StatA[2,k]>=0.05)
			{StatA[2,k]=c(">0.05")}
		else{
			if(StatA[2,k]<0.05&StatA[2,k]>0.0001)
				{StatA[2,k]=StatA[2,k]}	
			else{		
				if(StatA[2,k]<0.0001)
					{StatA[2,k]=c("<0.0001")
		}}}
 		if(StatA[8,k]>=0.05)
			{StatA[8,k]=c(">0.05")}
		else{
			if(StatA[8,k]<0.05&StatA[8,k]>0.0001)
				{StatA[8,k]=StatA[8,k]}	
			else{		
				if(StatA[8,k]<0.0001)
					{StatA[8,k]=c("<0.0001")
		}}}
	}
	StatB<- matrix(rep(NA,(dim(x)[2]*5)-5),nrow=5)
	for(j in 2:dim(x)[2])
	{
		som.q<-sum((x[,1]-med[1,j])^2)
		mod.l<-lm(x[,j]~x[,1],data=x)
		gl<-length(x[,1])-2
		test.coef<-((coef(mod.l)[2]-coef.A[j])*sqrt(gl))/(sqrt(sum(residuals(mod.l)^2)/som.q))
		p.1<-2*pt(test.coef,df=gl,lower.tail=F)
		StatB[1,c(j-1)]=round(test.coef,3)
		StatB[2,c(j-1)]=round(p.1,5)
		StatB[3,c(j-1)]=coef(mod.l)[2]
		StatB[4,c(j-1)]=coef.A[j]
	}
##Gráficos diagnósticos para o Modelo Gaussiano 
	if(diag.allom==T)
	{
		for(j in 2:dim(x)[2])
		{
			win.graph(width = 9, height = 9)
			par(mfrow=c(2,2))
			mlw<-lm(x[,j]~x[,1],data=x)
			plot(mlw,main=paste("ln",names(x[j]),nomes[1]))
		}
		par(mfrow=c(1,1))
	}
##Regressões 
	g=1
	for(i in 1:dim(x)[2])
	{
		g=g+1
		if(g<=dim(x)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			mlw<-lm(x[,g]~x[,1],data=x)
			gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
			mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
			abline(mlw,col=2)
			g=g+1
			if(g<=dim(x)[2]){
				mlw<-lm(x[,g]~x[,1],data=x)
				gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
				mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
				abline(mlw,col=2)
				g = g+1
				if(g<=dim(x)[2]){
					mlw<-lm(x[,g]~x[,1],data=x)
					gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
					mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
					abline(mlw,col=2)
					g = g+1
					if(g<=dim(x)[2]){
						mlw<-lm(x[,g]~x[,1],data=x)
						gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
						mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
						abline(mlw,col=2)
		}}}}
	}
##Dados normalizados pela média ou mediana
	R.all<-matrix(rep(NA,dim(x.3)[1]*dim(x.3)[2]),nrow=dim(x.3)[1],ncol=dim(x.3)[2])
	for(l in 1:dim(x.3)[1])		
	{
		for(j in 2:dim(x.3)[2])
		{
			if(TL0media==T)
				{raz<-x.3[l,j]*((mm[1,1]/x.3[l,1])^StatB[3,j-1])
				R.all[l,j-1]<-raz}
			else{
				{raz<-x.3[l,j]*((mm[2,1]/x.3[l,1])^StatB[3,j-1])
				R.all[l,j-1]<-raz}}
		}
	}
##Resultados - dados matriz A	
	nomes.0<-c(nomes[1],"","","","","","","")
	nomes.1<-c("-->TESTE-T","p/ avaliar","alometria","","","-->ANOVA","p/ avaliar","sig. regressão")
	nomes.2<-c("t","Pr(>|t|)","b obs","b exp","cresc","r2","F","Pr(>F)")
	StatC<-data.frame(nomes.0,nomes.1,nomes.2,StatA)
	nomes.componentes<-c("sp/pop.","testes","comp.")
	nomes.col<-matrix(rep(NA,(dim(x)[2])-1),nrow=1)
	for(j in 2:dim(x)[2])
	{
		nomes.col[1,j-1]<-names(x[j])
	}
	MF<-R.all[,-c(dim(R.all))]
	colnames(MF)<-nomes.col
	names(StatC)<-c(nomes.componentes,nomes.col)

##Teste de normalidade- dados matriz A
	if(diag.norm==T)
		{x11()
		par(tcl=0.3)
		par(mfrow=c(2,2),pty="s")
		devAskNewPage(TRUE)
		for(m in 1:dim(MF)[2])
		{
			qqnorm(MF[,m],main=paste(nomes[1],nomes.col[1,m]))
			qqline(MF[,m],col=2)
			hist(MF[,m],prob=T,main="",xlab=nomes.col[1,m])
			par(new=T)
			curve(expr=dnorm(x,mean(MF[,m]),sd(MF[,m])),main="",col=2,ann=F,axes=F)
		}
		par(mfrow=c(1,1))
	}
	devAskNewPage(F)

##Matriz B
##Remoção das linhas com NA, transformação ln e calculo da media e mediana
	y.1<-na.omit(B)
	y<-log(y.1)
	MMy<- matrix(rep(NA,dim(y)[2]*2),nrow=2)
	for(j in 1:dim(y)[2])
	{
		MMy[1,j]=mean(y[,j])
		MMy[2,j]=median(y[,j])
	}
##Remoção das linhas com NA e calculo da media e mediana
	y.3<-na.omit(B)
	MM2y<- matrix(rep(NA,dim(y.3)[2]*2),nrow=2)
	for(j in 1:dim(y.3)[2])
	{
		MM2y[1,j]=mean(y.3[,j])
		MM2y[2,j]=median(y.3[,j])
   	}
##Teste de hipótese (test-t) para crescimento alométrico 
	mmy=MM2y
	medy=MMy
	Staty<- matrix(rep(NA,(dim(y)[2]*8)-8),nrow=8)
	for(j in 2:dim(y)[2])
	{
		som.qy<-sum((y[,1]-medy[1,j])^2)
		mod.ly<-lm(y[,j]~y[,1],data=y)
		gly<-length(y[,1])-2
		test.coefy<-((coef(mod.ly)[2]-coef.B[j])*sqrt(gly))/(sqrt(sum(residuals(mod.ly)^2)/som.qy))
		p.1y<-2*pt(test.coefy,df=gly,lower.tail=F)
		Staty[1,c(j-1)]=round(test.coefy,3)
		Staty[2,c(j-1)]=round(p.1y,5)
		Staty[3,c(j-1)]=round(coef(mod.ly)[2],3)
		Staty[4,c(j-1)]=coef.B[j]
		Staty[6,c(j-1)]<-round(summary(mod.ly)$r.squared,3)
		Staty[7,c(j-1)]<-round(anova(mod.ly)$"F value"[1],3)
		Staty[8,c(j-1)]<-round(anova(mod.ly)$"Pr(>F)"[1],7)
		StatAy<-data.frame(Staty)
	}
	for(i in 1:dim(StatAy)[2])
	{
		if(StatAy[2,i]>=0.05)
			{StatAy[5,i]=c("isom")}
		else{
			if(StatAy[2,i]<0.05&StatAy[3,i]>StatAy[4,i])
				{StatAy[5,i]=c("alom+")}
			else{
				if(StatAy[2,i]<0.05&StatAy[3,i]<StatAy[4,i])
					{StatAy[5,i]=c("alom-")
		}}}
	}
	for(k in 1:dim(StatAy)[2])
	{
		if(StatAy[2,k]>=0.05)
			{StatAy[2,k]=c(">0.05")}
		else{
			if(StatAy[2,k]<0.0499999&StatA[2,k]>0.0001)
				{StatAy[2,k]=StatAy[2,k]}	
			else{		
				if(StatAy[2,k]<0.0001)
					{StatAy[2,k]=c("<0.0001")
		}}}
		if(StatAy[8,k]>=0.05)
			{StatAy[8,k]=c(">0.05")}
		else{
			if(StatAy[8,k]<0.0499999&StatA[8,k]>0.0001)
				{StatAy[8,k]=StatAy[8,k]}	
			else{		
				if(StatAy[8,k]<0.0001)
					{StatAy[8,k]=c("<0.0001")
		}}}
	}
	StatBy<- matrix(rep(NA,(dim(y)[2]*5)-5),nrow=5)
	for(j in 2:dim(y)[2])
	{
		som.qy<-sum((y[,1]-medy[1,j])^2)
		mod.ly<-lm(y[,j]~y[,1],data=y)
		gly<-length(y[,1])-2
		test.coefy<-((coef(mod.ly)[2]-coef.B[j])*sqrt(gly))/(sqrt(sum(residuals(mod.ly)^2)/som.qy))
		p.1y<-2*pt(test.coefy,df=gly,lower.tail=F)
		StatBy[1,c(j-1)]=round(test.coefy,3)
		StatBy[2,c(j-1)]=round(p.1y,5)
		StatBy[3,c(j-1)]=coef(mod.ly)[2]
		StatBy[4,c(j-1)]=coef.B[j]
	}
##Gráficos diagnósticos para o Modelo Gaussiano 
	if(diag.allom==T)
	{
		for(j in 2:dim(y)[2])
		{
			win.graph(width=9,height=9)
			par(mfrow=c(2,2))
			mlwy<-lm(y[,j]~y[,1],data=y)
			plot(mlwy,main=paste("ln",names(x[j]),nomes[2]))
		}
		par(mfrow=c(1,1))
	}
##Regressões		
	g=1
	for(i in 1:dim(y)[2])
	{
		g=g+1
		if(g<=dim(y)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			mlwy<-lm(y[,g]~y[,1],data=y)
			grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
			mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
			abline(mlwy,col=2)
			g=g+1
			if(g<=dim(y)[2]){
				mlwy<-lm(y[,g]~y[,1],data=y)
				grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
				mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
				abline(mlwy,col=2)
				g = g+1
				if(g<=dim(y)[2]){
					mlwy<-lm(y[,g]~y[,1],data=y)
					grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
					mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
					abline(mlwy,col=2)
					g = g+1
					if(g<=dim(y)[2]){
						mlwy<-lm(y[,g]~y[,1],data=y)
						grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
						mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
						abline(mlwy,col=2)
		}}}}
	}
##Dados normalizados pela média ou mediana
	R.ally<-matrix(rep(NA,dim(y.3)[1]*dim(y.3)[2]),nrow=dim(y.3)[1],ncol=dim(y.3)[2])
	for(l in 1:dim(y.3)[1])		
	{
		for(j in 2:dim(y.3)[2])
		{
			if(TL0media==T)
				{razy<-y.3[l,j]*((mmy[1,1]/y.3[l,1])^StatBy[3,j-1])
				R.ally[l,j-1]<-razy}
			else{
				{razy<-y.3[l,j]*((mmy[2,1]/y.3[l,1])^StatBy[3,j-1])
				R.ally[l,j-1]<-razy}}
		}
	}	
	nomes.0y<-c(nomes[2],"","","","","","","")
	nomes.1y<-c("-->TESTE-T","p/ avaliar","alometria","","","-->ANOVA","p/ avaliar","sig. regressão")
	nomes.2y<-c("t","Pr(>|t|)","b obs","b exp","cresc","r2","F","Pr(>F)")
	StatCy<-data.frame(nomes.0y,nomes.1y,nomes.2y,StatAy)
	nomes.componentesy<-c("sp/pop.","testes","comp.")
	nomes.coly<-matrix(rep(NA,(dim(y)[2])-1),nrow=1)
	for(j in 2:dim(y)[2])
	{
		nomes.coly[1,j-1]<-names(y[j])
	}
	MFy<-R.ally[,-c(dim(R.ally))]
	colnames(MFy)<-nomes.coly
	names(StatCy)<-c(nomes.componentesy,nomes.coly)

##Teste de normalidade- dados matriz B
	if(diag.norm==T)
		{x11()
		par(tcl=0.3)
		par(mfrow=c(2,2),pty="s")
		devAskNewPage(TRUE)
		for(m in 1:dim(MFy)[2])
		{
			qqnorm(MFy[,m],main=paste(nomes[2],nomes.coly[1,m]))
			qqline(MFy[,m],col=2)
			hist(MFy[,m],prob=T,main="",xlab=nomes.coly[1,m])
			par(new=T)
			curve(expr=dnorm(x,mean(MFy[,m]),sd(MFy[,m])),main="",col=2,ann=F,axes=F)
		}
		par(mfrow=c(1,1))}
		devAskNewPage(F)

##Geração dos dados normalizados e análises comparativas
##Geração do fator
	MFA=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			MFA[i,j]=nomes[1]
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			MFA[i+dim(MF)[1],j]=nomes[2]
			as.factor(MFA)
			MFB<-MFA[,1]
	}}
##Matriz p/ as análises
	mtx=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mtx[i,j]=MF[i,j]
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mtx[i+dim(MF)[1],j]=MFy[i,j]
	}}	
##Exportação das matrizes com dados normalizados
	
	if(WT==T)
		{write.table(MF,file="MFA.txt",row.names=T,col.names=T,sep=",")
		write.table(MFy,file="MFB.txt",row.names=T,col.names=T,sep=",")
	}
	
##Box-plots
	g=0
	for(i in 1:dim(mtx)[2])
	{
		g=g+1
		if(g<=dim(mtx)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
			g=g+1
			if(g<=dim(mtx)[2]){
				boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
				g = g+1
				if(g<=dim(mtx)[2]){
					boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
					g = g+1
					if(g<=dim(mtx)[2]){
						boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
		}}}}
	}

##ANOVA
	
##SS Total
	media.geral=matrix(rep(NA,dim(MF)[2]),nrow=1)
	for(i in 1:dim(MF)[2])
	{
		media.geral[1,i]=mean(c(MF[,i],MFy[,i]))
	}
	mtx1=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mtx1[i,j]=(MF[i,j]-media.geral[1,j])^2
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mtx1[i+dim(MF)[1],j]=(MFy[i,j]-media.geral[1,j])^2
	}}	
	ss.total=(apply(mtx1,2,sum))

##Somatória Quadráticas Intra grupo
	mediaA.q<-matrix(rep(NA,(dim(MF)[1]*dim(MF)[2])),nrow=dim(MF)[1],ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mediaA.q[i,j]<-(MF[i,j]-mean(c(MF[,j])))^2
	}}
	ss.mediaA.q=apply(mediaA.q,2,sum)
	mediaB.q<-matrix(rep(NA,(dim(MFy)[1]*dim(MFy)[2])),nrow=dim(MFy)[1],ncol=dim(MFy)[2])
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mediaB.q[i,j]<-(MFy[i,j]-mean(c(MFy[,j])))^2
	}}
	ss.mediaB.q=apply(mediaB.q,2,sum)
	ss.intra=ss.mediaA.q+ss.mediaB.q
##Somatória Quadráticas Entre grupos
	mediaA<-apply(MF,2,mean)
	mediaB<-apply(MFy,2,mean)
	ss.entre=dim(MF)[1]*((mediaA-media.geral)^2+(mediaB-media.geral)^2)
##Cálculo do F
	GL1=1
	GL2=(dim(mtx)[1]-2)
	ms.entre=ss.entre/GL1
	ms.intra=ss.intra/GL2
	F.anova=ms.entre/ms.intra
	
	format.pval(c(0.1, 0.001, 1e-27))
	
##Cálculo do P
	pf.anova=matrix(rep(NA,dim(MF)[2]*1),ncol=dim(MF)[2],nrow=1)
	for(i in 1:dim(MF)[2])
	{		
		pf.anova[1,i]=pf(F.anova[1,i],GL1,GL2,lower.tail=FALSE)
	}
##Tabela ANOVA
	PPato=round(c(pf.anova),3)
	SSE=round(c(ms.entre),3)
	MSI=round(c(ms.intra),3)
	FI=round(c(F.anova),3)
	SSI=round(c(ss.intra),3)
	GL1m=rep(1,dim(MF)[2])
	GL2m=rep(GL2,dim(MF)[2])
	GL1.r=round(GL1m,1)
	GL2.r=round(GL2m,1)	
	Tabela.ANOVA<-matrix(c(SSE,SSI,GL1.r,GL2.r,SSE,MSI,FI,PPato),ncol=dim(MF)[2],nrow=8,byrow=T)
	colnames(Tabela.ANOVA)=nomes.col
	rownames(Tabela.ANOVA)=c("Sq.entre","Sq.res","gl1","gl2","Mq.entre","Mq.res","F","Pr(>F)")
##Resultado 
	cat("\t\nResultado teste-t - avaliação de desvios do crescimento isométrico \nResultado ANOVA - teste de associação entre as variáveis resposta e independente \nTabela ANOVA\n")
	return(StatC,StatCy,Tabela.ANOVA)
	
}	

</code>
======Arquivo da Função======
{{:bie5782:01_curso2009:alunos:trabalho_final:RemAllom.r|RemAllom.r}}