## Lendo os dados dados <- read.table("MalhaQuadBiss81.txt",header=T,sep="",quote="\"'",dec=",") head(dados) # Selecionando as covariaveis dados1 <- (dados[,c(4:10,14:16,20)]) head(dados1) dim(dados1) #[1] 81 11 nrow <- dim(dados1)[2] nrow #[1] 11 ncolu <- dim(dados1)[1] ncolu #[1] 81 ## Padronizando os dados medias=colMeans(dados1) medias # ph.CaCl2 mo p k ca mg H.Al # 4.6456790 1.7753086 7.7777778 0.1419753 2.2543210 0.9098765 7.4197531 # dg dp Ptotal Cota # 1.6059259 2.6120988 38.3856790 11.2756173 desvios=apply(dados1,2,sd) desvios # ph.CaCl2 mo p k ca mg H.Al #0.33393464 0.27043061 5.70745127 0.06565858 0.64769773 0.24269577 3.05362161 # dg dp Ptotal Cota #0.08204843 0.11018979 4.67799074 1.45096959 z=(t(dados1)-medias)/desvios z=t(z) head(z) # ph.CaCl2 mo p k ca mg #[1,] 1.3605087 0.09130386 -0.4866932 -1.40081168 0.5337042 -0.45273366 #[2,] 0.7615891 0.09130386 -0.3114837 -0.18238756 0.2249182 -0.04069516 #[3,] 1.0610489 -0.27847676 -0.8371123 -1.40081168 0.5337042 -0.04069516 #[4,] 1.6599685 -1.01803801 -0.3114837 -0.03008455 0.5337042 0.78338185 #[5,] 0.7615891 -0.27847676 -0.8371123 -1.24850867 0.0705252 -0.04069516 #[6,] 2.2588881 -2.12737987 -0.1362741 -0.94390264 1.3056693 1.19542035 # H.Al dg dp Ptotal Cota #[1,] -0.8906647 1.0246884 0.1624582 -0.59334855 -0.8446885 #[2,] -0.5304367 -0.3159832 2.3405184 1.43102485 -0.6172543 #[3,] -0.7269247 -2.0222926 -0.5635619 1.09541067 -0.5000913 #[4,] -1.0544047 -1.2910171 -0.3820569 0.66146368 -0.4518477 #[5,] -0.7269247 -0.3159832 -0.4728094 -0.01617767 -0.4794155 #[6,] -1.4146327 -1.9004133 0.7069732 1.63837883 -0.7619852 ## Calculando matriz de correlacao R <- cor(dados1,use='pairwise.complete.obs') R ## Calculando autovalores e autovetores autos <- eigen(R) autos #$values # [1] 3.414496869 2.877694991 1.847746566 0.960360329 0.611308441 0.406551478 # [7] 0.362102768 0.256434559 0.158020165 0.103826863 0.001456971 ## Calculando a matriz de cargas fatoriais n.fatores <- 3 L <- matrix(nrow=nrow,ncol=n.fatores) L for(i in 1:n.fatores) { L[,i] <- sqrt(autos$value[i])*autos$vectors[,i] } L # cargas iniciais: # [1,] 0.777975602 -0.2273087 0.4406629 # [2,] -0.088855420 -0.7742079 -0.3988603 # [3,] -0.007214629 -0.7417165 -0.2457221 # [4,] -0.379594781 -0.3246010 -0.3269165 # [5,] 0.788861067 -0.3915128 -0.2149584 # [6,] 0.710149620 -0.5121871 -0.2544245 # [7,] -0.737244591 0.3439813 -0.3527378 # [8,] -0.398828171 -0.6998604 0.3248647 # [9,] 0.463415012 0.2716178 -0.6225293 #[10,] 0.526592808 0.6184030 -0.5691822 #[11,] -0.579612006 -0.3019513 -0.5242216 ## Escores fatoriais obtidos pelo Metodo dos Minimos Quadrados Ponderados ## com as cargas fatoriais obtidas pelo Metodo das Componentes Principais f=solve(t(L)%*%L)%*%t(L)%*%t(z) t(f) # escores iniciais: # [,1] [,2] [,3] # [1,] 0.64005884 -0.219820433 1.33544235 # [2,] 1.03015074 0.599423516 -0.79348952 # Escores fatoriais obtidos pelo Metodo de Regressao f2=t(L)%*%solve(R)%*%t(z) t(f2) # [,1] [,2] [,3] # [1,] 0.64005884 -0.219820433 1.33544235 # [2,] 1.03015074 0.599423516 -0.79348952 ## Rotacionando as cargas fatoriais Lrot <- varimax(L) Lrot #$loadings # #Loadings: # [,1] [,2] [,3] [,1] [,2] [,3] # [1,] 0.876 0.276 # [2,] 0.177 -0.849 0.120 # [3,] 0.269 -0.711 0.182 # [4,] -0.258 -0.531 # [5,] 0.806 -0.254 -0.328 # [6,] 0.783 -0.386 -0.263 # [7,] -0.873 -0.125 # [8,] -0.409 0.765 # [9,] 0.133 -0.807 #[10,] 0.224 -0.965 #[11,] -0.487 -0.681 # # [,1] [,2] [,3] #SS loadings 3.221 2.502 2.417 #Proportion Var 0.293 0.227 0.220 #Cumulative Var 0.293 0.520 0.740 # #$rotmat # [,1] [,2] [,3] #[1,] 0.8628286 0.2202128 -0.4550090 #[2,] -0.4482546 0.7493778 -0.4873404 #[3,] 0.2336550 0.6244511 0.7452960 ## Selecionando apenas as cargas fatoriais Lrot2=matrix(as.vector(loadings(Lrot)),nrow=11,ncol=3) Lrot2 # [,1] [,2] [,3] [,1] [,2] [,3] # [1,] 0.87611484 0.27615257 0.08521516 0.87611484 # [2,] 0.17717956 -0.84881014 0.12046381 -0.84881014 # [3,] 0.26883866 -0.71085607 0.18161549 -0.71085607 # [4,] -0.25840704 -0.53098380 0.08726058 -0.53098380 # [5,] 0.80592317 -0.25390466 -0.32834649 0.80592317 # [6,] 0.78288003 -0.38631325 -0.26313658 0.78288003 # [7,] -0.87272584 -0.12484628 -0.09507716 -0.87272584 # [8,] 0.04550157 -0.40942482 0.76466105 0.76466105 # [9,] 0.13263669 -0.08314482 -0.80719694 -0.80719694 #[10,] 0.04416507 0.22395355 -0.96518644 -0.96518644 #[11,] -0.48724174 -0.68126440 0.02018143 -0.68126440 fchap=t(solve(t(Lrot2)%*%Lrot2)%*%t(Lrot2)%*%t(z)) fchap #escores iniciais # [,1] [,2] [,3] #[1,]0.96282937 0.810139056 0.81119476 #[2,]0.43474635 0.180551723 -1.35223571 #[3,]0.61302482 1.189843709 -0.66129312 plot(fchap) par(pty='s') eqscplot(Lrot2)