Data Preparation – Part II

Published by:

This time i will talk about how to deal with large text files in chuncks with R. Just to provide some real data to work with download Airlines data, relative to 1988; from now on i will work with this file.

To work with this data i will use  iterators package. This package allow you pass the file, line by line, or chunck by chunk, without really load all file to memory. As you can feel the idea try this code:

 

OK, now you have a connection to your file. Let’s create a iterator:

As you can see you are printing line by line. So, you can work with one line, or a chunk of data even with a large file. If you want to read line by line till the end of the file you can use something like this:

that returns a FALSE at the end of the file. This a very useful trick in data preparation with large text files.

Share Button

MOOCs and courses to learn R

Published by:

Inspired by this article i thought about gather here all multimedia resources that i know to learn use R. Today there is a lot of online courses, some MOOC’s too, that offer reasonable resources to start with R.

I will just list the materials in sequence and offer my evaluation about them. Of course your evaluation can be different; this case fell free to comment. In the future i can update the material. Let’s begin:

1. Computing for Data Analysis

This course was offered multiple times from 2012 to 2014. I did the course at 2013 and the course was very well formatted, with good exercises and lots of resources. It was offered on Coursera platform, that i particularly think excellent, and take an average effort to finish. There is quiz questions and program assignments. The course is free.

2. R Programming

This course is offered at coursera too, from the same instructor as Computing for Data Analysis. According to course syllabus: “The course covers practical issues in statistical computing which includes programming in R, reading data into R, accessing R packages, writing R functions, debugging, and organizing and commenting R code. Topics in statistical data analysis and optimization will provide working examples.”

Well, once that there is now the coursera specializations, i thinks that this course is “Computing for Data Analysis” rebuilded. SO, if you already did the first course i don’t see any advantage doing this one, unless you want to get the specialization certificate.

3. Data Analysis

This course was about the practical aspect of data analysis with R. All activities was in R, but the course wasn’t about R itself. But it was very good and you can learn a lot of R with it. The course is free and offered on Cousera platform.

OBS: Both, Computing for Data Analysis and Data Analysis, can be replaced for other courses on coursera specialization. I won’t comment now about these new courses because some of them are being offered now and other will be offered on the  next months. Anyway both course materials are avaliable on youtube.

4. Introduction to Data Analysis with R

This course is being offered at bigdatauniversity platform. The course is a good starting point but i don’t this platform like so much. But the pros of this course is that you can do it on your pace and it have less material to cover. It’s a quick introduction  to R.

5. Data Mining with R

This course will be too offered at bigdatauniversity, but it’s not avaliable right now. It will be recorded, and you can participate through live stream from bigdatauniversity, and released on that platform. The course doesn’t have a syllabus available, but i think that the content of this book, will probably be the topics discussed.

Rattle is a GUI for datamining that uses R as backend. It’s very intuitive and resembles Weka interface. While it’s not as flexible as use R directly it provides a quick way to explore and buid models with R. With Rattle() every step taken is saved on a log that you can use as scripts to automate tasks.

6. Using R with Databases

This course is too being offered at bigdatauniversity. It’s supposing you have some R skills and is about the use of databases with R. It uses IBM DB2 specifically, but you can  apply the concepts to others databases as well. It’s free and i liked it.

7. Statistical Learning

I’m taking this course right now and i will write a post specifically about it. BUT, just to clarify, it’s a course about machine learning (or statiscical learning if you want) based on the introductory book Introduction to Statistical Learning

The course isn’t about R itself but all the techniques are implemented using R. This course is a easy way to read the book. The course is free and is offered at openedX platform, that is a wonderful platform for MOOCs.


8. Introduction to R

This course is offered at Udemy platform. Udemy is a great platform both to create your own course as to take a course. The course is not free, but it has good reviews of users and you have 30 days to evaluate the material and be refounded.

The author claim that you can learn the basics with 2 days of course.

9. Explore Statistics with R

This course will be released at september on edX. So no further comments.

10. Build Web Apps in R with Shiny

This is another Udemy course. Worth of checking.

Obs: Coursera will offer Developing Data Products for free next june. This course is about the very same topics. So maybe it’s worth wait.

Share Button

Genetic data, large matrices and glmnet()

Published by:


Recently talking to a colleague, had contact with a problem that I had never worked with before: modeling with genetic data. I have no special knowledge of the subject, but taking a look at some articles in the area knew that one of the most used techniques for this type of data was the lasso.

In R, one of the most used packages for the lasso is glmnet, which unlike most other packages like lm accepts as input a data.frame. So, before you start modeling, you must perform a pre-processing step passing the data to matrix format. Done it is possible to pass a formula or even passing an array with the response variable, plus another with data for the variables.

The problem with the formula approach is that, in general, genomic data has more columns than observations. The data that I worked in that case had 40,000 columns and only 73 observations. In order to create a small set of test data, run the following code:

So, with this data set we will try to fit a model with glmnet ():

And if you do not have a computer with more RAM than mine, you will probably leak memory and give a crash in R. The solution? My first idea was to try sparse.model.matrix() that creates a sparse matrix model using the same formula. Unfortunately did not work, because even with sparse matrix, the final model is still too big! Interestingly, this dataset occupies only 24MB from RAM, but when you use the model.matrix the result is an array with more than 1Gb.

The solution I found was to build the matrix on hand. To do this we encode the array with dummy variables, column by column, and store the result in a sparse matrix. Then we will use this matrix as input to the model and see if it will not leak memory:

NOTE: Pay attention to how we are using a sparse matrix the Matrix package is required. Also note that the columns are connected using cBind () instead of cbind ().

The matrix thus generated was much lower: less than 70 Mb when I tested. Fortunately glmnet () supports a sparse matrix and you can run the model:

So you can create models with this type of data without blowing the memory and without use R packages for large datasets like bigmemory and ff.

Share Button

Data Preparation – Part I

Published by:

mpg


The R language provides tools for modeling and visualization, but is still an excellent tool for handling/preparing data. As C++ or python, there is some tricks that bring performance, make the code clean or both, but especially with R these choices can have a huge impact on performance and the “size” of your code. A seasoned R user can manage this effectively, but this can be a headache to a new user. SO, in this series of posts i will present some data preparation techniques that anyone should know about, at least the ones i know!

1. Using apply, lappy, tapply

Sometimes the apply’s can make your code faster, sometimes just cleaner. BUT the fact is that, at least in R, is recommended avoid for loops. So, instead of using loops, you can iterate over matrixes, lists and vectors using these functions. As an example see this code:

matriz <- matrix(round(runif(9,1,10),0),nrow=3)
apply(matriz, 1, sum) ## sum by row
apply(matriz, 2, sum) ## sum by column

Particularly in this example there is no gain on performance, but you get a cleaner code.

Talking about means, sometimes tapply can be very usefull in this regard. Let’s say you want to get means by group, you can have this with one line too. For example, considering the mtcars dataset:

so

tapply(mtcars$hp, mtcars$cyl, mean)

and you can have the mean power by cylinder capacity. This function is very usefull on descriptive analysis. BUT sometimes you have lists, not vectors. In this case just use lappy or sapply (simplify the output). Let’s generate some data:

lista <- list(a=c('one', 'tow', 'three'), b=c(1,2,3), c=c(12, 'a')) 

Each element of this list is a vector. Let’s say you want to know how many elements there is in each vector:

lapply(lista, length) ## return a list
sapply(lista, length) ## coerce to a vector

2. Split, apply and recombine

This technique you must know. Basically we split the data, apply a function and combine the results. There is a package created with this in mind. But we will use just base R functions: split, *apply and cbind() ou rbind() when needed. Looking again at mtcars dataset, let’s say we want fit a model of mpg against disp, grouped by gears,  and compare the regression coefficients.

mpg

data <- split(mtcars, mtcars$gear) ## split 
fits <- lapply(data, function(x) return(lm(x$mpg~x$disp)$coef)) ## apply
do.call(rbind, fits) ## recombine

This technique is powerfull. You can use at different contexts.

Next part i will talk about some tricks with dates.

Share Button

Dados genéticos, grandes matrizes e o glmnet()

Published by:


Recentemente, conversando com um colega, tive contato com um problema com o qual eu nunca tinha trabalhado antes: modelagem com dados genéticos. Não tenho nenhum conhecimento especial do assunto, mas dando uma olhada em alguns artigos da área soube que uma das técnicas mais utilizadas para esse tipo de dado era o lasso.

No R, um dos pacotes mais utilizados para  o lasso é o glmnet, que diferente da maioria dos outros pacotes como o lm() não aceita como entrada um data.frame. Assim, antes de iniciar a modelagem, é necessário realizar uma etapa de pré-processamento passando os dados para formato matricial. Feito isso é possível passar uma fórmula ou mesmo passar um vetor com a variável resposta, mais uma matriz com os dados das variáveis.

O problema com a abordagem da fórmula é que, em geral, dados de genoma tem muito mais colunas do que observações. Os dados que eu trabalhei nesse caso tinham 40 mil colunas e somente 73 observações. De forma a criar um pequeno conjunto de dados para teste, rode o código a seguir:

data <- matrix(rep(0,50*49000), nrow=50)
for(i in 1:50) {
   x = rep(letters[2:8], 7000)
   y = sample(x=1:49000, size=49000)
   data[i,] <- x[y]
}

data <- as.data.frame(data)
x = c(rep('A', 20), rep('B', 15), rep('C', 15))
y = sample(x=1:50, size=50)
class = x[y]
data <- cbind(data, class)

Assim, com esse conjunto de dados vamos tentar ajustar um modelo com o glmnet():

formula <- as.formula(class ~ .)
X = model.matrix(formula, data)
model <- cv.glmnet(X, class, standardize=FALSE, family='multinomial', alpha=1, nfolds=10)

E se você não tiver um computador com muito mais memória RAM do que o meu, provavelmente você vai estourar a memória e dar um crash no R. Qual a solução? Minha primeira ideia foi tentar o sparse.model.matrix() que cria uma matriz de modelo esparsa usando a mesma fórmula. Infelizmente também não funcionou, pois mesmo sendo esparsa a matriz de modelo final ainda é muito grande! O interessante é que esse conjunto de dados, na memória, tem somente 24Mb, mas sempre que se usa o model.matrix o resultado é uma matriz com mais de 1Gb.

A solução que eu encontrei foi fazer a matriz na unha pois nesse caso usei a codificação de 1 000…000. Para isso vamos decodificar a matriz com variáveis dummy, coluna por coluna e armazenar o resultado em uma matriz esparsa. Em seguida vamos utilizar essa matriz como entrada para o modelo e ver se não vai estourar a memória:

## Cria a matriz de entrada usando a primeira coluna
X <- sparse.model.matrix(~data[,1]-1)

## Verifica se a coluna tem mais de um nível (tipo factor naturalmente!)
for (i in 2:ncol(data)) {

## Se tiver mais de um nível aplica codificação dummy
if (nlevels(data[,i])>1) {
  coluna <- sparse.model.matrix(~data[,i]-1)
  X <- cBind(X, coluna)
}
## Se tiver um número transforma em factor e depois em numérico.
else {
  coluna <- as.numeric(as.factor(data[,i]))
  X <- cBind(X, coluna)
}

OBS: Preste atenção que como estamos utilizando matrizes esparsas o pacote Matrix é necessário. Também observe que as colunas são ligadas usando o cBind() ao invés de cbind().

A matriz gerada dessa forma ficou muito menor: menos de 70Mb quando eu testei. Felizmente o glmnet() tem suporte a matrizes esparsas e é possível rodar o modelo com:

mod.lasso <- cv.glmnet(X, class, standardize=FALSE, family='multinomial', alpha=1, nfolds=10)

Assim é possível criar modelos com esse tipo de dado sem estourar a memória e sem precisar utilizar pacotes para conjuntos grandes de dados como o bigmemory ou o ff.

Share Button

Compiling R 3.0.1 with MKL support

Published by:

mkl-large-product-plain


Before you begin, be aware that there is others excellent posts about the issue, as:

1. Compiling 64-bit R 2.10.1 with MKL in Linux

2. Speeding up R with Intel’s Math Kernel Library (MKL)

3. Performance benefits of linking R to multithreaded math libraries

4. Using MKL-Linked R in Eclipse

But, as recently i faced some problems to compile R 3.0.1 com MKL support on Debian Wheezy, i would like to share this update.

First the acronyms: what is MKL? MKL is an acronym for Math Kernel Library, i.e., a library that provides a series of mathematical functions, as matrix product or even a Fourier transform. But what  is special and why I would need a library of functions like that on R?

mkl-large-product-plain

First it is important to note that regardless of your purpose  to use the MKL, this is a library specific to Intel processors. If you use other architectures is likely that MKL and this post can be useless to you. BUT, specifically in relation to R, you may already have heard that it can be very slow. One of the origins of a slow code, especially with those new processors like the Intel Core i5 or i7, is that when you install the standard R (binaries available for download from CRAN) you are installing a series of numerical libraries that perform mathematical operations where performance is critical. These functions can be written in C, C + + or Fortran in order to obtain superior performance. HOWEVER, the numerical libraries that comes with the R (generic BLAS, LAPACK) do not have multithreading capabilities. In practice this means that if you have a processor with 8 cores, using the R default you will be performing several critical operations with a single processor!

How to solve the problem? One possible solution is to install the MKL on R. The MKL, beyond the capabilities of multithreading, also has specific optimizations for Intel processors, with performance far superior to traditional libraries. So from now on I will show you how you can download and install this library and link it to R. But first a few caveats:

1) If you are using Windows, can be more interesting for you to use Revolution Analytics R, that already comes standard with MKL installed.

2) If you use Linux Red Hat and derivatives (CentOS) the Revolution R should also be the best option.

3) If you use Ubuntu you can install R from repositories and then install the packages r-revobase-revolution, revolution-r and revolution-mkl.

and at this point you may be asking yourself: if there are these options why I would compile R? Well, there is several reasons: you might want to do a custom install; use R for commercial purposes without having to pay a license for the Revolution; learn how to compile R. That said, I’ll have to say that the official documentation of the R is a bit confusing and outdated in relation to the installation of MKL. That’s why I’m writing this post.

1. Downloading and Installing MKL

I will not dwell in that part, but the installation process is actually very simple: just enter the site and register. Then you will receive an email with download instructions. Extract the downloaded file and you will find scripts and RPM packages. If you are trying to install on Debian may be necessary to install support for RPM packages. Install using the install.sh script. From this moment I will assume that your installation is at / opt/intel/mkl/11.0.

2. Compiling R

To compile the R you will need the source code. Extract to the directory of your choice:

> tar -zxvf R-3.0.1.tar.gz

Before starting the compilation is necessary to set two system variables: MKLROOT and LD_LIBRARY_PATH, which can be done at once with the following command:

> source /opt/intel/mkl/bin/mklvars.sh intel64

note the location of the file mklvars.sh that may be in a different place depending on your installation of MKL.

Once defined system variables just run the usual Linux compilation  with . /configure && make && make install. The moment to “tell” to the compiler where the MKL is installed and allow R has access to it, is at the . /configure step. So:

> ./configure --enable-R-shlib --enable-threads=posix --with-lapack --with-blas="-fopenmp -m64 -I$MKLROOT/include -L$MKLROOT/lib/intel64 -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lpthread -lm"
> make
> make install

And from that moment your R already has support for multithreading libraries from Intel. To ensure that libraries are installed correctly, run this benchmark script  before and after installation, and you can check differences.

source(R-benchmark-25.R)

If you already have an installation of R on your computer, you may need to reinstall all the packages again. An easy way to do this is with the following command:

update.packages(checkBuilt=TRUE)

Just as a reference, in my notebook (Core i7, 8GB RAM) a product of two matrices 5000×5000 takes around 19 seconds with MKL. Prior to the installation was more than 1 minute!

Share Button

Recompilando o R 3.0.1 com suporte à MKL

Published by:


Antes de você começar a leitura gostaria de destacar que existem outros excelentes posts sobre o assunto:

1. Compiling 64-bit R 2.10.1 with MKL in Linux

2. Speeding up R with Intel’s Math Kernel Library (MKL)

3. Performance benefits of linking R to multithreaded math libraries

4. Using MKL-Linked R in Eclipse

no entanto, como eu fui tentar instalar a MKL no R 3.01 (Debian Wheezy) e tive alguns contratempos, gostaria de compartilhar essa atualização.

Primeiramente, vamos às siglas: o que é MKL? MKL é uma sigla para Math Kernel Library, isto é, uma biblioteca que oferece uma séries de funções matemáticas como as que são usadas para efetuar um produto matricial ou mesmo uma transformada de Fourier. Mas o que ela tem de especial e porque eu precisaria de um biblioteca de funções como essa no R?

Primeiramente é importante salientar que independente do objetivo que se tenha ao utilizar a MKL, esta é uma biblioteca específica para processadores da Intel. Se você utiliza outras arquiteturas é provável que a MKL e esse post sejam inúteis para você. MAS, especificamente em relação ao R, você pode já ter ouvido falar que ele pode ser muito lento. Uma das origens de um código lento, principalmente nesses processadores novos da Intel como o core i5 ou i7, é que quando você instala o R padrão (os binários disponíveis para download no CRAN) você está instalando uma séries de bibliotecas numéricas que executam operações matemáticas onde o desempenho é crítico. Essas funções podem ser escritas em C, C++ ou Fortran de forma obter um desempenho superior. NO ENTANTO, as bibliotecas numéricas que vem com o R (generic BLAS; LAPACK) não tem recursos de multithreading. Na prática isso significa que se você têm um processador com 8 cores, ao usar o R padrão você estará realizando várias operações críticas com um só processador!

Como resolver o problema? Uma solução possível é a instalação da MKL no R. A MKL, além dos recursos de multithreading, também tem otimizações específicas para processadores da Intel, tendo um desempenho muito superior às bibliotecas tradicionais. Assim, a partir de agora vou mostrar como você pode baixar e instalar essa biblioteca junto ao R no seu computador. Mas antes algumas ressalvas:

1) Se você estiver utilizando o Windows, talvez seja mais interessante para você utilizar o R da Revolution Analytics que já vem por padrão com a MKL instalada. O R da Revolution tem um licença gratuita para uso acadêmico.

2) Se você usa o Linux Red Hat e derivados (CentOS) o R da Revolution também deve ser a melhor opção.

3) Se você usa o Ubuntu você pode instalar o R dos repositórios e depois instalar os pacotes r-revolution-revobase, revolution-r e revolution-mkl.

e nesse momento você pode estar se perguntando: se existem essas opções porque eu recompilaria o R?? Bom, exitem vários motivos: você pode querer fazer uma instalação customizada; usar o R para fins comerciais sem ter que pagar um licença para a Revolution; aprender como compilar o R. Dito isto, gostaria se salientar também que a documentação oficial do R é um pouco confusa e desatualizada em relação à instalação da MKL. É por isso que estou escrevendo esse post.

1. Baixando e Instalando a MKL

Não vou me alongar nessa parte, mas o processo de instalação é realmente muito simples: basta entrar no site da MKL e se cadastrar. Em seguida você vai receber um e-mail com as instruções para o download. Extraia o arquivo de download e você vai encontrar scripts e pacotes  RPM. Assim se você estiver tentando instalar no Debian pode ser necessário o suporte ao RPM. Instale utilizando o script install.sh. A partir desse momento vou supor que sua instalação da MKL está em /opt/intel/mkl/11.0.

2. Recompilando o R

Para recompilar o R você vai precisar do código fonte. Extraia para o diretório de sua escolha com:

> tar -zxvf R-3.0.1.tar.gz

Antes de iniciar a compilação é necessário definir valores corretos para duas variáveis do sistema: MKLROOT e LD_LIBRARY_PATH, que pode ser feito de uma vez com o seguinte comando:

> source /opt/intel/mkl/bin/mklvars.sh intel64

observe a localização do arquivo mklvars.sh que pode estar em uma lugar diferente dependendo da sua instalação da MKL.

Definidas as variáveis do sistema basta executar a compilação canônica no Linux, com ./configura && make && make install. O momento de “dizer” ao compilador onde está instalada a MKL e permitir que o R tenha acesso a ela é no ./configure. Assim:

> ./configure --enable-R-shlib --enable-threads=posix --with-lapack --with-blas="-fopenmp -m64 -I$MKLROOT/include -L$MKLROOT/lib/intel64 -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lpthread -lm"
> make
> make install

E a partir desse momento o seu R já tem suporte às bibliotecas multithreading da Intel. Para garantir que as bibliotecas estão instaladas corretamente rode este script de benchmark antes e depois da instalação, e compare as diferenças.

source(R-benchmark-25.R)

Caso você já tenha uma instalação do R no seu computador, pode ser necessário reinstalar todos os pacotes novamente. Uma forma fácil de fazer isso é com o seguinte comando:

update.packages(checkBuilt=TRUE)

Só como uma referência, no meu notebook (core i7, 8Gb de RAM) um produto de duas matrizes 5000×5000 leva em torno de 19 segundos com a MKL. Antes com a instalação padrão era mais de 1 minuto!

Share Button

ANOVA and Tukey’s test on R

Published by:


OBS: This is a full translation of a portuguese version.

In many different types of experiments, with one or more treatments, one of the most widely used statistical methods is analysis of variance or simply ANOVA . The simplest ANOVA can be called “one way” or “single-classification” and involves the analysis of data sampled from more then one population or data from experiments with more than two treatments.

It’s not my intent to study in depth the ANOVA, but to show how to apply the procedure in R and apply a “post-hoc” test called Tukey’s test. When we are conducting an analysis of variance, the null hypothesis considered is that there is no difference in treatments mean, so once rejected the null hypothesis, the question is what treatment differ.

To illustrate the procedure we consider an experimental situation where a company is applying a sensory test for a set of 15 panelists in three different brands of chocolate. Three brands are compared, one being the reference, and the goal is to verify the difference of marks with the control. In this experiment we have two factors, the brand and the tasters, and we hope that there is no significant effect of tasters. At each assessment the assessor must determine the difference on a scale 0-7.

data.frame(
  Sabor =
  c(5, 7, 3,
    4, 2, 6,
    5, 3, 6,
    5, 6, 0,
    7, 4, 0,
    7, 7, 0,
    6, 6, 0,
    4, 6, 1,
    6, 4, 0,
    7, 7, 0,
    2, 4, 0,
    5, 7, 4,
    7, 5, 0,
    4, 5, 0,
    6, 6, 3
  ),
Tipo = factor(rep(c("A", "B", "C"), 15)),
Provador = factor(rep(1:15, rep(3, 15))))

The average grade for types A, B and C were respectively 5.33 5.27 and 1.53. Clearly C, the control, was not different for most of the tasters, as expected. One way to easily obtain these averages per group is using tapply.

tapply(chocolate$Sabor, chocolate$Tipo, mean)

Finally we run ANOVA and assess whether there are differences between brands and tasters.

ajuste <- lm(chocolate$Sabor ~ chocolate$Tipo + chocolate$Provador)
summary(ajuste)
anova(ajuste)

what results:

From the output we see that the p-value is 0.8175 for tasters indicating that the assessors have no significant effect on the response. This is desirable since it is expected that the testers can discern correctly the marks of  chocolate. Also in the table we see that the ANOVA p-value for the type of chocolate is highly significant, indicating the difference between the marks. So from now on we can make the Tukey test to see where the differences lie.

a1 <- aov(chocolate$Sabor ~ chocolate$Tipo + chocolate$Provador)
posthoc <- TukeyHSD(x=a1, 'chocolate$Tipo', conf.level=0.95)

that results:

This output indicates that the differences C-A and C-B are significant , while B-A is not significant. A more “easy” way to interpret this output is visualizing the confidence intervals for the mean differences.

plot(a1)

Rplot

One can see that only the confidence interval for B-A contain 0. Thus, it appears that the marks A and B do not differ among themselves, but are different from brand control.

Finally an alternative way to evaluate the differences in a way more similar to the SAS is using the package agricolae. With it we will apply the same procedure, but the output is slightly different.

library(agricolae)
HSD.test(ajuste, 'chocolate$Tipo')

where the final output indicates that both A and B belong to the same group, the ‘a’, and C is different from the other two, belongs to the group ‘b’.

Share Button

ANOVA e teste de Tukey no R

Published by:

Rplot


Em muitos tipos diferentes de experimentos, com um ou mais de um fator, um dos procedimentos estatísticos mais utilizados é a análise de variância, ou simplesmente ANOVA. O ANOVA mais simples pode ser chamado “one way” ou mesmo “single-classification” e envolve a análise de dados amostrados de mais de uma população ou dados de experimentos com mais do que dois tratamentos.

Nesse post não é o meu objetivo estudar a fundo o ANOVA, mas sim mostrar como aplicar o procedimento no R e fazer uma análise “post-hoc” chamada teste de Tukey. Quando estamos realizando uma análise de variância, a hipótese nula considerada é que não há diferença na média dos tratamentos, assim, uma vez rejeitada a hipótese nula, a questão é saber quais tratamentos se diferenciam.

Para exemplificar o procedimento vamos analisar uma situação experimental onde uma empresa está aplicando um teste sensorial para um conjunto de 15 provadores em três marcas diferentes de chocolate. Três marcas são comparadas, sendo uma delas a referência, e o objetivo é verificar se existe diferença das marcas com o controle. Nesse experimento temos dois fatores, a marca e os provadores,  e esperamos que os provadores não tenham um efeito significativo. Em cada avaliação o provador tem que determinar essa diferença em uma escala que vai de 0 a 7.

data.frame(
  Sabor =
  c(5, 7, 3,
    4, 2, 6,
    5, 3, 6,
    5, 6, 0,
    7, 4, 0,
    7, 7, 0,
    6, 6, 0,
    4, 6, 1,
    6, 4, 0,
    7, 7, 0,
    2, 4, 0,
    5, 7, 4,
    7, 5, 0,
    4, 5, 0,
    6, 6, 3
  ),
Tipo = factor(rep(c("A", "B", "C"), 15)),
Provador = factor(rep(1:15, rep(3, 15))))

A média das notas para os tipos A, B e C, foi respectivamente 5,33 5,27 e 1,53. Claramente a marca C, o controle, não foi diferente dela mesma para a maioria dos provadores, como esperado. Uma forma de obter essas médias facilmente por grupo é utilizando o tapply.

tapply(chocolate$Sabor, chocolate$Tipo, mean)

Por fim vamos rodar o ANOVA e avaliar se há diferença entre marcas e provadores.

ajuste <- lm(chocolate$Sabor ~ chocolate$Tipo + chocolate$Provador)
summary(ajuste)
anova(ajuste)

o que resulta:

Pela saída do programa vemos que o p-valor para provadores é 0.8175 indicando que os provadores não tem efeito significativo na resposta. Isso é desejável, uma vez que se espera que os provadores sejam capazes de discernir corretamente as marcas de chocolate. Também na tabela ANOVA vemos que o p-valor para o tipo de chocolate é altamente signigicativo, indicando a diferença entre as marcas. Assim, a partir de agora podemos fazer o teste de Tukey para verificar onde estão as diferenças.

a1 <- aov(chocolate$Sabor ~ chocolate$Tipo + chocolate$Provador)
posthoc <- TukeyHSD(x=a1, 'chocolate$Tipo', conf.level=0.95)

que resulta:

Essa saída indica que as diferenças C-A e C-B são significativas, ao passo que B-A não é significativa. Uma forma mais “fácil” de interpretar essa saída é visualizando os intervalos de confiança para as diferenças das médias.

plot(a1)

Rplot

Pode-se ver que somente o intervalo de confiança para B-A contêm o 0. Assim, verifica-se que as marcas A e B não se diferenciam entre si, mas se diferenciam da marca controle.

Por fim uma forma alternativa para avaliar as diferenças com uma saída mais parecida com a do SAS é utilizando o pacote agricolae do R. Com ele vamos realizar o mesmo procedimento anterior, mas a saída é um pouco diferente.

library(agricolae)
HSD.test(ajuste, 'chocolate$Tipo')

onde a saída final indica que tanto A quanto B pertencem ao mesmo grupo, o ‘a’, e C é diferente dos outros dois, pertence ao grupo ‘b’. Espero que esse exemplo permita que outros usuários utilizem o R para fazer esse tipo de análise.

Share Button

Plataforma edX liberada!

Published by:

edX


O edX é uma iniciativa de ensino à distância fundada por Harvard e MIT. Além da parceria com as Universidades para criar, gerenciar e oferecer gratuitamente os cursos, também foi criada uma plataforma com recursos Web 2.0 e inteligência artificial que permite a eles oferecerem cursos com milhares participantes e um staff reduzido. No oferecimento desses cursos chamados MOOC, a plataforma é uma peça fundamental.

edX

Como forma de fomentar esse modelo de educação superior e permitir o acesso a essas tecnologias em todas as partes do mundo, o edX se juntou a Stanford para construir uma plataforma de código aberto que pode ser utilizada por empresas, universidades e individuais. A partir de hoje, dia 1 de junho de 2013, a plataforma estará disponível para download. Qualquer pessoa pode instalar, criar cursos e oferece-los on-line. De acordo com as palavras do fundador do edX: “eu sempre digo que o edX vai se tornar o Linux da educação”, assim espero.

 Maiores detalhes em http://code.edx.org.

Share Button