source("config/setup.R")
#----------------------------------------------------------------------- # Carrega os pacotes necessários. library(lattice) library(latticeExtra) library(plyr) library(reshape2) library(car) library(nFactors) library(rpart) library(EACS)
#----------------------------------------------------------------------- # Análise exploratório dos dados. str(teca_qui) str(teca_crapar) str(teca_arv) # Juntar as inforções químicas, físico-hídricas e de produção. db <- merge(teca_qui, teca_crapar, all = TRUE) str(db) # Todos os locais tem pelo menos um registro. Usar o registro que # estiver disponível na ordem de prioridade das camadas: II, III e I. xtabs(~loc, data = na.omit(db)) # Número de missings por variável separado por camada. by(data = is.na(db[, -(1:2)]), INDICES = db$cam, FUN = colSums) #----------------------------------------------------------------------- # Fazendo a seleção de um valor para cada local. # Será mantido os valores na camada II de cada local. Caso a camada II # esteja incompleta, será usada a camada III e por fim a camada I. # Remove linhas com registros ausentes. dc <- na.omit(db) attr(dc, "na.action") <- NULL # Reordena os níveis (para usar na seleção). dc$cam <- factor(dc$cam, levels = levels(db$cam)[c(2, 3, 1)]) # Reordena. dc <- arrange(dc, loc, cam) # Pega o primeiro registro de cada local. dd <- do.call(rbind, by(dc, INDICES = dc$loc, FUN = function(x) { x[1, ]} )) str(dd) # Deixar o S positivo para facilitar a interpretação. dd$S <- -1 * dd$S
Na análise fatorial, a variável acc
(areia + cascalho + calhau) não
foi usada por ter alta correlação, por construção, com areia e
cascalho. A variável da curva de água do solo Ur
foi deixada de fora
para ser usado o cad = Ui - Ur
.
#----------------------------------------------------------------------- # Em cada linha o grupo de variáveis químicas, físicas e hídricas. j <- c("ph", "p", "k", "ca", "mg", "al", "ctc", "sat", "mo", "arg", "are", "cas", "alp", "n", "I", "Us", "Ui", "S", "cad") X <- dd[, j] fit0 <- factanal(X, factors = 4, # covmat = cov2cor(var(X)), # covmat = var(X), rotation = "varimax") print(fit0, digits = 2, cutoff = 0.5, sort = TRUE)
Com 4 fatores chegou-se a uma explicação superior à 70% da variância total. Baseado nos carregamentos de valor superior absoluto maior que 0.5, foi possível interpretar cada um dos fatores em termos de índices:
Algumas variáveis não apresentaram carregamento alto e por isso alta singularidade (uniquenesses), sendo por isso variáveis pouco explicadas pelos fatores. Essas variáveis serão removidas e o ajuste será refeito.
#----------------------------------------------------------------------- # Removendo variáveis de pouca importância ou sem explicação (com # alta singularidade/unicidade). # str(fit0) # Variáveis abandonadas. j[fit0$uniquenesses > 0.7] # Reajuste. X <- dd[, j[fit0$uniquenesses <= 0.7]] fit1 <- factanal(X, factors = 4, rotation = "varimax", scores = "regression") colnames(fit1$loadings) <- c("cation", "agua", "poros", "inclin") print(fit1, digits = 2, cutoff = 0.5, sort = TRUE) # Singularidade das variáveis. sort(fit1$uniquenesses) # Carregamentos do fator 1 contra 2. load <- fit1$loadings[, 1:2] plot(load, type = "n", xlab = "Cátions", ylab = "Água") abline(v = 0, h = 0, lty = 2) text(load, labels = names(X), cex = 0.8) # Pares de gráficos dos escores. sc <- fit1$scores pairs(sc)
Os resultados com o abandono das variáveis de alta singularidade não sofreu modificações substanciais. A interpretação dos fatores foi mantida. Com estes 4 fatores foi obtido uma explicação de 80% da variância total.
Apenas por precaução, foi confirmado por simulação que o número de fatores é de fato 4.
#----------------------------------------------------------------------- # Determinação por simulação do número de fatores. ev <- eigen(cor(X)) ap <- parallel(subject = nrow(X), var = ncol(X), rep = 100, cent = 0.05) nS <- nScree(x = ev$values, aparallel = ap$eigen$qevpea) plotnScree(nS)
library(psych) # Em cada linha o grupo de variáveis químicas, físicas e hídricas. j <- c("ph", "p", "k", "ca", "mg", "al", "ctc", "sat", "mo", "arg", "are", "cas", "alp", "n", "I", "Us", "Ui", "S", "cad") X <- dd[, j] fit0 <- fa(X, nfactors = 4, scores = TRUE, rotate = "varimax") print(fit0) biplot(fit0, labels = rownames(X), main = "")
Os escores ou índices definidos pela análise fatorial serão usados como variáveis explicativas da produção de madeira.
#----------------------------------------------------------------------- # Justar os escores com a variável de produção de madeira. de <- merge(cbind(loc = dd[, "loc"], as.data.frame(sc)), teca_arv[, c("loc", "prod")]) names(de)[2:5] <- colnames(fit1$loadings) str(de) xyplot(sqrt(prod) ~ cation + agua + poros + inclin, data = de, outer = TRUE, type = c("p", "r"), as.table = TRUE, scales = list(x = list(relation = "free")), xlab = "Escores da análise fatorial", ylab = "Produção de madeira") #----------------------------------------------------------------------- # Ajuste do modelo de regressão. # ATTENTION: Existe um maior atendimento dos pressupostos usando sqrt do # que a variável original. # m0 <- lm(prod ~ . - loc, data = de) m0 <- lm(sqrt(prod) ~ . - loc, data = de) summary(m0) par(mfrow = c(2, 2)) plot(m0) layout(1) # MASS::boxcox(m0) # abline(v = c(0.5, 1), col = 2) residualPlots(m0) im <- influence.measures(m0) summary(im) # Remover observação com maior alavancagem. which(im$is.inf[, "hat"]) r <- which.max(im$infmat[, "hat"]) m0 <- lm(prod ~ . - loc, data = de[-r, ]) # m0 <- lm(sqrt(prod) ~ . - loc, data = de[, 35]) summary(m0) par(mfrow = c(2, 2)) plot(m0) layout(1) residualPlots(m0) m1 <- update(m0, . ~ cation) summary(m1) anova(m0, m1) # Carregamento das variáveis no fator 1. cbind(fit1$loadings[, "cation"])
#----------------------------------------------------------------------- de <- merge(subset(dd, select = -cam), teca_arv[, c("loc", "prod")]) de <- subset(de, select = -loc) str(de) # Número de observações. nrow(de) # m0 <- lm(prod ~ ., data = de) m0 <- lm(sqrt(prod) ~ ., data = de) anova(m0) # im <- influence.measures(m0) # summary(im) # residualPlots(m0) # AIC. m1 <- step(m0) # Resumo do modelo selecionado. summary(m1) # residualPlots(m1) anova(m1, m0) im <- influence.measures(m1) summary(im) par(mfrow = c(2, 2)) plot(m1) layout(1)
#----------------------------------------------------------------------- # Árvore de regressão. layout(1) ar <- rpart(prod ~ ., data = de, method = "anova") summary(ar) rsq.rpart(ar) plot(ar, uniform = TRUE, main = "Árvore de regressão para produção") text(ar, use.n = TRUE, all = TRUE, cex = 0.8) pred <- factor(predict(ar)) plot(sat ~ mo, data = de, col = as.integer(pred), pch = NA) with(de, text(y = sat, x = mo, labels = rownames(de), col = as.integer(pred)))
cat(format(Sys.time(), format = "Atualizado em %d de %B de %Y.\n\n")) sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.