Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
pessoais:jcfaria [2007/03/01 11:34] jcfaria |
pessoais:jcfaria [2007/03/02 18:35] jcfaria |
||
---|---|---|---|
Linha 1: | Linha 1: | ||
- | ====== Página WIKI de José Cláudio Faria ====== | + | ====== José Cláudio Faria ====== |
- | {{pessoais:i_in_the_beach_2007.png|125X210}} | + | {{ pessoais:i_in_the_beach_2007.png}} |
- | Na Praia do Sul de Ilhéus/Bahia (janeiro de 2007) refletindo profundamente sobre o R!!! | + | Eu na Praia do Sul de Ilhéus/BA, em janeiro de 2007, refletindo profundamente sobre estatística computacional e o R!!! |
+ | |||
+ | Brincadeiras a parte... | ||
+ | |||
+ | **1. Quem sou** | ||
+ | |||
+ | - Engenheiro Agrônomo; | ||
+ | |||
+ | - Mestrado e Doutorado em Produção Vegetal pela Universidade Federal de Viçosa - UFV/MG. | ||
+ | |||
+ | |||
+ | **2. O que tenho feito profissionalmente** | ||
+ | |||
+ | - Professor de estatística e pesquisador da Universidade Estadual de Santa Cruz - UESC/BA; | ||
+ | |||
+ | - Coordenador e desenvolvedor do projeto Tinn-R (GUI/editor para o ambiente R). | ||
+ | |||
+ | |||
+ | **3. Sobre o R** | ||
+ | |||
+ | - Gostaria de tê-lo encontrado desde o início de minha carreira na área de estatística computacional. | ||
+ | |||
+ | |||
+ | **4. Sobre o futuro** | ||
+ | |||
+ | - Desejo aprofundar os conhecimentos em análise multivariada de dados no ambiente R; | ||
+ | |||
+ | - Aprimorar o Tinn-R e disponibilizá-lo também para a plataforma Linux; | ||
+ | |||
+ | - Trocar experiências com pessoas e equipes envolvidas nestas áreas. | ||
===== Tinn-R ===== | ===== Tinn-R ===== | ||
Linha 21: | Linha 50: | ||
==== Scripts ==== | ==== Scripts ==== | ||
+ | Todos os usuários estão automaticamente convidados a darem sugestões e alterarem contrutivamente todas as funções e scripts desta página. Solicito a gentileza de me enviar um [[joseclaudio.faria@terra.com.br | email]] comunicando as alterações. | ||
+ | |||
=== Introdução ao R === | === Introdução ao R === | ||
Linha 559: | Linha 590: | ||
} | } | ||
</code> | </code> | ||
+ | |||
+ | ==== Funções úteis ==== | ||
+ | === Tabelas e histogramas === | ||
+ | == Função tb.table == | ||
+ | |||
+ | Função simples, flexível mas poderosa para descrever, via tabela de distribuição de freqüências e histogramas, vetores e data.frames. | ||
+ | |||
+ | <code> | ||
+ | #=============================================================================== | ||
+ | # Name : tb.table | ||
+ | # Original author: José Cláudio Faria, Gabor Gothendievisk and Enio Jelihovschi | ||
+ | # Date (dd/mm/yy): 1/3/07 11:06:02 | ||
+ | # Version : v24 | ||
+ | # Aim : To make tables of frequency distribution and associated | ||
+ | # histogram | ||
+ | #=============================================================================== | ||
+ | # Arguments: | ||
+ | # breaks : Method to determine number of classes= c('Sturges', 'Scott', 'FD') | ||
+ | # by : Variable to group | ||
+ | # end : Last class (high value) | ||
+ | # h : Classes extent | ||
+ | # k : Class number | ||
+ | # right : Intervals right open (default = FALSE) | ||
+ | # start : First class (small value) | ||
+ | # x : A R object (vector or data.frame) | ||
+ | # histogram : Plot histogram (default = TRUE) | ||
+ | # title.histogram: Title of histogram c('auto', 'none') | ||
+ | #=============================================================================== | ||
+ | |||
+ | # Common functions | ||
+ | tb.make.table.I <- function(x, | ||
+ | start, | ||
+ | end, | ||
+ | h, | ||
+ | right, | ||
+ | histogram, | ||
+ | titleH) | ||
+ | { | ||
+ | f <- table(cut(x, br=seq(start, end, h), right=right)) # Absolut freq | ||
+ | fr <- f/length(x) # Relative freq | ||
+ | frP <- 100*(f/length(x)) # Relative freq, % | ||
+ | fac <- cumsum(f) # Cumulative freq | ||
+ | facP <- 100*(cumsum(f/length(x))) # Cumulative freq, % | ||
+ | fi <- round(f, 2) | ||
+ | fr <- round(as.numeric(fr), 2) | ||
+ | frP <- round(as.numeric(frP), 2) | ||
+ | fac <- round(as.numeric(fac), 2) | ||
+ | facP <- round(as.numeric(facP),2) | ||
+ | res <- data.frame(fi, fr, frP, fac, facP) # Make final table | ||
+ | names(res) <- c('Class limits', 'fi', 'fr', 'fr(%)', 'fac', 'fac(%)') | ||
+ | |||
+ | # Making the histogram: With Benilton suggestions | ||
+ | if (histogram) { | ||
+ | hist(x, | ||
+ | breaks = seq(start, end, h), | ||
+ | freq = T, | ||
+ | right = right, | ||
+ | xlab = 'Class limits', ylab='Frequency', | ||
+ | col = 'LightYellow', | ||
+ | main = titleH, | ||
+ | xlim = c(start, end), ylim=c(0, max(fi)), | ||
+ | las = 1, | ||
+ | xaxt = 'n') | ||
+ | axis(1, at=round(seq(start, end, h), 2)) | ||
+ | } | ||
+ | return(res) | ||
+ | } | ||
+ | |||
+ | tb.make.table.II <- function (x, | ||
+ | k, | ||
+ | breaks=c('Sturges', 'Scott', 'FD'), | ||
+ | right=FALSE, | ||
+ | histogram, | ||
+ | titleH) | ||
+ | { | ||
+ | x <- na.omit(x) | ||
+ | |||
+ | # User defines only x and/or 'breaks' | ||
+ | # (x, {k,}[breaks, right]) | ||
+ | if (missing(k)) { | ||
+ | brk <- match.arg(breaks) | ||
+ | switch(brk, | ||
+ | Sturges = k <- nclass.Sturges(x), | ||
+ | Scott = k <- nclass.scott(x), | ||
+ | FD = k <- nclass.FD(x)) | ||
+ | tmp <- range(x) | ||
+ | start <- tmp[1] - abs(tmp[2])/100 | ||
+ | end <- tmp[2] + abs(tmp[2])/100 | ||
+ | R <- end-start | ||
+ | h <- R/k | ||
+ | } | ||
+ | |||
+ | # User defines 'x' and 'k' | ||
+ | # (x, k,[breaks, right]) | ||
+ | else { | ||
+ | tmp <- range(x) | ||
+ | start <- tmp[1] - abs(tmp[2])/100 | ||
+ | end <- tmp[2] + abs(tmp[2])/100 | ||
+ | R <- end-start | ||
+ | h <- R/abs(k) | ||
+ | } | ||
+ | tbl <- tb.make.table.I(x, start, end, h, right, histogram, titleH) | ||
+ | return(tbl) | ||
+ | } | ||
+ | |||
+ | # With Gabor Grotendieck suggestions (thanks Gabor, very much!) | ||
+ | tb.table <- function(x, ...) UseMethod("tb.table") | ||
+ | |||
+ | # Table form vectors | ||
+ | tb.table.default <- function(x, | ||
+ | k, | ||
+ | start, | ||
+ | end, | ||
+ | h, | ||
+ | breaks=c('Sturges', 'Scott', 'FD'), | ||
+ | right=FALSE, | ||
+ | histogram=TRUE, | ||
+ | title.histogram=c('auto', 'none')) | ||
+ | { | ||
+ | # User defines nothing or not 'x' isn't numeric -> stop | ||
+ | stopifnot(is.numeric(x)) | ||
+ | x <- na.omit(x) | ||
+ | |||
+ | # User defines only 'x' | ||
+ | # (x, {k, start, end, h}, [breaks, right]) | ||
+ | if (missing(k) && missing(start) && missing(end) && missing(h) ) { | ||
+ | brk <- match.arg(breaks) | ||
+ | switch(brk, | ||
+ | Sturges = k <- nclass.Sturges(x), | ||
+ | Scott = k <- nclass.scott(x), | ||
+ | FD = k <- nclass.FD(x)) | ||
+ | tmp <- range(x) | ||
+ | start <- tmp[1] - abs(tmp[2])/100 | ||
+ | end <- tmp[2] + abs(tmp[2])/100 | ||
+ | R <- end-start | ||
+ | h <- R/k | ||
+ | } | ||
+ | |||
+ | # User defines 'x' and 'k' | ||
+ | # (x, k, {start, end, h}, [breaks, right]) | ||
+ | else if (missing(start) && missing(end) && missing(h)) { | ||
+ | stopifnot(length(k) >= 1) | ||
+ | tmp <- range(x) | ||
+ | start <- tmp[1] - abs(tmp[2])/100 | ||
+ | end <- tmp[2] + abs(tmp[2])/100 | ||
+ | R <- end-start | ||
+ | h <- R/abs(k) | ||
+ | } | ||
+ | |||
+ | # User defines 'x', 'start' and 'end' | ||
+ | # (x, {k,} start, end, {h,} [breaks, right]) | ||
+ | else if (missing(k) && missing(h)) { | ||
+ | stopifnot(length(start) >= 1, length(end) >=1) | ||
+ | tmp <- range(x) | ||
+ | R <- end-start | ||
+ | k <- sqrt(abs(R)) | ||
+ | if (k < 5) k <- 5 # min value of k | ||
+ | h <- R/k | ||
+ | } | ||
+ | |||
+ | # User defines 'x', 'start', 'end' and 'h' | ||
+ | # (x, {k,} start, end, h, [breaks, right]) | ||
+ | else if (missing(k)) { | ||
+ | stopifnot(length(start) >= 1, length(end) >= 1, length(h) >= 1) | ||
+ | } | ||
+ | |||
+ | else stop('Please, see the function sintax!') | ||
+ | |||
+ | if (histogram) { | ||
+ | x11() | ||
+ | par(mfrow=c(1, 1)) | ||
+ | title.histogram <- match.arg(title.histogram) | ||
+ | switch(title.histogram, | ||
+ | auto = titleH <- 'x', | ||
+ | none = titleH <- '') | ||
+ | } | ||
+ | tbl <- tb.make.table.I(x, start, end, h, right, histogram, titleH) | ||
+ | return(tbl) | ||
+ | } | ||
+ | |||
+ | # Table form data.frames | ||
+ | tb.table.data.frame <- function(df, | ||
+ | k, | ||
+ | by, | ||
+ | breaks=c('Sturges', 'Scott', 'FD'), | ||
+ | right=FALSE, | ||
+ | histogram=TRUE, | ||
+ | title.histogram=c('auto', 'none')) | ||
+ | { | ||
+ | stopifnot(is.data.frame(df)) | ||
+ | tmpList <- list() | ||
+ | nameF <- character() | ||
+ | nameY <- character() | ||
+ | |||
+ | # User didn't defines a factor | ||
+ | if (missing(by)) { | ||
+ | logCol <- sapply(df, is.numeric) | ||
+ | nHist <- length(logCol[logCol]) | ||
+ | if (histogram) { | ||
+ | count = 0 | ||
+ | if (nHist > 1) { | ||
+ | x11() | ||
+ | par(mfrow=c(4, 1)) | ||
+ | } | ||
+ | } | ||
+ | for (i in 1:ncol(df)) { | ||
+ | if (logCol[i]) { | ||
+ | count <- (count + 1) | ||
+ | if (count == 5) { | ||
+ | x11() | ||
+ | par(mfrow=c(4, 1)) | ||
+ | count <- 1 | ||
+ | } | ||
+ | title.histogram <- match.arg(title.histogram) | ||
+ | switch(title.histogram, | ||
+ | auto = titleH <- names(logCol[i]), | ||
+ | none = titleH <- '') | ||
+ | x <- as.matrix(df[ ,i]) | ||
+ | tbl <- tb.make.table.II(x, k, breaks, right, histogram, titleH) | ||
+ | tmpList <- c(tmpList, list(tbl)) | ||
+ | } | ||
+ | } | ||
+ | valCol <- logCol[logCol] | ||
+ | names(tmpList) <- names(valCol) | ||
+ | return(tmpList) | ||
+ | } | ||
+ | |||
+ | # User defines one factor | ||
+ | else { | ||
+ | namesdf <- names(df) | ||
+ | pos <- which(namesdf == by) | ||
+ | stopifnot(is.factor((df[[pos]]))) | ||
+ | nF <- table(df[[pos]]) | ||
+ | logCol <- sapply(df, is.numeric) | ||
+ | nHist <- length(logCol[logCol]) | ||
+ | nDisGraph <- round((length(nF) * nHist) / 12) # 12 is the maximum easily visible | ||
+ | if (histogram) { | ||
+ | count <- 0 | ||
+ | x11() | ||
+ | par(mfrow=c(4, 3)) | ||
+ | } | ||
+ | for(i in 1:length(nF)) { | ||
+ | tmpdf <- subset(df, df[[pos]] == names(nF[i])) | ||
+ | logCol <- sapply(tmpdf, is.numeric) | ||
+ | for (j in 1:ncol(tmpdf)) { | ||
+ | if (logCol[j]) { | ||
+ | count <- (count + 1) | ||
+ | if (count == 13) { | ||
+ | x11() | ||
+ | par(mfrow=c(4, 3)) | ||
+ | count <- 1 | ||
+ | } | ||
+ | nameF <- names(nF[i]) | ||
+ | nameY <- names(logCol[j]) | ||
+ | nameFY <- paste(nameF,'.', nameY, sep="") | ||
+ | title.histogram <- match.arg(title.histogram) | ||
+ | switch(title.histogram, | ||
+ | auto = titleH <- nameFY, | ||
+ | none = titleH <- '') | ||
+ | x <- as.matrix(tmpdf[ ,j]) | ||
+ | tbl <- tb.make.table.II(x, k, breaks, right, histogram, titleH) | ||
+ | newFY <- list(tbl) | ||
+ | names(newFY) <- sub(' +$', '', nameFY) | ||
+ | tmpList <- c(tmpList, newFY) | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | return(tmpList) | ||
+ | } | ||
+ | </code> | ||
+ | |||
+ | == Testar função tb.table == | ||
+ | O script abaixo possibilita testar e aprender a usar a função tb.table. | ||
+ | |||
+ | <code> | ||
+ | #=============================================================================== | ||
+ | # Name : tb.table_test | ||
+ | # Original author: Jose Cláudio Faria | ||
+ | # Date (dd/mm/yy): 1/3/07 11:06:02 | ||
+ | # Version : v24 | ||
+ | # Aim : To learn how to use the function tb.table | ||
+ | #=============================================================================== | ||
+ | # Observation : Test it line by line | ||
+ | #=============================================================================== | ||
+ | # 1.Tables | ||
+ | # 1.1. Tables from vectors | ||
+ | #=============================================================================== | ||
+ | |||
+ | ## To debug | ||
+ | # mtrace.off() | ||
+ | # mtrace(tb.make.table.I) | ||
+ | # mtrace(tb.make.table.II) | ||
+ | # mtrace(tb.table.default) | ||
+ | # mtrace(tb.table.data.frame) | ||
+ | |||
+ | # Make a vector | ||
+ | set.seed(1) | ||
+ | x=rnorm(150, 5, 1) | ||
+ | |||
+ | tb.table(x, his=F) | ||
+ | tb.table(x) | ||
+ | tb.table(x, title.his='none') | ||
+ | tb.table(x, k=10, his=T) | ||
+ | |||
+ | #Title | ||
+ | tb.table(x, title.his='teste') #error! | ||
+ | tb.table(x, title.his='none') | ||
+ | tb.table(x, title.his='auto') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(x, breaks='Sturges') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(x, breaks='St') | ||
+ | |||
+ | tb.table(x, breaks='Scott') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(x, b='Sc') | ||
+ | |||
+ | tb.table(x, breaks='FD') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(x, breaks='F') | ||
+ | |||
+ | tb.table(x, breaks='F', right=T) | ||
+ | |||
+ | # Will make a error! | ||
+ | tb.table(x, breaks='S') #('S'turges) and ('S'cott) | ||
+ | |||
+ | tb.table(x, k=4) | ||
+ | |||
+ | tb.table(x, k=20) | ||
+ | |||
+ | # Partial | ||
+ | tb.table(x, start=4, end=6) # Will make error! | ||
+ | tb.table(x, start=4, end=6, his=F) | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(x, s=4, e=6, his=F) | ||
+ | |||
+ | # Partial | ||
+ | tb.table(x, start=4.5, end=5.5, his=F) | ||
+ | |||
+ | # Partial | ||
+ | tb.table(x, start=5, end=6, h=.5, his=F) | ||
+ | |||
+ | # Nonsense | ||
+ | tb.table(x, start=0, end=10, h=.5) | ||
+ | |||
+ | # First and last class forced (fi=0) | ||
+ | tb.table(x, start=1, end=9, h=1) | ||
+ | |||
+ | tb.table(x, start=1, end=10, h=2) | ||
+ | |||
+ | |||
+ | #=============================================================================== | ||
+ | # 1.2. Tables from data.frames | ||
+ | #=============================================================================== | ||
+ | # Make a data.frame | ||
+ | mdf=data.frame(X1 =rep(LETTERS[1:4], 25), | ||
+ | X2 =as.factor(rep(1:10, 10)), | ||
+ | Y1 =c(NA, NA, rnorm(96, 10, 1), NA, NA), | ||
+ | Y2 =rnorm(100, 60, 4), | ||
+ | Y3 =rnorm(100, 50, 4), | ||
+ | Y4 =rnorm(100, 40, 4)) | ||
+ | |||
+ | tb.table(mdf) | ||
+ | |||
+ | tb.table(mdf, title.his='none') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(mdf, breaks='Sturges') | ||
+ | |||
+ | # Equal to above | ||
+ | tb.table(mdf, breaks='St') | ||
+ | |||
+ | tb.table(mdf, breaks='Scott') | ||
+ | |||
+ | tb.table(mdf, breaks='FD') | ||
+ | |||
+ | tb.table(mdf, k=4) | ||
+ | |||
+ | tb.table(mdf, k=10) | ||
+ | |||
+ | levels(mdf$X1) | ||
+ | tbl = tb.table(mdf, k=5, by='X1') | ||
+ | length(tbl) | ||
+ | names(tbl) | ||
+ | tbl | ||
+ | |||
+ | tb.table(mdf, breaks='FD', by='X1') | ||
+ | |||
+ | # A 'big' result: X2 is a factor with 10 levels! | ||
+ | tb.table(mdf, breaks='FD', by='X2') | ||
+ | |||
+ | tb.table(mdf, breaks='FD', k=5, by='X2') | ||
+ | |||
+ | tb.table(iris, k=5) | ||
+ | |||
+ | tb.table(iris, k=10) | ||
+ | |||
+ | levels(iris$Species) | ||
+ | tbl=tb.table(iris, k=5, by='Species') | ||
+ | length(tbl) | ||
+ | names(tbl) | ||
+ | tbl | ||
+ | |||
+ | tb.table(iris, k=5, by='Species', right=T) | ||
+ | |||
+ | tb.table(iris, breaks='FD', by='Species') | ||
+ | |||
+ | library(MASS) | ||
+ | levels(Cars93$Origin) | ||
+ | tbl=tb.table(Cars93, k=5, by='Origin') | ||
+ | names(tbl) | ||
+ | tbl | ||
+ | |||
+ | tb.table(Cars93, breaks='FD', by='Origin') | ||
+ | </code> | ||
+ | |||
+ | === Superfície de resposta === | ||
+ | == Função plotlm3d == | ||
+ | The simple, power and very flexible function **plotlm3d** enables you to plot 3d points and/or surfaces obtained from linear methods. It was adapted from scatter3d [[http://socserv.socsci.mcmaster.ca/jfox/Misc/Rcmdr/index.html | Rcmdr package]] of John Fox and some [[ http://www.stat.wisc.edu/~deepayan | Deepayan Sarkar]] ideas. | ||
+ | |||
+ | It requires the **rgl** package that you can download from [[http://cran.r-project.org|CRAN]]. | ||
+ | |||
+ | <code> | ||
+ | #=============================================================================== | ||
+ | # Name : plotlm3d | ||
+ | # Original author: John Fox (scatter3d from package Rcmdr) | ||
+ | # Changes : Jose Claudio Faria and Duncan Murdoch | ||
+ | # Date (dd/mm/yy): 12/8/06 19:44:37 | ||
+ | # Version : v18 | ||
+ | # Aim : To plot 3d scatter, an or, surfaces with rgl package | ||
+ | #=============================================================================== | ||
+ | |||
+ | # Arguments: | ||
+ | # x variable for horizontal axis. | ||
+ | # y variable for out-of-screen axis. | ||
+ | # z variable for vertical axis (response). | ||
+ | # surface plot surface(s) (TRUE or FALSE). | ||
+ | # model one or more linear model to fit ('z ~ x + y' is the default). | ||
+ | # groups if NULL (the default), no groups are defined; if a factor, | ||
+ | # a different surface or set of surfaces is plotted for each | ||
+ | # level of the factor; in this event, the colours in plane.col | ||
+ | # are used successively for the points and surfaces. | ||
+ | # model.by.group if TRUE the function will adjust one model for each level | ||
+ | # of groups; the order of the models must be the same of the | ||
+ | # level of the. | ||
+ | # model.summary print summary or summaries of the model(s) fit (TRUE or FALSE). | ||
+ | # simple.axes whether to draw sinple axes (TRUE or FALSE). | ||
+ | # box whether to draw a box (TRUE or FALSE). | ||
+ | # xlab, | ||
+ | # ylab, | ||
+ | # zlab axis labels. | ||
+ | # surface.col vector of colours for regression planes, used in the order | ||
+ | # specified by fit. | ||
+ | # point.col colour of points. | ||
+ | # grid.col colour of grid lines on the regression surface(s). | ||
+ | # grid plot grid lines on the regression surface(s) (TRUE or FALSE). | ||
+ | # grid.lines number of lines (default, 26) forming the grid, in each of | ||
+ | # the x and z directions. | ||
+ | # sphere.factor relative size factor of spheres representing points; the | ||
+ | # default size is dependent on the scale of observations. | ||
+ | # threshold if the actual size of the spheres is less than the threshold, | ||
+ | # points are plotted instead. | ||
+ | # speed revolutions of the plot per second. | ||
+ | # revolutions number of full revolutions of the display. | ||
+ | |||
+ | plotlm3d <- function (x, y, z, | ||
+ | surface = T, | ||
+ | model = 'z ~ x + y', | ||
+ | groups = NULL, | ||
+ | model.by.group = F, | ||
+ | model.summary = F, | ||
+ | simple.axes = T, | ||
+ | box = F, | ||
+ | xlab = deparse(substitute(x)), | ||
+ | ylab = deparse(substitute(y)), | ||
+ | zlab = deparse(substitute(z)), | ||
+ | surface.col = c('blue', 'orange', 'red', 'green', | ||
+ | 'magenta', 'cyan', 'yellow', 'gray', 'brown'), | ||
+ | point.col = 'yellow', | ||
+ | grid.col = material3d("color"), | ||
+ | grid = T, | ||
+ | grid.lines = 26, | ||
+ | sphere.factor = 1, | ||
+ | threshold = 0.01, | ||
+ | speed = 0.5, | ||
+ | revolutions = 0) | ||
+ | { | ||
+ | require(rgl) | ||
+ | require(mgcv) | ||
+ | summaries <- list() | ||
+ | |||
+ | if ((!is.null(groups)) && model.by.group) | ||
+ | if (!nlevels(groups) == length(model)) | ||
+ | stop('Model number is different of the number of groups') | ||
+ | |||
+ | if ((!is.null(groups)) && (nlevels(groups) > length(surface.col))) | ||
+ | stop('Number of groups exceeds number of colors') | ||
+ | |||
+ | if ((!is.null(groups)) && (!is.factor(groups))) | ||
+ | stop('groups variable must be a factor.') | ||
+ | |||
+ | xlab; ylab; zlab | ||
+ | |||
+ | valid <- if (is.null(groups)) | ||
+ | complete.cases(x, y, z) | ||
+ | else | ||
+ | complete.cases(x, y, z, groups) | ||
+ | |||
+ | x <- x[valid] | ||
+ | y <- y[valid] | ||
+ | z <- z[valid] | ||
+ | | ||
+ | if (!is.null(groups)) | ||
+ | groups <- groups[valid] | ||
+ | |||
+ | levs <- levels(groups) | ||
+ | size <- max(c(x,y,z))/100 * sphere.factor | ||
+ | |||
+ | if (is.null(groups)) { | ||
+ | if (size > threshold) | ||
+ | spheres3d(x, y, z, color = point.col, radius = size) | ||
+ | else | ||
+ | points3d(x, y, z, color = point.col) | ||
+ | } | ||
+ | else { | ||
+ | if (size > threshold) | ||
+ | spheres3d(x, y, z, color = surface.col[as.numeric(groups)], radius = size) | ||
+ | else | ||
+ | points3d(x, y, z, color = surface.col[as.numeric(groups)]) | ||
+ | } | ||
+ | |||
+ | aspect3d(c(1, 1, 1)) | ||
+ | |||
+ | if (surface) { | ||
+ | xvals <- seq(min(x), max(x), length = grid.lines) | ||
+ | yvals <- seq(min(y), max(y), length = grid.lines) | ||
+ | | ||
+ | dat <- expand.grid(x = xvals, y = yvals) | ||
+ | |||
+ | for (i in 1:length(model)) { | ||
+ | if (is.null(groups)) { | ||
+ | mod <- lm(formula(model[i])) | ||
+ | |||
+ | if (model.summary) | ||
+ | summaries[[model[i]]] <- summary(mod) | ||
+ | |||
+ | zhat <- matrix(predict(mod, newdata = dat), grid.lines, grid.lines) | ||
+ | surface3d(xvals, yvals, zhat, color = surface.col[i], alpha = 0.5, lit = F) | ||
+ | |||
+ | if (grid) | ||
+ | surface3d(xvals, yvals, zhat, color = grid.col, alpha = 0.5, | ||
+ | lit = F, front = 'lines', back = 'lines') | ||
+ | } | ||
+ | else { # groups is not NULL | ||
+ | if (!model.by.group) { | ||
+ | for (j in 1:length(levs)) { | ||
+ | mod <- lm(formula(model[i]), subset = (groups == levs[j])) | ||
+ | |||
+ | if (model.summary) | ||
+ | summaries[[paste(model[i], '.', levs[j], sep = '')]] <- summary(mod) | ||
+ | |||
+ | zhat <- matrix(predict(mod, newdata = dat), grid.lines, grid.lines) | ||
+ | surface3d(xvals, yvals, zhat, color = surface.col[j], alpha = 0.5, lit = F) | ||
+ | |||
+ | if (grid) | ||
+ | surface3d(xvals, yvals, zhat, color = grid.col, alpha = 0.5, | ||
+ | lit = F, front = 'lines', back = 'lines') | ||
+ | |||
+ | texts3d(min(x), min(y), predict(mod, newdata = data.frame(x = min(x), y = min(y), | ||
+ | groups = levs[j])), paste(levs[j], ' '), adj = 1, color = surface.col[j]) | ||
+ | } | ||
+ | } | ||
+ | else { # model.by.group is TRUE | ||
+ | mod <- lm(formula(model[i]), subset = (groups == levs[i])) | ||
+ | |||
+ | if (model.summary) | ||
+ | summaries[[paste(model[i], '.', levs[i], sep = '')]] <- summary(mod) | ||
+ | |||
+ | zhat <- matrix(predict(mod, newdata = dat), grid.lines, grid.lines) | ||
+ | |||
+ | surface3d(xvals, yvals, zhat, color = surface.col[i], alpha = 0.5, lit = F) | ||
+ | |||
+ | if (grid) | ||
+ | surface3d(xvals, yvals, zhat, color = grid.col, alpha = 0.5, | ||
+ | lit = F, front = 'lines', back = 'lines') | ||
+ | |||
+ | texts3d(min(x), min(y), predict(mod, newdata = data.frame(x = min(x), y = min(y), | ||
+ | groups = levs[i])), paste(levs[i], ' '), adj = 1, color = surface.col[i]) | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | if(simple.axes) { | ||
+ | axes3d(c('x', 'y', 'z')) | ||
+ | title3d(xlab = xlab, ylab = ylab, zlab = zlab) | ||
+ | } | ||
+ | else | ||
+ | decorate3d(xlab = xlab, ylab = ylab, zlab = zlab, box = box) | ||
+ | |||
+ | if (revolutions > 0) { | ||
+ | start <- proc.time()[3] | ||
+ | startMatrix <- par3d("userMatrix") | ||
+ | while ((theta <- speed*(proc.time()[3] - start))/2/pi < revolutions) { | ||
+ | rgl.viewpoint(userMatrix = rotate3d(startMatrix, theta, 0, 0, 1)) | ||
+ | } | ||
+ | } | ||
+ | if (model.summary) | ||
+ | return(summaries) | ||
+ | else | ||
+ | return(invisible(NULL)) | ||
+ | } | ||
+ | </code> | ||
+ | |||
+ | == Testar função tb.table == | ||
+ | <code> | ||
+ | #=============================================================================== | ||
+ | # Name : Script to test plotlm3d | ||
+ | # Author : Jose Claudio Faria and Duncan Murdoch | ||
+ | # Date (dd/mm/yy): 23/7/06 07:21:23 | ||
+ | # Version : v17 | ||
+ | # Aim : To plot 3d scatter, an or, surfaces with rgl package | ||
+ | #=============================================================================== | ||
+ | |||
+ | # mtrace(plotlm3d) | ||
+ | # mtrace.off | ||
+ | |||
+ | # Example 1 | ||
+ | open3d() | ||
+ | rgl.bringtotop(stay = T) | ||
+ | with(iris, plotlm3d(Sepal.Length, Sepal.Width, Petal.Length, | ||
+ | surface = F, | ||
+ | groups = Species, | ||
+ | xlab = 'SL', | ||
+ | ylab = 'PL', | ||
+ | zlab = 'SW', | ||
+ | grid = F, | ||
+ | sphere.factor = 1)) | ||
+ | |||
+ | # Example 2 | ||
+ | open3d() | ||
+ | rgl.bringtotop(stay = T) | ||
+ | with(iris, plotlm3d(Sepal.Length,Sepal.Width, Petal.Length, | ||
+ | model = c('z ~ x + y', | ||
+ | 'z ~ x + y + I(x^2) + I(y^2) + I(x*y)'), | ||
+ | surface = T, | ||
+ | groups = Species, | ||
+ | simple.axes = F, | ||
+ | box = T, | ||
+ | xlab = 'SL', | ||
+ | ylab = 'PL', | ||
+ | zlab = 'SW', | ||
+ | grid = F, | ||
+ | sphere.factor = 1)) | ||
+ | |||
+ | # Example 3 | ||
+ | open3d() | ||
+ | rgl.bringtotop(stay = T) | ||
+ | with(iris, plotlm3d(Sepal.Length,Sepal.Width, Petal.Length, | ||
+ | model = c('z ~ x + y', | ||
+ | 'z ~ x + y + I(x^2) + I(y^2) + I(x*y)'), | ||
+ | surface = T, | ||
+ | xlab = 'SL', | ||
+ | ylab = 'PL', | ||
+ | zlab = 'SW', | ||
+ | grid = F, | ||
+ | sphere.factor = 1)) | ||
+ | |||
+ | # Example 4 | ||
+ | open3d() | ||
+ | rgl.bringtotop(stay = T) | ||
+ | with(iris, plotlm3d(Sepal.Length, Sepal.Width, Petal.Length, | ||
+ | model = c('z ~ x + y', # to setosa | ||
+ | 'z ~ x + y + I(x^2) + I(y^2) + I(x*y)', # to versicolor | ||
+ | 'z ~ I(x^3) + I(y^3)'), # to virginica | ||
+ | groups = Species, | ||
+ | model.by.group = T, | ||
+ | simple.axes = F, | ||
+ | box = F, | ||
+ | xlab = 'SL', | ||
+ | ylab = 'PL', | ||
+ | zlab = 'SW', | ||
+ | grid = F, | ||
+ | sphere.factor = 1)) | ||
+ | |||
+ | # Example 5: Netter | ||
+ | x = c( 274, 180, 375, 205, 86, 265, 98, 330, 195, 53, | ||
+ | 430, 372, 236, 157, 370) | ||
+ | y = c(2450, 3254, 3802, 2838, 2347, 3782, 3008, 2450, 2137, 2560, | ||
+ | 4020, 4427, 2660, 2088, 2605) | ||
+ | z = c( 162, 120, 223, 131, 67, 169, 81, 192, 116, 55, | ||
+ | 252, 232, 144, 103, 212) | ||
+ | |||
+ | mreg = lm(z ~ x + y) | ||
+ | ndata = data.frame(x = c(150, 274, 220, 370), y = c(4000, 2800, 3500, 3100)) | ||
+ | zpred = predict(mreg, newdata = ndata, se.fit = F) | ||
+ | |||
+ | open3d() | ||
+ | rgl.bringtotop(stay = T) | ||
+ | plotlm3d(x, y, z, | ||
+ | surface = T, | ||
+ | model = 'z ~ x + y', | ||
+ | xlab = 'x', | ||
+ | ylab = 'y', | ||
+ | zlab = 'z') | ||
+ | spheres3d(x = c(150, 274, 220, 370), y = c(4000, 2800, 3500, 3100), zpred, | ||
+ | col = 'red', radius = 60) | ||
+ | </code> | ||
+ |