A Data Scientist's blog

Top Machine Learning Algorithms for Data Scientists

2016-04-19T03:31:00.004+05:30

Machine Learning is a first-class ticket to the most exciting careers in data analysis today. As data sources proliferate along with the computing power to process them, going straight to the data is one of the most straightforward ways to quickly gain insights and make predictions.

Machine learning brings together computer science, mathematics and statistics to harness that predictive power. It’s a must-have skill for all aspiring data analysts and data scientists, or anyone else who wants to wrestle all that raw data into refined trends and predictions.

Have you ever thought of what is the end-to-end process of investigating data through a machine learning lens? Have you ever extracted and identified useful features that best represent your big data? Have you ever gone to a complex process of dealing with big data predictions and evaluated the performance of your machine learning algorithms?

Leading companies such as Amazon, Google and Facebook uses efficient algorithms fit for big data, indexation, attribution modeling, collaborative filtering, and recommendation engines. Here is the top and the most important machine learning algorithms for data analysts and data scientists.

Algorithm 1: Gradient Descent

The Gradient Descent algorithm is as used as the optimization algorithm at the core of so many machine learning algorithms.

Linear Algorithms

Algorithm 2: Linear Regression

Algorithm 3: Logistic Regression

Algorithm 4: Linear Discriminant Analysis

Nonlinear Algorithms

Algorithm 5: Classification and Regression Trees

Algorithm 6: Naive Bayes

Algorithm 7: K-Nearest Neighbours

Algorithm 8: Learning Vector Quantization

Algorithm 9: Support Vector Machines

Ensemble Algorithms

Algorithm 10: Bagged Decision Trees and Random Forest

Algorithm 11: Boosting and AdaBoost

Drawbacks of Some Algorithms

Naive Bayes

Variables are almost never uncorrelated

Linear Discriminant Analysis

Clusters are almost never separated by hyperplanes

Linear Regression

Numerous model assumptions - including linearity - are almost always violated in real data.

Data Science Architecture - What Data Scientist thinks when it comes to solving real world problems?

2016-02-27T03:32:00.001+05:30

When we solve real world data problems and when it comes to building something really valuable and can change the entire decision making process you have never dreamed of. We always think of a process, architecture and methodology that can help us, guide us to achieve such remarkable moments. There are many components of data science, machine learning, statistical analysis, mathematical modeling projects as well as the cross dependencies between these components that help us to solve the problems we have been dreaming about.

Here is a vary standard data science architecture that may lead us to reach beyond it.

Problem

When comes to identifying the real problem that we need to solve, we ask the business people, what the hack is happening out there? It can be anything from building a prediction models, building a market segmentation, building a recommendation engine, association rule discovery for fraud detection, minimizing production costs, minimizing advertisement costs, maximizing ROI, right rewards to right users, right offers from offline to online to right users, best deals and offers to right users, right gift to right users, or simulations to predict extreme events such as floods.

Data

Are you sure the data you gave is correct? I have been giving you incorrect data for years. This is the first time you have asked. What? I said the data is totally accurate. Think about it, dive into it. It comes in many shapes: transactional, real-time, sensor data (IoT), unstructured data (tweets), structured data, big data, images or videos, pdfs, news, press releases, public hearings, and so on. Typically raw data needs to be identified or even built and put into databases (NoSQL or traditional), then cleaned and aggregated using EDA (exploratory data analysis). The process can include selecting and defining metrics. Check outliers, check missing values, think how to treat them. Impute it? drop them? Check data sources? It totally depends on you how you want to take it.

Algorithms

Stay focused. Think what you have learned. Can you list steps how you used to go to schools? Can you write how you can cross a river? Tough eh? Ok can you write how you can cross a road? Looks easy task right? But know life is not so easy specially when it comes a data science problem. We need to make our mind sharp and start thinking some cutting edge techniques. Examples include agent based modeling, clustering, Hidden Markov model, structural equation modeling, decision trees, indexation algorithm, Bayesian networks, attribution modeling, Monte Carlo simulation rule based machine learning or support vector machines. A rather big list can be found here.

Models

Before reading next few lines, do not get scared, we are born to be trained, we are born to learn and teach. By models, it means testing algorithms, selecting, fine-tuning, and combining the best algorithms using techniques such as model fitting, model blending, data reduction, feature selection, and assessing the yield of each model, over the baseline. It also includes calibrating or normalizing data, imputation techniques for missing data, outliers processing, cross-validation, over-fitting avoidance, robustness testing and boosting, and maintenance. Criteria that make a model desirable include robustness or stability, scalability, simplicity, speed, portability, adaptability (to changes in the data), and accuracy (sometimes measured using R-squared, though I recommend this alternative instead). Take a long breathe and say all is well vision that we have still whole life to learn all of them :)

Programming

Think healthy, eat healthy, stay healthy, and program healthy. There is almost always some code involved, even if you use a black-box solution. Generally we use R, Python, Matlab, Java, SQL, Julia, Spark, Scala, Haskel, Prolog, etc. to design and implement our thoughts, models, algorithms, techniques, even steps we lists.

Environments

Some call it libraries, some call it packages. It can be anything such as a bare Unix box accessed remotely combined with scripting languages and data science libraries such as Pandas (Python), Shiny (R), Rattle (R) or something more structured such as Hadoop, Cassandra, MongoDB, Apache, or Hive. Or it can be an integrated database system from Teradata, Pivotal or other vendors, or a package like SPSS, SAS, RapidMiner or MATLAB, or typically, a combination of these.

Presentation

We may call it business intelligence dashboards (BID), Report, Application, Presentation, PPT, Deck etc. By presentation, BID, report, etc., it means presenting the results. Not all data science projects run continuously in the background, for instance to automatically buy stocks or predict the weather. Some are just ad-hoc analyses that need to be presented to decision makers, using Excel, Tableau and other tools. In some cases, the data scientist must work with business analysts to create dashboards, or to design alarm systems, with results from analysis e-mailed to selected people based on priority rules.

Installing latest version of R-base in ubuntu

2014-06-09T10:37:00.002+05:30

apt-cache showpkg r-base

which would give something like

Package: r-base
Versions:
3.0.2-1raring0 (/var/lib/apt/lists/cran.rstudio.com_bin_linux_ubuntu_raring_Packages)
Description Language:
                 File: /var/lib/apt/lists/cran.rstudio.com_bin_linux_ubuntu_raring_Packages
                  MD5: 5787ca79ed716232c4cc2087ed9b425b

3.0.1-6raring0 (/var/lib/apt/lists/cran.rstudio.com_bin_linux_ubuntu_raring_Packages)
Description Language:
                 File: /var/lib/apt/lists/cran.rstudio.com_bin_linux_ubuntu_raring_Packages
                  MD5: 5787ca79ed716232c4cc2087ed9b425b

and the I did a

sudo apt-get install -f r-base=3.0.2-1raring0

and done.

Second Method

Uninstall old R:
sudo apt-get remove r-base-core
Then:
sudo gedit /etc/apt/sources.list
Add the following to the file:
deb http://cran.rstudio.com/bin/linux/ubuntu precise/
and exit gedit.
Then copy/paste these commands into the command line:

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E084DAB9
sudo add-apt-repository ppa:marutter/rdev
sudo apt-get update
sudo apt-get upgrade
sudo apt-get install r-base

Ubuntu boot repair for windouws boot loader not working

2014-06-08T23:07:00.004+05:30

Boot repair
https://help.ubuntu.com/community/Boot-Repair

OTHERS

http://askubuntu.com/questions/450631/broadcom-b43-wifi-not-working-in-ubuntu-14-04-lts

http://itsfoss.com/things-to-do-after-installing-ubuntu-14-04/

http://www.howopensource.com/2012/10/install-broadcom-b43-legacy-wireless-driver-in-ubuntu-12-10-12-04/

http://www.howopensource.com/2012/10/install-broadcom-sta-wireless-driver-in-ubuntu-12-10-12-04/

http://en.community.dell.com/dell-blogs/direct2dell/b/direct2dell/archive/2008/10/03/linux-driver-available-for-dell-wireless-cards.aspx

sudo dpkg --configure -a
sudo apt-get install -fy
sudo apt-get purge -y --force-yes grub* shim-signed linux-signed*

Boot successfully repaired.

Please write on a paper the following URL:
http://paste.ubuntu.com/7613845/

In case you still experience boot problem, indicate this URL to:
boot.repair@gmail.com or to your favorite support forum.

You can now reboot your computer.

The boot files of [The OS now in use - Ubuntu 14.04 LTS] are far from the start of the disk. Your BIOS may not detect them. You may want to retry after creating a /boot partition (EXT4, >200MB, start of the disk). This can be performed via tools such as gParted. Then select this partition via the [Separate /boot partition:] option of [Boot Repair]. (https://help.ubuntu.com/community/BootPartition)

A broken Wubi has been detected. Please fix it this way:
https://wiki.ubuntu.com/WubiGuide#Cannot_boot_into_Ubuntu

10 reasons why you should learn R

2014-06-03T14:23:00.003+05:30

10. Can't crack that hard Sudoku problem?? Use R!

9. Want to pick a skill that will give you an early adopter advantage?? Learn R! It is the leading open source statistical and data analysis programming language, and is heating up!

8. Need to run statistical calculations in your software application?? Deploy R! It integrates with many programming languages like Java, Ruby, C++, Python

7. Looking for reusable libraries to solve a complex problem?? Get R! It has 2000+ free libraries to use in areas of finance, natural language processing, cluster analysis, optimization, prediction, high performance computing etc.

6. No Windows, No Doors - R runs on all the platforms. Just name it and you got it!! Windows PC, Mac, Linux to name a few

5. Did you know how much fun stats can be- Try R!!

4. Are you updated with the current trends?? Leading firms like NY Times, Google, Facebook, Bank of America, Pfizer, Merck are all using R, where are you??

3. Need to run your own analysis?? Need to solve an optimization problem?? Struggling with Excel or SQL in your model??..... just few statements away - Try R!!

2. Want to create a compelling chart?? Try R!

1. Want the coolest job in 2014?? Learn Statistics. It is the future. Data Scientists will be the sexy job in 2018.

Logistic Regression in structural equation modeling SEM framework

2013-11-14T17:46:00.002+05:30

Logistic Regression in structural equation modeling SEM framework

Dear lavaan users,

as far as i understand regression models are a special case of the more general structural equation models. i wonder, whether it is possible to model logistic regression in lavaan? Wouldn't one just need to transform the left-hand-side of the formula with the logit? Probably it's much more complicated, isn't it?

It is. Lavaan can not handle logistic regression (yet). It can handle
probit regression (using the probit link instead of the logit link) for
binary outcomes. But parameters can only be estimated by using WLS(MV),
not ML (at least not in 0.5-9), but the results are typically almost
identical. The only downside is that you can not directly interpret the
(exponentiated) regression coefficients as odds (ratios).

This would be an example (if y is binary):

model <- br="" x1="" x2="" x3="" y="">fit <- data="myData," model="" nbsp="" ordered="y" sem="" span="">

Beyond logistic regression: structural equations modelling for binary variables and its application to investigating unobserved confounders

2013-11-14T17:44:00.003+05:30

Beyond logistic regression: structural equations modelling for binary variables and its application to investigating unobserved confounders

Emil Kupek

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Structural equation modelling (SEM) has been increasingly used in medical statistics for solving a system of related regression equations. However, a great obstacle for its wider use has been its difficulty in handling categorical variables within the framework of generalised linear models.

Methods

A large data set with a known structure among two related outcomes and three independent variables was generated to investigate the use of Yule's transformation of odds ratio (OR) into Q-metric by (OR-1)/(OR+1) to approximate Pearson's correlation coefficients between binary variables whose covariance structure can be further analysed by SEM. Percent of correctly classified events and non-events was compared with the classification obtained by logistic regression. The performance of SEM based on Q-metric was also checked on a small (N = 100) random sample of the data generated and on a real data set.

Results

SEM successfully recovered the generated model structure. SEM of real data suggested a significant influence of a latent confounding variable which would have not been detectable by standard logistic regression. SEM classification performance was broadly similar to that of the logistic regression.

Conclusion

The analysis of binary data can be greatly enhanced by Yule's transformation of odds ratios into estimated correlation matrix that can be further analysed by SEM. The interpretation of results is aided by expressing them as odds ratios which are the most frequently used measure of effect in medical statistics.

Background

Statistical problems that require going beyond standard logistic regression

Although logistic regression has become the cornerstone of modelling categorical outcomes in medical statistics, separate regression analysis for each outcome of interest is hardly challenged as a pragmatic approach even in the situations when the outcomes are naturally related. This is common in process evaluation where the same variable can be an outcome at one point in time and a predictor of another outcome in future. For example, preterm delivery is both an important obstetric outcome and a risk factor for low birthweight, which in turn can adversely affect future health. Sequential nature of these outcomes is not encompassed by repeated measures models which deal with the same outcome at different time points. Another example of a research problem difficult to handle by logistic regression model is when an outcome is determined not only by direct influences of the predictor variables but also by their unobserved common cause. For example, survival time since the onset of an immune system disease may be adversely affected by concomitant occurrence of various markers of disease progression indicating immunosupression as an underlying common factor, the latter being an unobserved latent variable whose estimation requires solving a system of related regression equations.

Structural equation modelling (SEM) is a very general statistical framework for dealing with above issues. In recent years, it has been increasingly used in medical statistics. In addition to traditional areas such as psychometric properties of health questionnaires and tests, behavioural genetics [1], measurement errors [2] and covariance structure in mixed regression models [3] have received particular attention. In addition to specific applications, important research methodology issues in SEM have been given more space in medical statistics, among which a comparison with multiple regression [4], the relevance of latent variable means in clinical trials [5] and power of statistical tests [6] deserve special attention.

However, a great obstacle for wider use of SEM has been its difficulty in handling categorical variables. The aim of this paper is to briefly review main aspects of this difficulty and to demonstrate a new approach to this problem based on a simple transformation. Two examples with both simulated and real data are provided to illustrate this approach.

SEM includes both observed and unobserved (latent) variables such as common factors and measurement errors. The Linear Structural Relationships (LISREL) model [7] was the first to spread in psychometric applications due to the availability of software. Other formulations of SEM and corresponding software emerged (see [8] for an overview). The details of these models, as well as important issues regarding their identifiability, estimation and robustness, are beyond the scope of this work but an illustration of the situations where SEM is needed is presented instead (Figure 1). As a general rule, SEM is indicated when more than one regression equation is necessary for statistical modelling of the phenomena under investigation.

Figure 1. Statistical problems needing SEM approach.

The left part of Figure 1 shows a situation where two outcomes, denoted Y₁and Y₂, are mutually related (a feed-back loop) and influenced by two predictors, denoted X₁and X₂. For example, the outcomes could be demand and supply of a particular health service or risk perception and incidence of a particular health problem. The predictor variables' error terms, denoted e₁and e₂, may be correlated (r) if an important variable influencing both predictors is omitted, i.e. in the case of bias in exposure measures. The terms d₁and d₂indicate disturbances of the two outcomes. The right part of Figure 1 illustrates a combination of common factors and regression model. In this case, it is of interest to test whether the outcome Y is determined not only by direct influences of the predictor variables, denoted X₁, X₂, X₃and X₄, but also by their latent determinant as indicated by the regression coefficient b.

SEM has received many criticisms, most of which have been concerned with vulnerability of complex models relying on many assumptions, as well as with uncritical use and interpretation of SEM. These are well placed concerns but are not intrinsic to SEM; even well known and widely applied techniques such as regression share the same concerns. Complex phenomena require complex models whose inferential aspects are more prone to error as the number of parameters increases. SEM is often the only statistical framework by which many of these issues can be addressed by testing and comparing the models obtained [9].

Handling categorical variables in SEM

Specific criticism regarding the treatment of categorical and ordinal variables in SEM has been a strong deterrent for its wider use. Naive treatment of binary and ordered categorical variables as if they were normally distributed in some SEM applications was partly due to the lack of viable alternatives in its early days. Inadequate use of standardized regression coefficients as the measures of effect in some SEM applications was also criticised [10]. Even when distributional properties of categorical variables were taken into account, the interpretation of SEM parameter estimates in terms of impact measures such as attributable risk was not applied. Standard errors and confidence limits – rarely used in SEM – are generally underestimating structural model uncertainties such as selection of relevant variables and correct specification of their influences.

A recent review of handling categorical and other non-normal variables in SEM [11] listed four main strategies: a) asymptotic distribution free (ADF) estimators adjusting for non-normality by taking into account kurtosis in joint multivariate distribution [12], b) the use of robust maximum likelihood estimation or resampling techniques such as jacknife or bootstrap to obtain the standard errors of SEM parameters as these are most affected by departure from multivariate normality [13], c) calculating polyserial, tetrachoric or polychoric correlations for pairs of variables with non-normal joint distribution by assuming that these have an underlying (latent) continuous scale whose large sample joint distribution is bivariate normal, then using these correlations as the input for SEM [14], and d) estimating probit or logit model scores for observed categorical variables as the first level, then proceeding with SEM based on these scores as the second-level [15]. The ADF estimation generally requires large samples to keep the type II error at a reasonable level and extremely non-normal variables such as binary may be difficult to handle with sufficient precision. The last two strategies critically depend on how well the first-level model fits the data.

A review of statistical models for categorical data reveals the lack of a method capable of handling more than one regression equation [16]. Although log-linear models for contingency tables may analyse related categorical outcomes and their relationship with independent variables, possibly complex interactions between the variables in the model do not indicate the direction of influences as in regression models. This underlines the need for a SEM framework for categorical data analysis in order to handle both dimensionality reduction and regression techniques within the same model (cf. the right part of Figure 1).

Two major recent developments in handling categorical data include Muthen's extension of SEM to the 'latent variable modeling' approach [17] and an extension of generalized linear models to latent and mixed variables under GLLAMM (Generalized Linear Latent And Mixed Models) framework [18]. Despite coming from different statistical backgrounds, both Muthén's Mplus software [19] and GLLAMM are capable of modelling a mixture of continuous, ordinal and nominal scale variables, multiple groups (including clusters) and hierarchical (multi-level) data, random effects, missing data, latent variables (including latent classes and latent growth models) and discrete-time survival models. Both of these developments are based on the vision of generalized linear models as a unifying framework for both continuous and categorical variables, where the latter are first transformed into continuous linear functions and subsequently modelled by SEM. This paper follows the same line but proposes a different transformation for categorical variables, so far unused in SEM. A simulated and a real data example with a latent confounding variable are presented.

Methods

Data generation and transformation

This work illustrates the application of SEM for binary variables using Yule's transformation to approximate the matrix of Pearson's correlation coefficients from odds ratio (OR) by a well known formula (OR-1)/(OR+1). The first example is based on known data generating processes to avoid uncertainty about true model, virtually inevitable for empirical data. A data set with 5000 observations was generated to allow normal theory approximation. First, three continuous random variables, denominated x₁to x₃, were created from the uniform distribution. The variables were uncorrelated in the population. Their binary versions, denominated BIN₁to BIN₃, were obtained by coding the values above the mean as one versus zero otherwise. Two continuous dependent variables were created by the following equations: m = 1.5 x₁+ 2 x₂+ e₁and y = 0.5 x₂- 2.5 x₃+ 1.3 m + e₂, with e₁and e₂being normally distributed random errors (N~0,1), generated from different seeds. The binary versions of the dependent variables, denominated MBIN and YBIN, were created by applying the logistic regression classification rule, i.e. score 1 if exp(m)/(1+exp(m)) and exp(y)/(1+exp(y)) exceed 0.5 versus 0 otherwise, where 'exp' stands for 'exponentiation'.

Observed odds ratios between the variables of interest in the generated data sets are reported in table 1. The structural relationships among the variables in the second data set are depicted in Figure 2.

Table 1. Simulated data: Observed odds ratios (OR), associated 95% confidence intervals (CI) and SEM regression coefficients with corresponding standard errors (SE) obtained via ML estimation (N = 5000)

Figure 2. Simulated model.

In addition, a random sample of 100 observations was taken from the generated data set with 5000 observations in order to illustrate small sample performance of the SEM based on Yule's transformation compared to logistic regression. Finally, a real data example with related binary obstetric outcomes, including premature birth, lower segment Caesarian section, low birthweight (<2500 10574="" 1="" a="" and="" applied="" as="" baby="" between="" care="" compare="" data.="" data="" delivered="" extracted="" from="" g="" logistic="" multiparous="" obstetric="" of="" pregnancies="" records="" regression="" sem="" singleton="" special="" standard="" sup="" technique="" the="" this="" to="" type="" unit="" used="" utilization="" was="" were="" who="" with="" women="">st

August 1994 and 31^stJuly 1995 in nine maternity units in England and Wales [20].

Yule's transformation was used to estimate the matrix of Pearson's correlation coefficients for both simulated and real obstetric data. The correlations were used as input for SEM. For the simulated data, both logistic and SEM analysis were repeated for a random subset of 100 observations taken from the original data set. Maximum likelihood (ML) estimation was used.

SEM raw regression coefficients were back-transformed from Q-metric into odds metric by (1+Q)/(1-Q) to get an impact measure for the binary predictor variables. SAS software procedures CALIS and LOGISTIC were used for SEM and logistic analysis, respectively [21].

Evaluation of classification performance

Raw data residuals were calculated as the difference between observed and SEM-predicted values for both data sets. The predicted values were calculated by multiplying the raw regression parameters obtained in SEM with corresponding observed values of the predictor variables. The back-transformation from SEM parameters, denoted S, to the odds metric is given by (1+S)/(1-S) and provides the odds of being the case for each independent variable; summing these odds over the independent variables gives the odds of being the case for each profile of independent variables. The odds greater than one were classified as SEM predicted cases versus otherwise.

For logistic regression, the percent of correctly classified outcomes was calculated using the cut-off point of 0.5 for the estimated probability of outcome variables.

The classification performance of SEM and logistic regression was compared on a real data set with several obstetric outcomes of interest [20] and on a small random sample of 100 observations taken from the simulated data set of 5000 observations.

Power analysis

Statistical power analysis used a calculation based on non-central chi-squared distribution, providing the number of observations required to achieve the 90% power (beta or type II error of 0.10), denoted as N [22,23]. If n denotes the number of observations used in SEM, k denotes the multiplying factor for a chosen power level, degrees of freedom and alpha (type I error), and d denotes the chi-square difference between the SEM with and without the parameter(s) of interest, then N = k*n/d gives the required sample size. Releasing one parameter at a time (one degree of freedom), with fixed type I error of 5% and type II error of 10%, point to the tabulated k-value of 10.51 [23]. This approach assumes that the model is correctly specified.

Results

Table 1 contains observed odds ratios for the simulated data set and their decomposition into regression effects based on SEM using Yule's transformation of odds ratios.

A standard approach to the analysis of binary variables using multivariate logistic regression for the simulated data is presented in Table 2.

Table 2. Multivariate logistic regression for generated data: parameter estimates (standard errors) for large (N = 5000) and small (N = 100) samples

The normal probability plot of raw data residuals between observed outcomes and the estimated probability of outcome based on SEM for simulated data showed some departure from the normal distribution (Figure 3). On the other side, the residuals fall within the normal range. Both SEM and logistic regression models for real obstetric data (Figure 4) showed satisfactory fit regarding individual data residuals.

Figure 3. Normal probability plots for raw data residuals. Normal probability plots for raw data residuals in the simulated data model with two related outcomes: YBIN (top) and MBIN (bottom). Asterisk may represent up to 30 residuals.

Figure 4. Comparison of SEM and logistic model estimates for the obstetric data example.

The comparison of classification performance for SEM versus logistic regression showed slightly better results with the latter for one outcome in a small sample analysis and very similar results for all other comparisons (Table 4). True positive fraction for events was always considerably higher for SEM compared to logistic regression, albeit at the expense of lower true negative fraction for non-events.

Table 3. Small sample (N = 100) parameter estimates and their standard errors (SE) for SEM using Q-statistic input (correlations estimated via Yule's transformation)

Table 4. Percentage of correctly classified events for logistic regression (LR) models in table 2 versus SEM in tables 1 and 3

Logistic regression showed better overall classification rate due to better prediction of non-events (Table 5). On the other hand, events were better predicted by SEM.

Table 5. Classification performance for the obstetric data example (N = 10574): logistic regression (LR) and SEM with Q-metric input (see Figure 4)

SEM permitted further investigation of the unobserved determinant of observed obstetric risk factors in predicting the need for specialised neonatal care through a latent variable. A model was tested assuming that a common cause of some of the risk factors is a latent confounding variable influencing both observed risk factors and the outcome of interest (special baby care unit) and adding predictive power over and above the observed risk variables (Figure 5). The estimation was possible upon solving the observed variables' parameters first (so-called path analysis) and fixing the factor loading for preterm delivery to the value of one – a convention allowing the comparison of the contribution of the other two observed risk variables to the unobserved latent risk using premature birth as unit risk. The factor loadings (standard errors) for Caesarian section and low birthweight were -0.3948 (0.003) and 0.8630 (0.001), respectively.

Figure 5. SEM with latent risk variable for the obstetric data example.

The relevance of the latent variable for predicting the use of special care baby unit was also tested by linear regression with raw data SEM residuals (observed minus SEM predicted probability of using special care baby unit) as the dependent variable and the latent variable scores as the predictor variable. The predictor was estimated at 0.0874 (standard error 0.0053) and was highly significant (p < 0.001).

The model suggested that propensity for premature birth resulting in low birthweight upon delivery which did not use Caesarian section increased the chances of special neonatal care utilization. The raw SEM coefficient representing this effect, denominated b₄on Figure 5, was estimated at 0.0956 with corresponding standard error of 0.016, leading to a highly significant t-value of 61.54. Transforming back to odds metric via (1+b₄)/(1-b₄) resulted in odds ratio of 1.21 and corresponding 95% confidence intervals from 1.14 to 1.29. Although a multivariate logistic regression model for the special baby care unit utilization did not find the above combination of risk factors statistically significant when it was added as interaction term to the risk factors themselves (odds ratio 1.16 with 95% confidence intervals from 0.72 to 1.86), it should be stressed that this is a model different from the above SEM.

Statistical power analysis found that only the b₃parameter in table 3 would require a larger sample size (N = 5918) than the one available to achieve the 90% power.

Discussion

The analysis demonstrated the viability of SEM using Yule's Q-transformation of odds ratio as input for binary variables models. On the level of individual data points, the raw data residuals were within the normal range and the discriminant rule for classification of outcomes into events and non-events based on SEM Q-scores performed slightly worse but still similarly to the results based on standard approach using logistic regression. The conclusion holds for the small sample example with generated data and for the real data set tested here. All these elements point out to the feasibility and utility of SEM using Yule's transformation for binary data, principally when complex relationships between the variables are present. For example, the investigation of the common cause of obstetric risk indicators on the outcome of interest identified a latent confounding variable which increased the chances of utilizing special neonatal care over and above the impact of the same risk indicators taken as independent predictors (Figure 5). The interpretation of the latent variable may lead to hypothesising a health service routine of treating premature births in a particular way (i.e. restraining from Caesarian section) or a biological propensity for birth complications, with both of these alternatives leading to an increased need for intensive neonatal care. This illustrates how SEM helps generating and investigating complex hypothesis not available by other methods. Yule's transformation may be helpful in preparing binary data for SEM. By using odds ratio both as a starting point and for the results presentation, the proposed transformation facilitates the interpretation of effects in the model.

For alpha level <0 .05="" and="" both="" em="" for="" likelihood="" ratio="" t-test="" test="" the="" univariate="">b

₃parameter being equal to zero indicated its statistical significance in SEM (details not shown) despite non-significance of observed odds ratio (table 3). However, the power of this test is less than the pre-established criterion of 90% and the impact of this parameter is clearly inferior to that of the other predictors in the model. The tendency to include extra parameters was also reported for SEM ML estimates where ordered categorical variables were treated as continuous [24] and may be expected for ADF estimates in SEM with raw binary data input. It should be noted that binary variables and the amount of noise introduced in the model analysed are serious obstacles to specifying the correct relationship between the variables for ADF estimation methods, typically applied to the data with smaller departure from the multivariate normal distribution. However, there has been some progress in developing both large sample and finite sample robustness of SEM parameters in handling non-normal data and outliers [25,26].

The advantage of SEM over separate logistic regression models for each outcome is twofold. First, SEM can model all regression equations simultaneously, thus providing a flexible framework for testing a range of possible relationships between the variables in the model, including mediating effects and possible latent confounding variables. Second, on a more general level, SEM parameters can quantify the contribution of each predictor to the covariance structure such as common factors model (Figure 5 is an example), whereas neither the interaction of continuous variables, defined as their crossproduct, nor the interaction terms for categorical independent variables in a regression model, can do this. The modelling of a common cause of observed risk factors and its influence on the outcome of interest is impossible outside SEM framework. Genetic propensity for various diseases is probably the most vivid example of the need for above model, enabling an investigation of the latent confounding variables frequently cited in the study design literature. This includes latent growth models with a relatively long sequence of indicators of an evolving process such as disease whose symptoms are typically binary indicators used for statistical modelling of the outcomes of interest. It is no coincidence that some recent developments in regression modelling have been marked by the efforts to integrate regression with a variety of covariance structure models [1-3].

Another advantage of SEM using Yule's Q-transformation of odds ratios for binary variables over two-level approach, based on probit or logit model or estimated correlations for non-normal variables as first level and SEM as second level modeling, may lay in the fact that the former is based on data transformation rather than estimation, thus avoiding the sources of error due to the latter. However, this view is not universally accepted and the discussion goes back to the beginning of the 20^thcentury when Karl Pearson and George Udney Yule argued whether a measure of association of two binary variables needs to assume underlying continuum and bivariate normal distribution [16]. While the former based his calculation of tetrachoric correlation on these assumptions, the latter disagreed, saying that some categorical variables are inherently discrete, so that the continuum assumption is tenuous and in fact unnecessary because a measure of association for such cases can be obtained directly from cell counts in a 2 by 2 table as in odds ratio and its transformation, today known as Yule's Q. Although the popularity of odds ratio over Pearson's correlation in medical statistic points to a prevailing tendency of embracing Yule's view in this field, an attempt to reconcile the two viewpoints has been made [16].

The fact that Yule's transformation is well known and allows an easy back-transformation of model parameters to odds metric makes it easier to interpret them as effect measures. Although SEM estimates based on already existing methods for handling categorical variables could be converted to an odds ratio metric for the purpose of interpretation, it has been used very rarely in the publications in the field and almost exclusively with GLLAMM.

Usual tools for evaluating SEM fit such as the analysis of residuals are available not only for input covariance matrix but also for individual data points. When classification of outcomes into events and non-events is of interest, sensitivity and specificity parameters can easily be obtained, thus making this approach applicable to a wide range of research problems.

Although other measures of comparative model fit, abundant in the SEM literature [9], may also be useful to assess various aspects of this important issue, classification performance is a preferred measure of predictive power in practice, particularly if cross-validated. For example, both data sets analysed here used saturated models which perfectly predicted the input correlation matrices, so the fit indices based on the discrepancy between observed and SEM-predicted correlation matrices obtained maximum values possible, but this was not particularly informative. On the other hand, SEM fit indices may be useful to select the best model in many other situations.

Despite the advantages of SEM mentioned above, there are several limitations of this work. First, Yule's Q is not exactly Pearson's correlation coefficient but rather an approximation to it which seems reasonable in large samples and for the types of models tested. Although the illustration of a small sample size performance seems satisfactory compared to logistic regression models, it is yet to be tested fully for a much wider range of dependency structures than presented here in order to evaluate the robustness of the parameters obtained. However, this requirement is a consequence of complex modelling issues which often arise in SEM as Yule's Q is no new estimator. Therefore, the findings about the properties of ML, ADF and least squares estimators in SEM, accumulated for almost three decades of research, apply here. This is the main reason why no attempt of a simulation study of SEM parameter estimates has been made in this work. Second, the lack of a simple rule for variable selection in SEM and the need to test a variety of models before selecting the acceptable ones can make it difficult to use this approach for quick decision making often favoured in routine applications of medical statistics. Model selection based on Bayes factors [27] may be helpful in this situation. Finally, although logit is the most popular transformation in modelling binary outcomes in medical statistics, there are many other link functions which may be more suitable for a particular model. GLLAMM [18] theory and software seem to be the most complete framework for such investigation up to date.

When the scale of SEM variables is not equal or their variances differ significantly, covariance matrix input should be preferred instead of correlation matrix input. Although SEM standard errors are less accurate with the latter even with the sample size of few hundreds, the data used here had much larger sample sizes and therefore are less influenced by the type of input matrix. In addition, the input of all SEM variables was on the same scale, i.e. in the odds metric. On the other hand, many SEM applications are performed on moderate and small samples, so the covariance matrix input would be preferable. With multivariate normal distribution, sample covariance matrix contains all the necessary information for SEM. However, with non-normal data, kurtosis was shown to be the most relevant parameter to be taken into account to correct the standard errors of SEM parameters, as in ADF estimators [12]. If means are of interest in SEM, input covariance matrix can be augmented with this information as well. Another way of dealing with SEM standard errors from non-normal data is bootstrapping, already included in several statistical packages with SEM module.

If the raw regression parameters from SEM exceed the domain of the inverse of Yule's transformation function, i.e. the interval from -1 to 1, then standardized SEM parameters can be used to get the odds metric via (1+Q)/(1-Q). Alternatively, a transformation mapping the raw SEM coefficients to this interval may be used, such as Yule's or logit, with corresponding back-transformation of the results to odds metric.

Although this work does not address the question of the association between continuous and dichotomous variables, extensions to include this case can be envisaged. One strategy would be to transform continuous variables into ordered categories with one of them serving as a baseline and then calculate odds ratios using logistic regression. Subsequently, Yule's transformation can be used to convert the odds into correlation metric to be analyzed by SEM. Another strategy would be to use polychoric or poliserial correlation for above situation and only substitute tetrachoric correlation by Yule's Q, particularly when the structural relationships of interest are between binary variables in the model and some exogenous variables are ordered or continuous.

Further research is needed to elucidate various aspects of the SEM based on Q-metric input, particularly small sample performance for a wide range of statistical models and their classification performance. In addition, the variance of odds ratios may be used to weight the estimated correlation matrix, so that Q-metric input for SEM takes into account the precision of the original scale and not only the magnitude of association between two binary variables. Relative fit measures such as those recently proposed by Agresti & Caffo [28] may help selecting among competing models of different kind.

Conclusion

SEM based on Q-transformation of odds ratios can be used to investigate complex dependency structures such as latent confounding factors and their influences on both observed risk factors and categorical outcome variables.

Competing interests

The author(s) declare that they have no competing interest.

References

Posthum D, de Geus EJ, Neale MC, Hulshoff Pol HE, Baare WEC, Kahn RS, Boomsma D: Multivariate genetic analysis of brain structure in an extended twin design.
Behav Genet 2000, 30:311-319. PubMed Abstract | Publisher Full Text
Plummer MT, Clayton DG: Measurement error in dietary assessment, an investigation using covariance structure models.
Stat Med 1993, 12:925-948. PubMed Abstract
Littell RC, Pendergast J, Natarajan R: Modelling covariance structure in the analysis of repeated measures data.
Stat Med 2000, 19:1793-1819. PubMed Abstract | Publisher Full Text
Wall MM, Li R: Comparison of multiple regression to two latent variable techniques for estimation and prediction.
Stat Med 2003, 22:3671-3685. PubMed Abstract | Publisher Full Text
Donaldson GW: General liner contrasts on latent variable means: structural equation hypothesis tests for multivariate clinical trials.
Stat Med 2003, 22:2893-2917. PubMed Abstract | Publisher Full Text
Miles J: A framework for power analysis using structural equation modelling procedure.
BMC Med Res Methodol 2003, 3:27. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text
Joreskog K, Sorbom D: LISREL VI. Uppsala: Uppsala University Press; 1981.
Structural Equation Modelling [http://www.gsu.edu/~mkteer] webcite
Tanaka JS: Multifaceted conceptions of fit in structural equation models. In Testing structural equation models. Edited by Bollen KA, Long J. Newbury Park, CA: Sage Publications; 1993::10-39.
Greenland S, Schlesselman JJ, Criqui MH: The fallacy of employing standardized regression coefficients and correlations as measures of effect.
Am J Epidemiol 1986, 2:203-208.
Kupek E: Log-linear transformation of binary variables: A suitable input for structural equation modeling.
Structural Equation Modeling: A Multidisciplinary Journal 2005, 12:35-47.
Browne MW: Asymptotically distribution-free methods for the analysis of covariance structures.
British Journal of Mathematical and Statistical Psychology 1984, 37:62-83.
Bollen KA, Stine RA: Bootstrapping goodness-of-fit measures in structural equation models. In Testing structural equation models. Edited by Bollen KA, Long J. Newbury Park, CA: Sage Publications; 1993::111-135.
Joreskog KG, Soerbom D: PRELIS: A preprocessor for LISREL. Chicago, IL: Scientific Software International; 1994.
Muthén BO: Goodness of fit with categorical and other nonnormal variables. In Testing structural equation models. Edited by Bollen KA, Long J. Newbury Park, CA: Sage Publications; 1993::205-234.
Agresti A: An introduction to categorical data analysis. New York: John Wiley & Sons; 1996.
Muthén BO: Beyond SEM: General latent variable modeling.
Behaviormetrika 2002, 29:81-117.
Skrondal A, Rabe-Hesketh S: Generalized latent variable modelling: multilevel, longitudinal and structural equation models. Boca Raton: Chapman & Hall/CRC; 2004.
Muthén B, Muthén B: Mplus, 2.0 [computer program]. Los Angeles, CA: Muthén & Muthén; 2001.
Vause S, Maresh M: Indicators of quality of antenatal care: a pilot study.
Br J Obstet Gynaecol 1999, 106:197-205. PubMed Abstract
SAS Institute: The CALIS procedure. In SAS/STATR User's Guide. Edited by SAS Institute. Cary, NC: SAS Institute; 1989::245-366.
Satorra A, Saris WE: The power of the likelihood ratio test in covariance structure analysis.
Psychometrika 1985, 50:83-90. Publisher Full Text
Dunn G, Everitt B, Pickles A: Modelling covariances and latent variables using EQS. London: Chapman & Hall; 1993.
Bentler PM, Newcomb MD: Linear structural modeling with nonnormal continuous variables. Application: Relations among social support, drug use, and health in young adults. In Statistical models for longitudinal studies of health. Edited by Dwayer JH, Feinlieb M, Lippert P. Oxford: Oxford University Press; 1992::132-160.
Yuan KH, Bentler PM: Robust mean and covariance structure analysis.
British Journal of Mathematical and Statistical Psychology 1998, 51:63-88.
Yuan KH, Bentler PM: Normal theory based test statistics in structural equation modeling.
British Journal of Mathematical and Statistical Psychology 1998, 51:289-309.
Raftery AE: Bayesian model selection in structural equation models. In Testing structural equation models. Edited by Bollen KA, Long J. Newbury Park, CA: Sage Publications; 1993::163-180.
Agresti A, Caffo B: Measures of relative model fit.
Computational Statistic and Data Analysis 2002, 39:127-136. Publisher Full Text

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2288/6/13/prepub

Interpreting Regression Output

2013-11-12T23:20:00.002+05:30

Interpreting Regression Output

Introduction

This guide assumes that you have at least a little familiarity with the concepts of linear multiple regression, and are capable of performing a regression in some software package such as Stata, SPSS or Excel. You may wish to read our companion page Introduction to Regression first. For assistance in performing regression in particular software packages, there are some resources at UCLA Statistical Computing Portal .

Brief review of regression

Remember that regression analysis is used to produce an equation that will predict a dependent variable using one or more independent variables. This equation has the form

Y = b1X1 + b2X2 + ... + A

where Y is the dependent variable you are trying to predict, X1, X2 and so on are the independent variables you are using to predict it, b1, b2 and so on are the coefficients or multipliers that describe the size of the effect the independent variables are having on your dependent variable Y, and A is the value Y is predicted to have when all the independent variables are equal to zero.

In the Stata regression shown below, the prediction equation is price = -294.1955 (mpg) + 1767.292 (foreign) + 11905.42 - telling you that price is predicted to increase 1767.292 when the foreign variable goes up by one, decrease by 294.1955 when mpg goes up by one, and is predicted to be 11905.42 when both mpg and foreign are zero.

Coming up with a prediction equation like this is only a useful exercise if the independent variables in your dataset have some correlation with your dependent variable. So in addition to the prediction components of your equation--the coefficients on your independent variables (betas) and the constant (alpha)--you need some measure to tell you how strongly each independent variable is associated with your dependent variable.

When running your regression, you are trying to discover whether the coefficients on your independent variables are really different from 0 (so the independent variables are having a genuine effect on your dependent variable) or if alternatively any apparent differences from 0 are just due to random chance. The null (default) hypothesis is always that each independent variable is having absolutely no effect (has a coefficient of 0) and you are looking for a reason to reject this theory.

P, t and standard error

The t statistic is the coefficient divided by its standard error. The standard error is an estimate of thestandard deviation of the coefficient, the amount it varies across cases. It can be thought of as a measure of the precision with which the regression coefficient is measured. If a coefficient is large compared to its standard error, then it is probably different from 0.

How large is large? Your regression software compares the t statistic on your variable with values in the Student's t distribution to determine the P value, which is the number that you really need to be looking at.The Student's t distribution describes how the mean of a sample with a certain number of observations (your n) is expected to behave. For more information on the t distribution, look at this web page .

If 95% of the t distribution is closer to the mean than the t-value on the coefficient you are looking at, then you have a P value of 5%. This is also reffered to a significance level of 5%. The P value is the probability of seeing a result as extreme as the one you are getting (a t value as large as yours) in a collection of random data in which the variable had no effect. A P of 5% or less is the generally accepted point at which to reject the null hypothesis. With a P value of 5% (or .05) there is only a 5% chance that results you are seeing would have come up in a random distribution, so you can say with a 95% probability of being correct that the variable is having some effect, assuming your model is specified correctly.

The 95% confidence interval for your coefficients shown by many regression packages gives you the same information. You can be 95% confident that the real, underlying value of the coefficient that you are estimating falls somewhere in that 95% confidence interval, so if the interval does not contain 0, your P value will be .05 or less.

Note that the size of the P value for a coefficient says nothing about the size of the effect that variable is having on your dependent variable - it is possible to have a highly significant result (very small P-value) for a miniscule effect.

Coefficients

In simple or multiple linear regression, the size of the coefficient for each independent variable gives you the size of the effect that variable is having on your dependent variable, and the sign on the coefficient (positive or negative) gives you the direction of the effect. In regression with a single independent variable, the coefficient tells you how much the dependent variable is expected to increase (if the coefficient is positive) or decrease (if the coefficient is negative) when that independent variable increases by one. In regression with multiple independent variables, the coefficient tells you how much the dependent variable is expected to increase when that independent variable increases by one, holding all the other independent variables constant. Remember to keep in mind the units which your variables are measured in.

Note: in forms of regression other than linear regression, such as logistic or probit, the coefficients do not have this straightforward interpretation. Explaining how to deal with these is beyond the scope of an introductory guide.

R-Squared and overall significance of the regression

The R-squared of the regression is the fraction of the variation in your dependent variable that is accounted for (or predicted by) your independent variables. (In regression with a single independent variable, it is the same as the square of the correlation between your dependent and independent variable.) The R-squared is generally of secondary importance, unless your main concern is using the regression equation to make accurate predictions. The P value tells you how confident you can be that each individual variable has some correlation with the dependent variable, which is the important thing.

Another number to be aware of is the P value for the regression as a whole. Because your independent variables may be correlated, a condition known as multicollinearity, the coefficients on individual variables may be insignificant when the regression as a whole is significant. Intuitively, this is because highly correlated independent variables are explaining the same part of the variation in the dependent variable, so their explanatory power and the significance of their coefficients is "divided up" between them.

Correlation and Causation

2013-11-12T23:19:00.002+05:30

Correlation and Causation

What are correlation and causation and how are they different?

Two or more variables considered to be related, in a statistical context, if their values change so that as the value of one variable increases or decreases so does the value of the other variable (although it may be in the opposite direction).

For example, for the two variables "hours worked" and "income earned" there is a relationship between the two if the increase in hours worked is associated with an increase in income earned. If we consider the two variables "price" and "purchasing power", as the price of goods increases a person's ability to buy these goods decreases (assuming a constant income).

Correlation is a statistical measure (expressed as a number) that describes the size and direction of a relationship between two or more variables. A correlation between variables, however, does not automatically mean that the change in one variable is the cause of the change in the values of the other variable.

Causation indicates that one event is the result of the occurrence of the other event; i.e. there is a causal relationship between the two events. This is also referred to as cause and effect.

Theoretically, the difference between the two types of relationships are easy to identify — an action or occurrence can cause another (e.g. smoking causes an increase in the risk of developing lung cancer), or it can correlate with another (e.g. smoking is correlated with alcoholism, but it does not cause alcoholism). In practice, however, it remains difficult to clearly establish cause and effect, compared with establishing correlation.

Why are correlation and causation important?

The objective of much research or scientific analysis is to identify the extent to which one variable relates to another variable. For example:

Is there a relationship between a person's education level and their health?
Is pet ownership associated with living longer?
Did a company's marketing campaign increase their product sales?

These and other questions are exploring whether a correlation exists between the two variables, and if there is a correlation then this may guide further research into investigating whether one action causes the other. By understanding correlation and causality, it allows for policies and programs that aim to bring about a desired outcome to be better targeted.

How is correlation measured?

For two variables, a statistical correlation is measured by the use of a Correlation Coefficient, represented by the symbol (r), which is a single number that describes the degree of relationship between two variables.

The coefficient's numerical value ranges from +1.0 to –1.0, which provides an indication of the strength and direction of the relationship.
If the correlation coefficient has a negative value (below 0) it indicates a negative relationship between the variables. This means that the variables move in opposite directions (ie when one increases the other decreases, or when one decreases the other increases).

If the correlation coefficient has a positive value (above 0) it indicates a positive relationship between the variables meaning that both variables move in tandem, i.e. as one variable decreases the other also decreases, or when one variable increases the other also increases.

Where the correlation coefficient is 0 this indicates there is no relationship between the variables (one variable can remain constant while the other increases or decreases).

While the correlation coefficient is a useful measure, it has its limitations:Correlation coefficients are usually associated with measuring a linear relationship.
For example, if you compare hours worked and income earned for a tradesperson who charges an hourly rate for their work, there is a linear (or straight line) relationship since with each additional hour worked the income will increase by a consistent amount.

If, however, the tradesperson charges based on an initial call out fee and an hourly fee which progressively decreases the longer the job goes for, the relationship between hours worked and income would be non-linear, where the correlation coefficient may be closer to 0.

Care is needed when interpreting the value of 'r'. It is possible to find correlations between many variables, however the relationships can be due to other factors and have nothing to do with the two variables being considered.
For example, sales of ice creams and the sales of sunscreen can increase and decrease across a year in a systematic manner, but it would be a relationship that would be due to the effects of the season (ie hotter weather sees an increase in people wearing sunscreen as well as eating ice cream) rather than due to any direct relationship between sales of sunscreen and ice cream.

The correlation coefficient should not be used to say anything about cause and effect relationship. By examining the value of 'r', we may conclude that two variables are related, but that 'r' value does not tell us if one variable was the cause of the change in the other.
How can causation be established?

Causality is the area of statistics that is commonly misunderstood and misused by people in the mistaken belief that because the data shows a correlation that there is necessarily an underlying causal relationship .

The use of a controlled study is the most effective way of establishing causality between variables. In a controlled study, the sample or population is split in two, with both groups being comparable in almost every way. The two groups then receive different treatments, and the outcomes of each group are assessed.

For example, in medical research, one group may receive a placebo while the other group is given a new type of medication. If the two groups have noticeably different outcomes, the different experiences may have caused the different outcomes.

Due to ethical reasons, there are limits to the use of controlled studies; it would not be appropriate to use two comparable groups and have one of them undergo a harmful activity while the other does not. To overcome this situation, observational studies are often used to investigate correlation and causation for the population of interest. The studies can look at the groups' behaviours and outcomes and observe any changes over time.

The objective of these studies is to provide statistical information to add to the other sources of information that would be required for the process of establishing whether or not causality exists between two variables.

Developing a Hedonic Regression Model For Camcorders In the U.S. CPI

2013-11-12T23:18:00.003+05:30

Developing a Hedonic Regression Model For Camcorders In the U.S. CPI

Background

The Bureau of Labor Statistics (BLS) has been conducting research into extending the use of hedonic regression models for quality adjustment purposes to additional items within the Consumer Price Index (CPI). Hedonic models estimate values for individual characteristics bundled together to form a good or service(2). This allows the CPI to calculate the value of quality change between two items. The CPI commodity analyst will use the parameter estimates obtained from the hedonic model to adjust the price change used in index calculations in instances where the new item and old item differ in quality. Consumer electronics manufacturers are constantly improving the quality of consumer electronic items in an effort to remain competitive. Quality change often occurs at the time manufacturers introduce new models to replace previous year's models. Video cameras (or camcorders) are one of the items chosen for hedonics research due to recent advances in camcorder technology.

Camcorders are included in the Other Video Equipment CPI item stratum (RA03) along with video cassette recorders (VCRs), DVD players, satellite video products and other miscellaneous video products. Camcorders have an estimated 31 percent of the weight within Other Video Equipment. During the time period from December 1997 to December 1999 the Other Video Equipment index decreased 26.5 percent. The average monthly decline was 1.1 percent. (The Other Video Equipment index was redefined for the January 1998 CPI revision so a longer term comparison is not available. Prior to the revision, camcorders were included in the Video Products Other than Televisions index. This index decreased 38.7 percent from December 1988 to December 1997 — an average monthly decline of 0.4 percent.)

Data and Regression Model

Part of the CPI hedonic initiative called for collecting additional data by CPI field economists. The CPI camcorder sample size was deemed insufficient for regression modeling purposes. Based on current CPI sampling procedures, CPI statisticians designed a supplemental sample for hedonic modeling purposes only. The new sample added 190 outlets with 2 observations assigned in each outlet. The final sample included 350 specially collected observations and 130 observations from the CPI sample. The field economists were unable to collect data for 8 percent of the supplemental sample.

Formulating a hedonic regression model for high tech goods such as camcorders is not a simple task. Determining which camcorder characteristics contribute to the price of a camcorder is difficult due to rapid technical improvements. Manufacturers and retailers of high tech goods further complicate the situation by using different names for the same feature, and retailers have limited information available at the camcorder displays. Consumers are left at the mercy of the retailers' sales personnel for technical information. For the purpose of formulating the camcorder hedonic regression model, data was taken from various sources. CPI field economists collected the primary data. The manufacturer model numbers obtained by the data collectors were matched with specifications provided by the manufacturer internet sites. (Those model numbers that did not match any manufacturer model numbers were dropped from the data set.) Further research helped to develop an a priori model — comparing retailer advertisements, determining which attributes are consistently reported by manufacturers, reviewing consumer magazines and websites, and reading the feedback provided by the data collectors.

Consumers researching camcorders for a purchase are typically advised to first select a format. There are currently five analog camcorder formats available: full-size VHS, 8 millimeter (8mm), Hi-8, VHS-C (compact VHS), and super VHS-C (S-VHS-C); and two digital formats available: mini digital (miniDV), and digital 8 millimeter (digital 8mm).

The digital formats are currently the most technologically advanced. Their picture and audio quality is superior to analog camcorders, but digital camcorders convert their signal to analog for playback on non-digital television sets. Digital computer editing and digital still capability are other technological advances popular with consumers. The miniDV is the smaller of the two formats. Its tapes (although relatively expensive) are the smallest and have immense storage capacity. The digital 8mm is Sony's latest format and enabled Sony to reach consumers shopping the low-end digital market(3). The digital 8mm tapes are compatible with the analog Hi-8 and 8mm formats.

Full-size VHS camcorders have been available the longest. This is the bulkiest camcorder format and usually has to be held on your shoulder to maintain a steady picture. The appeal to this camcorder is that it uses the standard VHS VCR tape for recording; it does not require any adapters for playback in your VCR.

VHS-C produces video comparable to the full-size VHS but is a much smaller camcorder. In order to make the camcorder smaller the manufacturers had to reduce the size of the tape. The tape is about the size of two packs of playing cards(4) and must be inserted into an adapter for playback in your VCR. These smaller tapes record less than half as long as the full-size VHS camcorder. Super VHS-C is similar to the VHS-C but produces a higher quality video.

The 8mm camcorder is about the same size as the VHS-C camcorders but can not be played back on a VHS VCR. Instead it is connected to the television set for playback. Its tape has a longer recording time than the VHS-C camcorder and is generally of higher quality than the VHS-C camcorders. The Hi-8 camcorder has a better picture and audio quality than the standard 8mm.

The data set used for estimating the regression model included 453 observations. Several observations in the data set were deleted due to incomplete data and lack of information from secondary sources. Secondary sources were used to verify the majority of the data. Other observations were deleted because the outlet eligibility for collection in the CPI was questionable. Several camcorders were collected at rental outlets. At these outlets a consumer pays a fee each month to rent an item. If they maintain these payments over a number of months or years the item is theirs to keep, although greatly overpriced. Since we could not be sure of the accuracy of the reported prices due to the nature of this type of outlet, these observations were deleted.

The natural log of the collected price was specified as the dependent variable. The "transaction" price was preferable to the regular price since it represented the price that a consumer would more likely pay. Sale prices were reported for 32 percent of the data set. The mean collected price for observations with a sale price was not much lower than the mean collected price for observations with a regular price (mean for sale priced observations was $674.33 compared to $676.48 for regular priced observations). However, comparing the mean sale price versus the mean regular price across the various formats found that the sale prices for almost all of the formats were lower than the regular prices (full-size VHS was the only exception). Since sale price is a price factor a dummy variable for sale price was included as an explanatory variable in the model.

Specifying the model with the variables mentioned thus far yielded the following results:

Variable Name	Parameter Estimate	Standard Error	T for H0: Parameter=0	Tolerance
Intercept	6.158586	0.01823837	337.672
Sale price	-0.063380	0.02151685	-2.946	0.96758929
Camcorder Format:
VHS	-0.076831	0.05650739	-1.360	0.93656894
8mm	-0.022241	0.02731889	-0.814	0.76197381
VHS-C
S-VHS-C	0.242112	0.07634319	3.171	0.94697938
HI-8	0.395480	0.03118752	12.681	0.80072111
Digital 8mm	0.692510	0.02949420	23.480	0.78800195
MiniDV	1.034472	0.03482691	29.703	0.83312266
R2 = 0.7645; Adjusted R2 = 0.7608; F statistic = 206.86; Number of observations = 453

Further examination of the data found that Sony and Canon camcorders had the highest mean prices. Sony comprised 45 percent of the data set and Canon less than 1 percent. Since Sony has a strong reputation among consumers as a high quality brand and accounts for the majority of brands within the data set, Sony was included in the model. Unfortunately this caused problems with the specification of the model. Sony is currently the only maker of digital 8mm camcorders; therefore, the Sony variable is highly correlated with the digital 8mm variable. Thus multicollinearity will be present in the model if both variables are included. Since both these variables are strong price determinants, not including either variable would result in a specification bias; however, including both variables would result in biased parameter estimates. To remedy the situation, combination variables were created for Sony and camcorder type:

Sony digital 8mm
Non Sony miniDV
Sony miniDV
Non Sony 8mm
Sony 8mm
Non Sony Hi8
Sony Hi8

Sony does not make camcorders in the other three formats (full size VHS, S-VHS-C and VHS-C) so it was not necessary to combine these variables with Sony. The remaining brands represented in the data set (Canon, Hitachi, JVC, Panasonic, Proscan, Quasar, RCA, Samsung, and Sharp) were not included in the model. These brands are highly correlated with other variables in the data set and lack of a priori information regarding their usefulness in explaining price prevented their inclusion.

Based on a priori expectations several other variables were included in the model: monitor size in inches, color viewfinder, image stability, and weight without battery is less than or equal to one pound. Most camcorders currently come with both a viewfinder and a monitor. In this data set 81 percent had a monitor. The monitor is much larger than the viewfinder and also functions as a way to playback the video. Some monitors are built into the side of the camcorder and some swivel out for viewing. Monitor size is the only continuous variable included in the model. The monitor size is measured along the diagonal and the values in this data set range from zero (or no monitor included) to 4 inches. Manufacturers and retail outlets heavily advertise image stability as a price factor. This feature compensates for movement in hand held camcorders so the video playback is not shaky. There are two different types of image stabilization, optical and electronic, but reliable data could not be found for the models in the data set to test if image stabilization type is a price factor. Weight without battery is less than one pound was thought to be price determining since manufacturers heavily advertise this fact. Several manufacturers advertised that their lightweight camcorder was the lightest model currently available to the consumer. The final non-control variable in the model is for Joint Photographers Experts Group or JPEG capability. Images are saved in the universal JPEG format, which allows you to view your images on virtually any computer running Windowsâ or Mac OSâ(5). This variable was included to account for an outlier observation. In this data set there was only one observation that had this feature.

Several variables were included in the model that control for type of business and area of the country where the data are collected. Most of these variables behaved as expected. Discount department stores and warehouse outlets usually have the lowest prices and thus a negative parameter estimate. The positive parameter estimate for furniture/appliance outlets is not unexpected since these are typically the smaller, more specialized and usually local outlets that are known for their customer service. Catalog outlets selling electronic items usually sell mostly high-end electronic items and also charge a steep shipping fee. This is accounted for in the large parameter estimate for catalog outlets.

Variable Name	Parameter Estimate	Standard Error	T for H0: Parameter=0	Tolerance
Intercept	5.850072	0.02210676	264.628	.
Sale price	-0.077808	0.01273242	-6.111	0.88223450
Camcorder Format:
Non-Sony 8MM	-0.166143	0.02551081	-6.513	0.51232970
VHS	-0.087593	0.03552068	-2.466	0.75304756
VHS-C
Non-Sony Hi 8MM	0.084923	0.04161441	2.041	0.90206421
Sony 8MM	0.088555	0.02412902	3.670	0.52296749
S-VHS-C	0.090691	0.04396110	2.063	0.90732597
Sony Hi 8 MM	0.379994	0.01894278	20.060	0.78249993
Sony Digital 8	0.586502	0.01753512	33.447	0.70856195
Non-Sony Digital	0.672284	0.02981434	22.549	0.65632066
Sony Digital	0.880286	0.03190287	27.593	0.62014672
Monitor and Viewfinder:
Monitor Size	0.121159	0.00539411	22.461	0.68787524
Black and White Viewfinder
Color Viewfinder	0.093201	0.01477815	6.307	0.60799750
Other Features:
Image Stability	0.071747	0.01882168	3.812	0.45864474
Weight w/o battery < = 1 lb	0.290191	0.04067331	7.135	0.71303600
JPEG file format capable	0.504677	0.12048018	4.189	0.95143846
Control Variables:
Discount Department Store	-0.081571	0.01776653	-4.591	0.80489622
Midwest Region	-0.041152	0.01368860	-3.006	0.89410748
Warehouse Store	-0.038919	0.02664068	-1.461	0.63459473
Appliance Store	0.050699	0.02090357	2.425	0.86476762
Furniture/Appliance Store	0.281176	0.03281548	8.568	0.88232259
Mail Order Catalog	0.369475	0.12367893	2.987	0.90286013
R2 = 0.9276; Adjusted R2 = 0.9241; F Statistic = 263.044; Number of Observations = 453

The parameter estimates in the model above conform with a priori expectations and the R-squared value indicates that almost 93 percent of the variation in the dependent variable is explained by the independent variables. This is quite high for a hedonic regression model calculated using CPI data.

The final model was examined for multicollinearity using the tolerance statistic and pairwise correlations. Two potential problems were found: a correlation of -0.56 between image stability and Sony 8mm; and a correlation of 0.47 between Non-Sony 8mm and warehouse outlet. However, all these variables were still included in the final model. High correlations did preclude the inclusion of another price factor. Stereo audio is the highest quality audio available for camcorders but stereo audio is highly correlated (0.66) with Sony digital 8mm. Including stereo audio in the model caused the magnitude of the parameter estimate for Sony digital 8mm to fall and its standard error to increase.

A much simpler experimental model was specified with number of pixels as the only explanatory variable. Number of pixels is a way to determine the picture resolution. This information was not specifically requested on the checklist but was sometimes reported by the field economists. Secondary source data was used as a supplement. The "pixel" model did yield good results. The parameter estimate for number of pixels was positive and the R-squared value was large. Number of pixels could not be included in the final camcorder regression model since it is highly correlated with several variables. The Circuit City web site warns against "pixel counting". They mention that the number of pixels reported by manufacturers is actually the potential number of pixels and that most camcorders have a similar picture resolution regardless of the number of pixels.

Variable Name	Parameter Estimate	Standard Error	T for H0: Parameter=0	Tolerance
Intercept	5.470421	0.03157592	173.247
Number of pixels	0.002621	0.00008169	32.088	1.0000000
R2 = 0.7001; Adjusted R2 = 0.6994; F statistic = 1029.614; Number of observations = 443

The final camcorder regression model proposed above may appear to be much simpler than one might expect for a complicated high tech good. The main price determinants are the format of the camcorder along with Sony versus other brands. Those features that were included in the final model appear to be the non-gimmick and non-technical features. Several of the features that manufacturers and retailers tout as the "must haves" were not found to be price determinants. For example, based on the results of preliminary regression models, a larger digital zoom ratio did not have any impact on price despite manufacturers' attempts to persuade consumers of its importance. Also, "technical" specifications were not found to be price contributors. For example, the number of pickup devices(6) and lux(7) did not have an impact on price. Most camcorders available to the consumer (opposed to professional camcorders) have the same number of pickup devices. As with number of pixels mentioned above, Circuit City warns consumers against using manufacturers' lux rating as a guide for purchasing a camcorder. They state that "not all manufacturers use the same scale nor do they use the same testing methods" when determining the lux rating. Based on that statement it is easy to see why consumers are confused over the more technical features.

Index Results

In order to determine the impact of using the camcorder hedonic model in the CPI, an experimental Other Video Equipment index was calculated for the six month time period between June and November 1999. The parameter estimates obtained from the model were applied to camcorder substitute items (an item chosen by CPI data collectors to replace the previously collected item when it is no longer available) with quality changes. During the time period examined there were 58 camcorder substitutions. In the published index, 38 percent had the price of the substitute item directly compared with the price of previous item. The price change for the remainder of the substitutions was imputed via the class-mean imputation method(8). For the purpose of calculating the experimental index, the camcorder substitutions were reassessed. Sixty-seven percent of the substitutions were determined to have changes in quality that could be adjusted using the hedonic model. Out of the remaining substitutions, 28 percent of the prices were directly compared and 5 percent of the price changes were imputed via the class-mean imputation method. The substitution comparability ratio (the ratio of directly compared and quality adjusted substitute quotes to the total number of substitute quotes) improved from 38 percent to 95 percent. Most of the quality adjustments were to adjust for changes in monitor size or changes in type and brand. The table below summarizes the specification changes that occurred with the substitutions.

Specification Change	Number of Occurrences
Unknown
Model number change only
Same item (was actually not a substitution)
Monitor size*
Camcorder brand and/or type*
Color viewfinder*
Image stability*
Low weight*
* Note: more than one of these specifications could have changed for a substitution.

From May 1999 to November 1999 the index using the direct hedonic quality adjustments decreased 6.9 percent compared to a 6.7 percent decrease for the published index, a difference of 0.2 percent(9). Chart 1 compares the published versus the quality adjusted index and chart 2 compares the one month changes of the published versus quality adjusted index.

The months that had the biggest differences between the two indexes were August and November. In August the quality adjusted one month index decreased 0.4 percent more than the published index and in November the quality adjusted one month index rose 0.5 percent more. Comparing the (unweighted) mean price changes for the camcorder substitutions in the published index versus the quality adjusted index shows that the quality adjustments caused the mean price change to be lower in the quality adjusted index. The mean price change for substitutions where the price of the old item was directly compared to the price of the new item increased in the experimental quality adjusted index. These substitutions were mainly those where only the model number changed. The mean price change for substitutions with quality changes fell 9.2 percent, bringing the overall substitution mean price change down more than in the published index. The table below compares the mean price changes of the camcorder substitutions.

Camcorder Substitutions from June 1999 to November 1999
	Published Index		Hedonics Index
	Number	Mean Price Change	Number	Mean Price Change
All substitutions		-3.60 %		-4.56 %
Directly compared substitutions		-1.07 %		3.35 %
Quality adjusted substitutions				-9.21 %
Class-mean imputed (noncomparable) substitutions		-5.14 %		13.66 %

The impact of applying the quality adjustments to the Other Video Equipment index was minimal since camcorders are only a part of the index and camcorder substitutions an even smaller part. These results are similar to studies done with the Televisions index and Video Products Other than Televisions index. The table below summarizes previous CPI hedonics studies.

Item	Difference in 12-Month Index
Personal-computers(10)	6.5% Lower
Televisions(11)	0.1% Lower
Video Products Other than Televisions (VCRs only)(12)	0.1% Higher
Audio Equipment(13)	1.4% Higher

Potentially, using hedonic models for camcorders, DVD players, and VCRs (the three most heavily weighted items in the Other Video Equipment index) should yield a larger impact. The application of the camcorder hedonic model does noticeably decrease the number of noncomparable (imputed) substitutions. This alone makes hedonics a worthwhile endeavor.

Notes

(1) The author wishes to thank Charles Fortuna, Paul Liegey, Bill Thompson, Lynn Reese, John Greenlees, and Mary Kokoski for helpful suggestions and Frank Joseph for his assistance in preparing the camcorder data.

(2) See Dennis Fixler, Charles Fortuna, John Greenlees, and Walter Lane, "The Use of Hedonic Regressions to Handle Quality Change: The Experience in the U.S. CPI," 1999, presented at the fifth meeting of the International Working Group on Price Indices

(3) Consumers Digest, November/December 1998, page 84

(4) Circuit City website, "Learn About Camcorders", http://www.circuitcity.com, visited September 30, 1999

(5) Sony website, http://www.sel.sony.com/SEL/consumer/ss5/office/camcorder/digitalvideoproducts/dcr-trv900_specs.shtml, visited August 12, 1999

(6) The part of the camcorder that translates the optical image from the lens into an electrical signal that can be recorded or viewed. Circuit City website, "Learn About Camcorders", http://www.circuitcity.com, visited September 30, 1999

(7) A rating for the amount of light needed to produce a recognizable image. Circuit City website, "Learn About Camcorders", http://www.circuitcity.com, visited September 30, 1999

(8) See Marshall B. Reinsdorf, Paul Liegey, and Kenneth J. Stewart, "New Ways of Handling Quality Change in the U.S. Consumer Price Index," BLS working paper no. 276 (Bureau of Labor Statistics. 1996).

(9) In addition to quality adjusting the camcorder substitutions, some of the imputed price changes for the class-mean substitutions (the noncomparable substitutions) in the Other Video Equipment index were recalculated since the inclusion of the quality adjustments changed the information used in calculating the imputations.

(10) Kenneth J. Stewart and Stephen B. Reed, "CPI Research Series using current methods, 1978-98", Monthly Labor Review, June 1999, pp. 29-38.

(11) Brent R. Moulton, Timothy J. LaFleur, and Karin E. Moses, "Research on Improved Quality Adjustment in the CPI: The Case of Televisions," Proceedings of the fourth meeting of the International Working Group on Price Indices (U.S. Department of Labor, sponsored by the Bureau of Labor Statistics, January 1999), pp. 77-99.

(12) Paul Liegey and Nicole Shepler, "Adjusting VCR prices for quality change: a study using hedonic methods," Monthly Labor Review, September 1999, pp. 22-37.

(13) Mary Kokoski, Keith Waehrer, and Patricia Rozaklis, "Using Hedonic Methods for Quality Adjustment in the CPI: The Consumer Audio Products Component," BLS draft paper

Difference-in-Difference Estimation: Garbage Incinerators and Home Prices

2013-11-12T23:18:00.000+05:30

Difference-in-Difference Estimation: Garbage Incinerators and Home Prices

The price of a a home can be affected by current interest rates, unemployment rates, and a host of other macroeconomic factors. Home prices are also subject to microeconomic externalities in the form of neighborhood characteristics such as quality of schools, crime, and even the proximity to garbage incinerators. A garbage incinerator was built in North Andover, Massachusetts and the goal is to figure out what kind of impact the garbage incinerator had on home prices over the course of two years. Kiel and McClain estimated the effect of the garbage incinerator on home prices and found that there was only a negligible effect which was statistically insignificant.

What is the Difference-in-Difference Estimation?

A linear regression that is used in policy analysis when there exist a treatment and a control group and two time periods before and after.
A more accurate way of verifying that the average differences between treatment and control groups across time are really meaningful.
It is a way of eliminating unobserved heterogeneity, in other words it is a way of eliminating fixed factors that might have an impact between treatment and control groups.

Figure 1: Linear Regression Model with Difference-in-Difference Estimator

The regression for the average difference between Massachusetts home prices before and after the incinerator demonstrate that the incinerator didn’t impact the prices of homes in any significant way. The reason that home prices are lower is probably not because the incinerator was build, but because home prices were lower that maybe the incinerator was built.

Figure 2: This figure shows the difference -in- difference estimation for the treatment group post policy. In other words it shows the average treatment effect of home prices near the incinerators post-policy.

The interpretation of these coefficients are a little tricky, but one thing to keep in mind is that the numbers are in natural logarithmic form since it is better to get figures in percentages.

Coefficient Explanation-All coefficients are in natural logs and have been converted with the natural number e

y81-The change in the average price of homes between 1978 and 1981 that are away from the incinerator
nearinc- The effect of being near the incinerator in 1978.
y81nrinc-Difference in price from being near the incinerator in 1981 compared to 1978
_cons-Value of house in 1978 that is far from the incinerator, in natural logarithm, to convert to regular price =exp(11.28)= $79,221

The coefficient that we are interested in is the one y81nearinc coefficient of – 6.26% with a p-value of 45.3 percent under the hypothesis that y81nearinc is statistically insignificant from zero.
One can cannot reject the hypothesis that living near the newly build incinerator did not cause a decrease in home prices.
There appears to be other factors that are much more important in determining home prices than the presence of an incinerator.

The following regression shows how other factors are much more significant in determining the change in home prices than whether or not there is a new garbage incinerator near by.

Figure 3: After controlling for other factors that are important in determining home prices y8nrinc is still statistically insignificant.

Interpretation of Regression Coefficient Changes and Control Variables:

y81 – the time trend in home prices for the control group is much less pronounced 14% increase as opposed to a 19% increase when you don’t control for relevant variables.
y81nrinc – is slightly larger and still insignificant at the 10% level, indicating that the new incinerator probably didn’t have ANY affect on home prices in a 3 mile radius.
nearinc – goes from being highly statistically significant to becoming statistically insignificant and much smaller in this new regression.
Significant Control Variables- (bath) an additional bathroom adds 12% to a homes value, (area) an extra 100 feet in area adds about 1% to a persons home price, and (age) every 10 years of aging reduces a homes value by about 8.2%.

About these ads

6 Comments on “Difference-in-Difference Estimation: Garbage Incinerators and Home Prices ”

This is a good and accessible blog. The explanations are clear and suitable for both undergrads and graduate students. I have a question which i hope you will be able to help me with.
I am graduate student at the university of Reading in England. I am carrying out a policy evaluation and need help with running the regression for the difference in differences estimator/ effect of the treatment. I will try as best as possible to narrate the policy, what i am doing, and what i am struggling to understand.

The policy, lets call it X, which is the treatment was enacted in 2001 and the outcome Y is a binary variable , so i have decided to look at the effect of X by comparing the effect on the ‘treated- those affected by X’ and the ‘control- those not affected by X’ . In addition, i have controlled for differences between the control and treatment over time by controlling for fixed effects such as sex, earnings, marital status, job etc between the control and the treated. my research is two dimensional. First, i am looking at the pre policy period of 2000 and then comparing it to a post policy reform period of 2002, bearing in mind that the ‘treatment’ was enacted in 2001. Secondly, i am considering pooling together the observations in 2000 and 2001, treating this as the pre policy period, identify the treatment and control for this period, and then compare this to the post policy period of 2002 and 2003, treating this as the post policy period.
so my question is how do i carry out this two aspects of the analysis in a regression?
What i have done so far is reg Y on set of variables and a time dummy ( in the case of the first analysis t=0 if 2000 and 1 if 2002), on the another dummy variable =1 for the treatd and 0 for the control, this gives the difference between the control and treated over time) and then an interaction of t=1 and treated=1 which i generated by generating another variable by multiplying the treatment which equal 1 by the time period which equal 1 as well ) for the estimator coefficient, i.e treatment effect) however, when i tried to do this in Stata, the interaction term was dropped due to collinearity. As such, i am unable to get the difference estimator( the coefficient on the interaction ) by carrying out the regression this way. Is there a way i can solve this problem?
And since i am using a binary dependent variable, can i still use OLS or should i use the probit non ols estimator ?
I appreciate any insights you can offer on this!

Reply
- Hi Becky,
  
  Hope all is well with you. If you have multiple-years of data, you might want to use a fixed effects estimation that you can use on STATA, or a random effects model if you are familiar with them. This might eliminate the multicollinearity that is common when you are using alot of binary terms in a regresison.
  
  I think it would be quiet interesting to use a probit estimate since you have a binary dependent variable. I haven’t done such a thing when testing policy questions, but I am think it ould be done.Depending on how your data varies these non-linear aspect of this estimate might eliminate this multicollinearity problem.
  
  Last, but not least take a good look at the interaction term. You can export the data on Excel and then compare the values of the interaction term with both the treatment binary and the time binary variable. Do something like =a2=b2 and Excel will respond true if a2=b2, then sort to see if you have “True” for every row when you compare the treatment and time variables to the interaction term. If all are true there is your source of multicollinearity.
  
  Good luck and please let me know what you find out. You have sparked my curiosity.

5 functions to do Principal Components Analysis in R

2013-11-12T23:17:00.002+05:30

5 functions to do Principal Components Analysis in R

Principal Component Analysis (PCA ) is a multivariate technique that allows us to summarize the systematic patterns of variations in the data. From a data analysis standpoint, PCA is used for studying one table of observations and variables with the main idea of transforming the observed variables into a set of new variables, the principal components, which are uncorrelated and explain the variation in the data. For this reason, PCA allows to reduce a “complex” data set to a lower dimension in order to reveal the structures or the dominant types of variations in both the observations and the variables.

PCA in R
In R, there are several functions from different packages that allow us to perform PCA. In this post I’ll show you 5 different ways to do a PCA using the following functions (with their corresponding packages in parentheses):

prcomp (stats)
princomp (stats)
PCA (FactoMineR)
dudi.pca (ade4)
acp (amap)

Brief note: It is no coincidence that the three external packages (FactoMineR, ade4, and amap) have been developed by French data analysts, which have a long tradition and preference for PCA and other related exploratory techniques.

No matter what function you decide to use, the typical PCA results should consist of a set of eigenvalues, a table with the scores or Principal Components (PCs), and a table of loadings (or correlations between variables and PCs). The eigenvalues provide information of the variability in the data. The scores provide information about the structure of the observations. The loadings (or correlations) allow you to get a sense of the relationships between variables, as well as their associations with the extracted PCs.

The Data
To make things easier, we’ll use the dataset USArrests that already comes in R. It’s a data frame with 50 rows (USA states) and 4 columns containing information about violent crime rates by US State. Since most of the times the variables are measured in different scales, the PCA must be performed with standardized data (mean=0, variance=1). The good news is that all of the functions that perform PCA come with parameters to specify that the analysis must be applied on standardized data.

Option 1: using prcomp
The function prcomp comes with the default stats package, which means that you don’t have to install anything. It is perhaps the quickest way to do a PCA if you don’t want to install other packages.

 # PCA with function prcomp
 pca1 = prcomp(USArrests, scale. = TRUE)
 pca1$sdev       # sqrt of eigenvalues
 pca1$rotation   # loadings
 pca1$x          # PCs (aka scores)

Option 2: using princomp
The function princomp also comes with the default stats package, and it is very similar to her cousin prcomp. What I don’t like of princomp is that sometimes it won’t display all the values for the loadings, but this is a minor detail.

# PCA with function princomp
pca2 = princomp(USArrests, cor=TRUE)
pca2$sdev        # sqrt of eigenvalues
pca2$loadings    # loadings
pca2$scores      # PCs (aka scores)

Option 3: using PCA
A highly recommended option, especially if you want more detailed results and assessing tools, is the PCA function from the package FactoMineR . It is by far the best PCA function in R and it comes with a number of parameters that allow you to tweak the analysis in a very nice way.

# PCA with function PCA
require(FactoMineR)
pca3 = PCA(USArrests, graph=FALSE)
pca3$eig         # matrix with eigenvalues
pca3$var$coord   # correlations between variables and PCs
pca3$ind$coord   # PCs (aka scores)

Option 4: using dudi.pca
Another option is to use the dudi.pca function from the package ade4 which has a huge amount of other methods as well as some interesting graphics.

# PCA with function dudi.pca
require(ade4)
pca4 = dudi.pca(USArrests, nf=5, scannf=FALSE)
pca4$eig        # eigenvalues
pca4$c1         # loadings
pca4$co         # correlations between variables and PCs
pca4$li         # PCs

Option 5: using acp
A fifth possibility is the acp function from the package amap .

# PCA with function acp
require(amap)
pca5 = acp(USArrests)
pca5$sdev       # sqrt of eigenvalues
pca5$loadings   # loadings
pca5$scores     # scores

Of course these are not the only options to do a PCA, but I’ll leave the other approaches for another post.

PCA plots
Everybody uses PCA to visualize the data, and most of the discussed functions come with their own plot functions. But you can also make use of the great graphical displays of ggplot2. Just to show you a couple of plots, let’s take the basic results from prcomp

Plot of observations

# load ggplot2
require(ggplot2)

# create data frame with scores
scores = as.data.frame(pca1$x)

# plot of observations
ggplot(data=scores, aes(x=PC1, y=PC2, label=rownames(scores))) +
geom_hline(yintercept=0, colour="gray65") +
geom_vline(xintercept=0, colour="gray65") +
geom_text(colour="tomato", alpha=0.8, size=4) +
opts(title="PCA plot of USA States - Crime Rates")

Circle of correlations

# function to create a circle
circle <- 2="" aes="" and="" arrows="data.frame(x1=c(0,0,0,0)," b="" between="" center="c(0,0)," circles="" colour="gray65" coordinates="" corcir="circle(c(0,0)," correlations="as.data.frame(cor(USArrests," cos="" create="" data.frame="" data="correlations," do="" frame="" function="" geom_hline="" geom_path="" geom_segment="" geom_text="" geom_vline="" ggplot="" label="rownames(correlations)))" length="npoints)" npoints="100)" open="" pca1="" pcs="" pi="" r="" return="" sin="" style="border: 0px; font-family: inherit; font-size: inherit; font-style: inherit; font-variant: inherit; line-height: inherit; margin: 0px; padding: 0px; vertical-align: baseline;" tt="" variables="" will="" with="" x2="correlations$PC1," x="PC1," xend="x2," xintercept="0," xlim="" xx="center[1]" y1="c(0,0,0,0)," y2="correlations$PC2)" y="PC2," yend="y2)," yintercept="0," yy="center[1]">.

1,1.1) + ylim(-1.1,1.1) + labs(x="pc1 aixs", y="pc2 axis") + opts(title="Circle of correlations")

About these ads

R commands for PCA

2013-11-12T23:16:00.001+05:30

Here are some R commands for PCA

pcdat = princomp(data) - It does actual job and put the results to pcdat. It will use covariance matrix

pcdat = princomp(data,cor=TRUE) - It will use correlation matrix

summary(pcdat) - It will print standard deviation and proportion of variances for each component

screeplot(pcdat) - It will plot screeplt

biplot(pcdat) or biplot.princomp(pcdat,scale=1) - It will give you biplot

loadings(pcdat) - it will give information how much each variable contribute to each component. For principal components you can ignore loading subsection of the output from this command

pcdat$scores - It will plot scores of each observation for each variable

For further details about this and other R commands type

help.start()

MODEL FORMULAE for predictive analytics

2013-11-12T23:13:00.005+05:30

MODEL FORMULAE

This is a short tutorial on writing model formulae for ANOVA and regression analyses. It will be linked to from those tutorials, but you are welcome to read it just for kicks if you'd like.

R functions such as aov( ), lm( ), and glm( ) use a formula interface to specify the variables to be included in the analysis. The formula determines the model that will be built (and tested) by the R procedure. The basic format of such a formula is...

response variable ~ explanatory variables

The tilde should be read "is modeled by" or "is modeled as a function of." The trick is in how the explanatory variables are given.

A basis regression analysis would be formulated this way...

y ~ x

...where "x" is the explanatory variable or IV, and "y" is the response variable or DV. Additional explanatory variables would be added in as follows...

y ~ x + z

...which would make this a multiple regression with two predictors. This raises a critical issue that must be understood to get model formulae correct. Symbols used as mathematical operators in other contexts do not have their usual mathematical meaning inside model formulae. The following table lists the meaning of these symbols when used in a formula.

symbol	example	meaning
+	+ x	include this variable
-	- x	delete this variable
:	x : z	include the interaction between these variables
*	x * z	include these variables and the interactions between them
/	x / z	nesting: include z nested within x
\|	x \| z	conditioning: include x given z
^	(u + v + w)^3	include these variables and all interactions up to three way
poly	poly(x,3)	polynomial regression: orthogonal polynomials
Error	Error(a/b)	specify the error term
I	I(x*z)	as is: include a new variable consisting of these variables multiplied
1	- 1	intercept: delete the intercept (regress through the origin)

You may have noticed already that some formula structures can be specified in more than one way...

y ~ u + v + w + u:v + u:w + v:w + u:v:w
y ~ u * v * w
y ~ (u + v + w)^3

All three of these specify a model in which the variables "u", "v", "w", and all the interactions between them are included. Any of these formats...

y ~ u + v + w + u:v + u:w + v:w
y ~ u * v * w - u:v:w
y ~ (u + v + w)^2

...would delete the three way interaction.

The nature of the variables--binary, categorial (factors), numerical--will determine the nature of the analysis. For example, if "u" and "v" are factors...

y ~ u + v

...dictates an analysis of variance (without the interaction term). If "u" and "v" are numerical, the same formula would dictate a multiple regression. If "u" is numerical and "v" is a factor, then an analysis of covariance is dictated.

That ought to do if for now. Specific examples will appear in the tutorials devoted to specific analyses.

Structural Equation Modeling in R for the personality project

2013-11-12T23:13:00.002+05:30

Structural Equation Modeling in R

Structural equation models combine measurement models (e.g., reliability) with structural models (e.g., regression). The sem package, developed by John Fox, allows for some basic structural equation models. To use it, add the sem package by using the package manager.
Structural Equation Modeling may be thought of as regression corrected for attentuation. The sem package developed by John Fox uses the RAM path notation of Jack McCardle and is fairly straightforward. The examples in the package are quite straightforward. A text book, such as John Loehlin's Latent Variable Models (4th Edition) is helpful in understanding the algorithm.
For much more detail on using R to do structural equation modeling, see the course notes for sem (primarily using R) available at the syllabus for my sem course. Also see John Fox's notes that he has prepared as a brief description of SEM techniques as an appendix to his statistics text.

#Loehlin problem 2.5 
 obs.var2.5 = c('Ach1',  'Ach2',  'Amb1',  'Amb2',  'Amb3')
 R.prob2.5 = matrix(c(
  1.00 ,  .60  , .30,  .20,   .20,                                               
   .60,  1.00,   .20,   .30,   .10,                                                
   .30,   .20,  1.00,   .70,   .60 ,                                               
   .20,   .30,   .70,  1.00,   .50,                                                
   .20,   .10,   .60,  .50,  1.00), ncol=5,byrow=TRUE)    
  
  
  
  #correlated factors structure (ambition <-> Achievement) 
  model2.5=matrix(c(
        'Ambit ->  Amb1',      'a', NA,
        'Ambit -> Amb2' ,      'b', NA,
        'Ambit -> Amb3' ,      'c', NA,
        'Achieve -> Ach1',     'd', NA,
        'Achieve -> Ach2',     'e', NA,
        'Ambit <-> Achieve',   'f', NA,
        'Amb1 <-> Amb1' ,      'u', NA,
        'Amb2 <-> Amb2' ,      'v', NA,
        'Amb3 <-> Amb3' ,      'w', NA,
        'Ach1 <-> Ach1' ,      'x', NA,
        'Ach2 <-> Ach2' ,      'y', NA,
        'Achieve <-> Achieve',  NA, 1,
        'Ambit <-> Ambit',      NA, 1),
        ncol=3, byrow=TRUE)
 sem2.5= sem(model2.5,R.prob2.5,60, obs.var2.5)
  summary(sem2.5,digits=3)

 Model Chisquare =  9.74   Df =  4 Pr(>Chisq) = 0.0450
 Goodness-of-fit index =  0.964
 Adjusted goodness-of-fit index =  0.865
 RMSEA index =  0.120   90 % CI: (0.0164, 0.219)
 BIC =  -15.1 

 Normalized Residuals
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-5.77e-01 -3.78e-02 -2.04e-06  4.85e-03  3.87e-05  1.13e+00 

 Parameter Estimates
  Estimate Std Error z value Pr(>|z|)                   
a    0.920    0.0924   9.966 0.00e+00    Amb1 <--- 0.0955="" 0.0965="" 0.1138="" 0.1509="" 0.1762="" 0.356="" 0.652="" 0.683="" 0.761="" 0.879="" 1.45e-11="" 1.55e-15="" 1.76e-03="" 3.127="" 4.525="" 4.986="" 6.03e-06="" 6.16e-07="" 6.753="" 7.974="" ach1="" ach2="" achieve="" amb2="" amb3="" ambit="" b="" c="" d="" e="" f=""> Ambit
u    0.153    0.0982   1.557 1.20e-01     Amb1 <--> Amb1
v    0.420    0.0898   4.679 2.88e-06     Amb2 <--> Amb2
w    0.575    0.0949   6.061 1.35e-09     Amb3 <--> Amb3
x    0.228    0.2791   0.816 4.15e-01     Ach1 <--> Ach1
y    0.534    0.1837   2.905 3.67e-03     Ach2 <--> Ach2

 Iterations =  26 
>  #causal structure with errors in achievement
>   #ambition -> achievement
>   
>  model2.51=matrix(c(
+       'Ambit ->  Amb1',      'a', NA,
+       'Ambit -> Amb2' ,      'b', NA,
+       'Ambit -> Amb3' ,      'c', NA,
+       'Achieve -> Ach1',     'd', NA,
+       'Achieve -> Ach2',     'e', NA,
+       'Ambit -> Achieve',   'f', NA,
+       'Amb1 <-> Amb1' ,      'u', NA,
+       'Amb2 <-> Amb2' ,      'v', NA,
+       'Amb3 <-> Amb3' ,      'w', NA,
+       'Ach1 <-> Ach1' ,      'x', NA,
+       'Ach2 <-> Ach2' ,      'y', NA,
+       'Achieve <-> Achieve',  NA, 1,
+       'Ambit <-> Ambit',      NA, 1),
+       ncol=3, byrow=TRUE)




>    
>  sem2.51= sem(model2.51,R.prob2.5,100, obs.var2.5)
>  summary(sem2.51,digits=3)

 Model Chisquare =  9.74   Df =  4 Pr(>Chisq) = 0.0450
 Goodness-of-fit index =  0.964
 Adjusted goodness-of-fit index =  0.865
 RMSEA index =  0.120   90 % CI: (0.0164, 0.219)
 BIC =  -15.1 

 Normalized Residuals
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-5.77e-01 -3.78e-02 -6.33e-06  4.85e-03  3.06e-05  1.13e+00 

 Parameter Estimates
  Estimate Std Error z value Pr(>|z|)                   
a    0.920    0.0924   9.966 0.00e+00    Amb1 <--- 0.0955="" 0.0965="" 0.0982="" 0.1343="" 0.1394="" 0.153="" 0.1781="" 0.381="" 0.638="" 0.652="" 0.761="" 0.821="" 1.20e-01="" 1.45e-11="" 1.557="" 1.55e-15="" 2.01e-06="" 2.731="" 3.98e-06="" 4.613="" 4.752="" 6.32e-03="" 6.753="" 7.974="" ach1="" ach2="" achieve="" amb1="" amb2="" amb3="" ambit="" b="" c="" d="" e="" f="" u=""> Amb1
v    0.420    0.0898   4.679 2.88e-06     Amb2 <--> Amb2
w    0.575    0.0949   6.061 1.35e-09     Amb3 <--> Amb3
x    0.228    0.2791   0.816 4.15e-01     Ach1 <--> Ach1
y    0.534    0.1837   2.906 3.67e-03     Ach2 <--> Ach2

 Iterations =  27 
> #causal structure with errors in achievement
>   #ambition -> achievement
>   
>  model2.52=matrix(c(
+       'Ambit ->  Amb1',      'a', NA,
+       'Ambit -> Amb2' ,      'b', NA,
+       'Ambit -> Amb3' ,      'c', NA,
+       'Achieve -> Ach1',     'd', NA,
+       'Achieve -> Ach2',     'e', NA,
+       'Ambit -> Achieve',   'f', NA,
+       'Amb1 <-> Amb1' ,      'u', NA,
+       'Amb2 <-> Amb2' ,      'v', NA,
+       'Amb3 <-> Amb3' ,      'w', NA,
+       'Ach1 <-> Ach1' ,      'x', NA,
+       'Ach2 <-> Ach2' ,      'y', NA,
+       'Achieve <-> Achieve', 'z', NA,
+       'Ambit <-> Ambit',      NA, 1),
+       ncol=3, byrow=TRUE)
>    
>  sem2.52= sem(model2.52,R.prob2.5,100, obs.var2.5)
Warning message: 
Negative parameter variances.
Model is probably underidentified.
 in: sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars,  
>  summary(sem2.52,digits=3)

 Model Chisquare =  9.74   Df =  3 Pr(>Chisq) = 0.0209
 Goodness-of-fit index =  0.964
 Adjusted goodness-of-fit index =  0.82
 RMSEA index =  0.151   90 % CI: (0.0517, 0.261)
 BIC =  -8.9 

 Normalized Residuals
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-5.77e-01 -3.78e-02 -1.62e-06  4.85e-03  3.43e-05  1.13e+00 

 Parameter Estimates
  Estimate Std Error z value Pr(>|z|)                     
a    0.920    0.0924   9.966 0.00e+00      Amb1 <--- 0.0955="" 0.0965="" 0.0982="" 0.153="" 0.346="" 0.652="" 0.703="" 0.761="" 0.905="" 1.20e-01="" 1.45e-11="" 1.557="" 1.55e-15="" 6.753="" 7.974="" ach1="" ach2="" achieve="" amb1="" amb2="" amb3="" ambit="" b="" c="" d="" e="" f="" nan="" u=""> Amb1
v    0.420    0.0898   4.679 2.88e-06       Amb2 <--> Amb2
w    0.575    0.0949   6.061 1.35e-09       Amb3 <--> Amb3
x    0.228    0.2791   0.816 4.15e-01       Ach1 <--> Ach1
y    0.534    0.1837   2.905 3.67e-03       Ach2 <--> Ach2
z    0.824       NaN     NaN      NaN Achieve <--> Achieve

 Iterations =  27 

 Aliased parameters: d e f z 
Warning message: 
NaNs produced in: sqrt(diag(object$cov)) 
>   #causal structure with errors in achievement
>   #ambition -> achievement
>   #fix achievement to Loehlin answer
>   
>  model2.53=matrix(c(
+       'Ambit ->  Amb1',      'a', NA,
+       'Ambit -> Amb2' ,      'b', NA,
+       'Ambit -> Amb3' ,      'c', NA,
+       'Achieve -> Ach1',     'd', NA,
+       'Achieve -> Ach2',     'e', NA,
+       'Ambit -> Achieve',   'f', NA,
+       'Amb1 <-> Amb1' ,      'u', NA,
+       'Amb2 <-> Amb2' ,      'v', NA,
+       'Amb3 <-> Amb3' ,      'w', NA,
+       'Ach1 <-> Ach1' ,      'x', NA,
+       'Ach2 <-> Ach2' ,      'y', NA,
+       'Achieve <-> Achieve', NA, .873,
+       'Ambit <-> Ambit',      NA, 1),
+       ncol=3, byrow=TRUE)
>    
>  sem2.53= sem(model2.53,R.prob2.5,100, obs.var2.5)
>  summary(sem2.53,digits=3)\
Error: syntax error
> summary(sem2.53,digits=3)

 Model Chisquare =  9.74   Df =  4 Pr(>Chisq) = 0.0450
 Goodness-of-fit index =  0.964
 Adjusted goodness-of-fit index =  0.865
 RMSEA index =  0.120   90 % CI: (0.0164, 0.219)
 BIC =  -15.1 

 Normalized Residuals
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-5.77e-01 -3.78e-02 -1.14e-06  4.85e-03  3.19e-05  1.13e+00 

 Parameter Estimates
  Estimate Std Error z value Pr(>|z|)                   
a    0.920    0.0924   9.966 0.00e+00    Amb1 <--- 0.0955="" 0.0965="" 0.0982="" 0.1303="" 0.1437="" 0.153="" 0.1906="" 0.356="" 0.652="" 0.683="" 0.761="" 0.879="" 1.20e-01="" 1.45e-11="" 1.557="" 1.55e-15="" 2.01e-06="" 2.731="" 3.98e-06="" 4.612="" 4.752="" 6.32e-03="" 6.753="" 7.974="" ach1="" ach2="" achieve="" amb1="" amb2="" amb3="" ambit="" b="" c="" d="" e="" f="" u=""> Amb1
v    0.420    0.0898   4.679 2.88e-06     Amb2 <--> Amb2
w    0.575    0.0949   6.061 1.35e-09     Amb3 <--> Amb3
x    0.228    0.2791   0.816 4.15e-01     Ach1 <--> Ach1
y    0.534    0.1837   2.906 3.67e-03     Ach2 <--> Ach2

 Iterations =  27 
> 


ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ

 #Loehlin problem from table 2-12
 #Note that version a is a classic example of congeneric measurement.
 #Alternatively, this could be thought of as underidentified higher order model

 obs.var2.12 = c('a',  'b',  'c',  'd')
  R.prob2.12 = matrix(c(
  1.00 ,  .30,    .20,   .10,                                               
   .30,  1.00,   .20,   .20,                                               
   .20,   .20,  1.00,   .30,                                            
   .10,   .20,   .30,  1.00), 
   ncol=4,byrow=TRUE)

  model2.12a=matrix(c(
        'g ->  a',      'a1', NA,
        'g -> b' ,      'b1', NA,
        'g -> c' ,      'c1', NA,
        'g -> d',     'd1', NA,
        'a <-> a',      'e1', NA,
        'b <-> b',      'e2', NA,
        'c <-> c',      'e3', NA,
        'd <-> d',      'e4', NA,
        'g <-> g',       NA, 1),
        ncol=3, byrow=TRUE)
  sem2.12a= sem(model2.12a,R.prob2.12,120, obs.var2.12)
 summary(sem2.12a,digits=3)

#a 1 degree of freedom model
model2.12b=matrix(c(
        'g ->  a',      'a1', NA,
        'g -> b' ,      'b1', NA,
        'f -> c' ,      'c1', NA,
        'f -> d',     'd1', NA,
        'a <-> a',      'e1', NA,
        'b <-> b',      'e2', NA,
        'c <-> c',      'e3', NA,
        'd <-> d',      'e4', NA,
        'g <-> g',       NA, 1,
        'f <-> f',       NA, 1,
        'g <-> f',        'fg', NA),
        ncol=3, byrow=TRUE)
  sem2.12b= sem(model2.12b,R.prob2.12,120, obs.var2.12)
 summary(sem2.12b,digits=3)
 
#the following higher level model has 0 degrees of freedom
 model2.12c=matrix(c(
        'g ->  a',      'a1', NA,
        'g -> b' ,      'b1', NA,
        'f -> c' ,      'c1', NA,
        'f -> d',       'd1', NA,
        'a <-> a',      'e1', NA,
        'b <-> b',      'e2', NA,
        'c <-> c',      'e3', NA,
        'd <-> d',      'e4', NA,
        'g <-> g',       NA, 1,
        'f <-> f',       NA, 1,
        'h -> g',        'hg', NA,
        'h -> f',        NA,1,
        'h <-> h',        NA,1),
        ncol=3, byrow=TRUE)
  sem2.12c= sem(model2.12c,R.prob2.12,120, obs.var2.12)
 summary(sem2.12c,digits=3)
 
 
 
 ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ
 
 #Loehlin problem 2.9
 
  obs.var2.09 = c('w',  'x',  'y',  'z')
  R.prob2.09 = matrix(c(
  1.00 ,  .40,    .50,   .30,                                               
   .40,  1.00,   .55,   .35,                                               
   .50,   .55,  1.00,   .40,                                            
   .30,   .35,   .40,  1.00), 
   ncol=4,byrow=TRUE)

  model2.09=matrix(c(
        'g ->  w',      'a1', NA,
        'g -> x' ,      'b1', NA,
        'g -> y' ,      'c1', NA,
        'g -> z',     'd1', NA,
        'w <-> w',      'e1', NA,
        'x <-> x',      'e2', NA,
        'y <-> y',      'e3', NA,
        'z <-> z',      'e4', NA,
        'g <-> g',       NA, 1),
        ncol=3, byrow=TRUE)
  sem2.09= sem(model2.09,R.prob2.09,500, obs.var2.09)
 summary(sem2.09,digits=3)
 
 
  obs.var2.09b = c('w',  'x',  'y',  'z')
  
  R.prob2.09b = matrix(c(
  1.00 ,  .40,    .50,   .30,                                               
   .40,  1.00,   .55,   .35,                                               
   .50,   .55,  1.00,   .40,                                            
   .30,   .35,   .40,  1.00), 
   ncol=4,byrow=TRUE)

  model2.09b=matrix(c(
        'g ->  w',      NA,1,
        'g -> x' ,      'b1', NA,
        'g -> y' ,      'c1', NA,
        'g -> z',     'd1', NA,
        'w <-> w',      'e1', NA,
        'x <-> x',      'e2', NA,
        'y <-> y',      'e3', NA,
        'z <-> z',      'e4', NA,
        'g <-> g',       'e',NA),
        ncol=3, byrow=TRUE)
        
  sem2.09b= sem(model2.09b,R.prob2.09b,500, obs.var2.09b)
 summary(sem2.09b,digits=3)

part of a short guide to R
Version of February 15, 2007
William Revelle
Department of Psychology
Northwestern University

Useful R links

Readings and software:
Comprehensive R Archive Network (CRAN)
An introduction to R
R Studio

Structural Equation modelling:
sem
lavaan
psych for sem
EFA and factor extension (fa)

Multilevel modeling:
Multilevel
Linear and Non Linear Mixed Effects nlme
statsBy

Item Response Models:
Latent Trait Model (ltm)
mirt
mokken
irt by factor analysis (irt.fa)

Delete outliers from analysis or the data set

2013-11-12T23:12:00.002+05:30

Delete outliers from analysis or the data set

There are no specific R functions to remove outliers . You will first have to find out what observations are outliers and then remove them , i.e. finding the first and third quartile (the hinges) and the interquartile range to define numerically the inner fences.

One way of getting the inner fences is to use

id1 <- boxplot.stats="" gnpgrow="" td="">
id2 <- boxplot.stats="" coef="2)</td" gnpgrow="">	Uses a step of 2, instead of the default 1.5
id1$stats	Display the 5 values (see below)
id1$stats[1]	The lower adjacent value
id1$stats[5]	The upper adjacent value

The boxplot.stats function; is a ancillary function that produces statistics for drawing boxplots. It returns among other information a vector stats with five elements: the extreme of the lower whisker, the lower ‘hinge’, the median, the upper ‘hinge’ and the extreme of the upper whisker, the extreme of the whiskers are the adjacent values (last non-missing value, i.e. every value beyond is an outlier.

Outliers are then all values outside the interval id1$stats[1] < and id1$stats[5]

The adjacent values are in fact not the inner fences, they are further outside, but before an outlier. You can of course define the inner fences precisely, you only need the lower and upper hinger (first and third quartile).

lh <- gnpgrow="" probs="0.25)</td" quantile="">	Lower hinge (first quartile)
uh <- gnpgrow="" probs="0.75)</td" quantile="">	Upper hinge (third quartile)
step<- 1.5="" td="" uh-lh="">	Define the step as 1.5×IQR

Outliers are then all values outside the interval lh-step < and lh+step, the logical expression: gnpgrow lh+step can be used to select the outliers for further processing.

Analyzing Big Data with Twitter

2013-08-28T19:33:00.003+05:30

http://blogs.ischool.berkeley.edu/i290-abdt-s12/

Marti Hearst, the course instructor at UC Berkeley, introduces the main concepts for the course, and Gilad Mishne (@gilad) of Twitter describes his goals for the course and provides an introduction to Twitter. (slides for lecture1a) (slides for lecture1b).

Twitter Philosophy and Software Architecture

Othman Laraki (@othman), Twitter’s Vice President for Growth, International and Revenue, on Growing a Human-Scale Service, and Raffi Krikorian (@raffi), the Director of Twitter’s Platform Services group, on the Twitter Software Ecosystem. View the slides for Othman‘s and Raffi‘s talks.

Introduction to Hadoop

Bill Graham (@billgraham), who is active in the Hadoop community and a Pig contributor, gave a very clear and detailed intro to Hadoop and outlined how it is used at Twitter. His slides can be found here.

Introduction to Apache Pig

Jon Coveney (@jco) gives an in-depth tutorial on Apache Pig, including how it interacts with Hadoop. The log analysis group at Twitter uses Pig extensively. Jon’s slides can be found here: (pdf)

Coding to the Twitter API

Rion Snow (@rion) gave an introduction to the Twitter API, including the RESTful API and the streaming API for both Java and Python. See all the slides (no video).

Slide on Sampling the Streaming API: Twitter4J

Detecting Twitter Trends

If you’d like to know how Twitter computes its Trending Topics, Kostas Tsioutsiouliklis (@kostas) shared some of the secrets with the class. He also talked about MinHash algorithms. See his lecture notes.

Real-Time Twitter Search

Brian Larson (@larsonite), the tech lead for search and relevance at Twitter, gives a detailed technical talk about how real-time search works at Twitter.

Splunk’s Software Architecture and GUI for Analyzing Twitter Data

Stephen Sorkin of Splunk described alternative software architecture for processing large data. Splunk also has a sophisticated GUI for analyzing Twitter and other data sources in real time; be sure to watch the last 15 minutes of the video to see the demo. Stephen’s slides: pdf

Twitter’s Social Network

Learn about weak ties, triadic closures, and personal pagerank, and how they all relate to the Twitter social graph from Aneesh Sharma (@aneeshs) in this lecture. Slides here.

Big Learning with Graphs

Joey Gonzalez, a recent PhD from CMU and a postdoc at UC Berkeley, is working on GraphLab, the hot technology for processing huge graphs quickly. There is new a version called GraphChi (for chihuahua) that you can run on your personal computer; so you don’t even need access to EC2 to run it going forward. Slides here.

Twitter Recommendations

Alpa Jain (@alpa), who works on monetization algorithms at Twitter, described SVD and other recommendation algorithms used at Twitter. Alpa’s slides are here: pdf

Security at Twitter and Elsewhere

Kurt Thomas is a former Twitter engineer and a current PhD student at UC Berkeley who studies how the criminal underground conspires to make money via unintended uses of computer systems. He talked about fraud detection for Twitter and other online systems. See his lecture notes.

Information Diffusion on Twitter

Stan Nikolov (@snikolov) of the Twitter Search and Relevance team walked through one particular theoretical model of information diffusion which tries to predict under what conditions an idea stops spreading based on a network’s structure. The slides in his Lecture Notes let you see the Pig scripts in detail, and you can see the video simulatinos that Stan created on his blog.

Introduction to Scalding

On Thursday we learned about an alternative language for analyzing big data: Scalding. It’s built on Scala and is used extensively by the Twitter Revenue group. Oscar Boykin (@posco) presented a lecture that he and Argyris Zymnis (@argyris) put together. See the lecture notes for more details.

Spark: Making Big Data Analytics Interactive and Real-Time

Spark is the hot next thing for Hadoop / MapReduce, and Matei Zaharia (@matei_zaharia), a PhD student in UC Berkeley’s AMP Lab, described how it works and what’s coming next. The key idea is to make analysis of big data interactive and able to respond in real time. Next up in the research agenda is streaming data and blending real time and batch processing. Matei also gave a live demo. (slides here)

Course Wrapup

These lecture notes simply summarized the course at a high level. Thanks for the great semester!

Posted in Uncategorized | Comments Off

Project Talks and Demos: Tues Dec 11 4pm-6pm

Analyzing the Twitter Conversation and Interest Graphs

For assignment 3, students analyzed and compared a portion of the Twitter “conversation graph” and the “interest graph”. Conversations were found by looking for Twitter “@mentions” and interest graph by looking at the friend/follow graphs for a user (finding friends of friends, taking a k-core analysis, and closing the triangles). The attached document highlights many of the students’ work.
One of the most impressive graphs was made by Achal Soni. He used Java and the Twitter4J library to obtain 3000 tweets for 4 rappers (Drake, Kendrick Lamar, J Cole, and Big Sean). He extracted @mentions from these tweets, and created a graph recording edges were between the celebrities and who they were conversing with.

To compare the interest network with the conversation network for these celebrities he looked at the friends of each of these celebrities, and classified all of these relationships into 3 categories per celebrity: users who they followed but didn’t mention, users who they mentioned but didn’t follow, and users in the intersection. This was repeated for each of the celebrities.
To see what it would look like to compose the relationships for all four celebrities into a single network, he input the combined graph into Gephi. The results were surprising, and phenomenal!
The layout algorithm beautifully clusters each of the celebrity’s networks together, and keeps the shared components in the middle. Each network is composed of nodes of different colors (3 colors each). Each of these colors represent which category the node is associated with belongs to (whether it is a friend of the respective celebrity, somebody who he has conversed with, or both). For example, Big Sean’s network is primarily blue, as opposed to dark orange or purple, which means he talks to a lot of people whom he does not follow. It is also interesting to see which users are in the intersection of the different hip hop artists, and a quick examination of a few of these indicate that they are actually mainly other big figures in the hip hop scene (artists and dj’s, etc.).

More Election Day Twitter Analysis and Visualization

Finishing up with last week’s in-class assignment to respond to election-day twitter activity, some other great work from the class:
Priya Iyer processed a collection of tweets just before and after the election returns and pulled out the most frequent words. She then ran them through the wonderful Wordle program to create a visualization where size indicates frequency (color was chosen by the program, and coincidentally (?) choosing blue and red for “obama” and “romney”.
Before the election:

Just after the returns:

Below is a heatmap of the locations of people who tweeted with the hashtag #Ivoted on election day, based on their geolocation information, by the BandHype team.

Derrick Cheng and Katrina Rogan made a video of their app running that shows the proportion of tweets that talk solely about Obama and solely about Romney. Here is a screenshot:

Check it out!
I think this was a great exercise, as we all learned from each other a bit more about how to get our apps online and get faster at doing our coding, which should help with the final projects.

Video Lecture: Spark: Making Big Data Analytics Interactive and Real-Time by Matei Zaharia

Spark is the hot next thing for Hadoop / MapReduce, and yesterday Matei Zaharia, a PhD student in UC Berkeley’s AMP Lab, gave us a terrific lecture about how it works and what’s coming next. The key idea is to make analysis of big data interactive and able to respond in real time. Matei also gave a live demo. Watch here:
Next up in the research agenda is streaming data and blending real time and batch processing. This research team is doing truly amazing work work and we’re very grateful to @matei_zaharia for the clear and impressive lecture! (slides here)

SF Map of Obama Victory Tweets

Another great result of the in-class assignment to do something with tweets on election day (good thing I gave students 24 hours to turn it in!) since Arian Shams and Guarav Shetti managed to run their tweet capture code for 30 minutes right after Obama was declared the President Elect a second time. Great job, students!

Election Tweets From Restaurants

Students in the Twitter class were challenged yesterday to write some code in class to do something interesting with election day tweets. One group leveraged the code they are developing for their final project to quickly build an app that what people are saying about the election who are near (or claim to be near) Berkeley restaurants. Now it is changed to people contemplating or celebrating the results.
Give it a try!

Video Lecture: Intro to Scalding by @posco and @argyris

On Thursday we learned about an alternative language for analyzing big data: Scalding. It’s built on Scala and is used extensively by the Twitter Revenue group. Oscar Boykin presented a lecture that he and Argyris Zymnis put together for us:
Because scalding is built on the functional programming language Scala, it has an advantage over Pig in that you can have the equivalent of user-defined functions directly in your code. See the lecture notes for more details. Be sure watch the video to get all the details especially since Oscar managed to make us all laugh throughout his lecture. Thanks guys!

Video Lecture: Information Diffusion on Twitter by @snikolov

Today Stan Nikolov, who just finished his masters at MIT in studying information diffusion networks, walked us through one particular theoretical model of information diffusion which tries to predict under what conditions an idea stops spreading based on a network’s structure (from the popular Easley and Kleinberg Network book). Stan also gathered a huge amount of Twitter data, processed it using Pig scripts, and graphed the results using Gephi. The video lecture below shows you some great visualizations of the spreading behavior of the data!
The slides in his Lecture Notes let you see the Pig scripts in more detail.
You can see the videos that Stan created on his blog.
For those who want the details before watching the video, this is a threshold-based model for people choosing to do A or B based on what their neighbors did, modeled as a coordination game where if neighbors pick the same thing, they get a payoff. Even though spread of topics on Twitter is not quite the same kind of coordination game, Stan tells us the threshold model is very popular independent of the game-theoretic justification.
Easley and Kleinberg ask what is it about the structure of the network that cause something to keep spreading or stop spreading. They prove that clusters defined in terms of cluster density stop spreading and that, in fact, they are the only thing that stops spreading.
For his data, he got some trending hashtags (the type that are memes or word games, like #ThingsYouSayToYourBestFriend, not the type that are events, like #debates), and recorded how they spread before the hashtag actually becomes a trending topic (so that the dominant mode of spreading is from person to person, not from some exogenous event, or from the hashtag being on the trending topics list).
Information diffusion in networks is a really difficult topic to work on empirically, so Stan, thank you so much for this terrific work!

SAS Vs SPSS Vs R

2013-08-24T15:57:00.003+05:30

Any serious researcher will fight you to the death to convince you that their favourite stats program is THE best stats program. But, there are good and bad things about each of them.

SPSS: Great for people who rarely use statistics, who don’t remember code, who are scared by code, have no time to learn code, who just need to use some of the basic processes with no special adjustment, who work with normal data that never changes. It’s quick and easy to learn and use, and you don’t 4 feet of manuals to find what you’re looking for. You can just click through the menu to find all the basics you’ll ever need. Compared to SAS, SPSS is major el cheapo.

SAS: Great for people who love coding. SAS does have a menu version but if you like menus, then you should be using SPSS. SAS is great for data manipulation as in creating brand new variables, flipping cases into variables and vice versa, and re-running the same bit of code repeatedly with only a tiny change every time. The macros you can write are unlimited and awesome. And, if it doesn’t do the variation of a statistic you want, you can actually program that statistic into SAS. If you like collecting books, SAS can quickly contribute to that addiction. Buy them. You will need them.

R: What? A third choice? Oh yes. R is great for people who want to flaunt how anti-establishment they are. It’s open source and requires a lot of commitment to become competent. But it can do any statistic you’ve ever dreamed of and a billion more. And you can brag that you know hard core statistics programming. Do not attempt to learn it if Excel intimidates you. Otherwise, go and download it right now and get ready for a rip-roaring awesome 6 week holiday. It is free so prepare to salivate. Mmmmmm…. rrrrrrr. There are even a couple of self-help books now so you might want to find one of them. So, if you’re fresh out of school and no longer have the student version of SAS, R will be perfect for you.

What’s my preference? Well, I wish I was competent in R

List of Free Statistical Software and tool

2013-08-24T12:49:00.000+05:30

Description:

A list of links to free statistics programs, including bioinformatics, psycho metrics, econometrics, simulations, database, data mining and spreadsheets software; and some mathematical software (if it is useful in statistics).

The list shouldn’t include demo/student versions (unless highly functional) or temporary trial versions. Scripts to be run within a program shouldn’t be included in this list.

Search Icons
1	\|STAT	General stats package.
2	ADaMsoft	Data mining, data management and reporting. Has web based features.
3	ADePT	Developed to automate and standardize the production of analytical reports.
4	ADMB	Non-linear statistical modeling
5	AFGen	Fragment-based descriptors for chemical compounds (chemometrics).
6	AM	General stats package
7	AnSWR	A software system for coordinating and conducting large-scale, team-based analysis projects that integrate qualitative and quantitative techniques.
8	Antaeus	Data visualisation software.
9	Apophenia	General stats package but with more flexibility to be creative in model building.
10	Arc	Applied Regression Including Computing and Graphics
11	ARMiner	Data mining application specialising in finding association rules.
12	Assistat	General stats package
13	AutoClass	Data mining (clustering) program from NASA’s Ames Research Center.
14	BAMarray	Java software for microarray data. Implements Bayesian Analysis of Variance for Microarrays (BAM)
15	Bayesian Filtering Library	Bayesian software for use in advanced machine and robot control.
16	BiNGO	Bioinformatics software. Biological Network Gene Ontology.
17	BioEstat	Ecological statistics program (only available in Portuguese).
18	Biogeme	Stands for: Bierlaire's Optimization Toolbox for GEV (Generalized Extreme Value) Model Estimation.
19	Biomapper	A GIS-toolkit to model ecological niche and habitat suitability.
20	BKD	Bayesian Knowledge Discoverer
21	BlockTreat	A general frequentist Monte Carlo program for block and treatment tests, tests with matching, k-sample tests, and tests for independence of a random variable from another random variable.
22	BrightStat	General stats package
23	BV4.1	A procedure for the decomposition and seasonal adjustment of monthly and quarterly economic series. Used by the German Federal Statistical Office.
24	Calcugator	Calculator, plotting engine, and programming environment.
25	Caleydo	Bioinformatics data visualisation software.
26	Carrot-squared	Text mining program. An open source search results clustering engine.
27	CCOUNT	Designed for market research purposes, including: data cleaning, manipulation, cross tabulation and data analysis. Similar to SPSS.
28	CDC EZ-Text	Developed to assist researchers create, manage, and analyze semi-structured qualitative databases.
29	Chronux	Developed for the analysis of neural data.
30	Citation Statistics Scraper	An Excel macro to scrap citation statistics from certain websites.
31	CLUTO	Data mining software for clustering high-dimensional datasets.
32	CodeTwo Outlook Export	Tool for exporting data from Microsoft Outlook to CSV text files.
33	COIN-OR	COmputational INfrastructure for Operations Research
34	Conc	Text concordance program.
35	Concorder	Text concordance program. Allows you to take given texts and find out the frequency of words.
36	ConvertData	Converts structured data from one form to another.
37	CropStat	General stats package. Designed for analysis of agricultural field trials data. Formerly known as IRRIStat.
38	CSPro	Census and Survey Processing System (CSPro) is a questionnaire-oriented statistical package for Windows.
39	Cytoscape	Bioinformatics software platform for visualizing molecular interaction networks and integrating these interactions with gene expression profiles and other state data.
40	DAP	Alternative to SAS.
41	DASS-GUI	Calculates all significant closed sets of a given dataset containing the host sets. Also allows additional analyses of the identified closed sets.
42	Database Browser Portable	Connect to databases, browse or modify data, run SQL scripts, export and print data. Works directly with Oracle, MS Sql Server, ODBC, MySql, PostgreSQL, SQLite, Ole DB, Interbase and Firebird.
43	Dataplot	General stats package. Designed for science and engineering.
44	DAVID	Bioinformatics software. Database for Annotation, Visualization and Integrated Discovery.
45	dChip	Bioinformatics software. dChip Software: Analysis and visualization of gene expression and SNP microarrays
46	Deducer	Graphical user interface (GUI) for R
47	Demetra	Time series package.
48	Dia	Diagram drawing program. Allows the drawing of entity relationship diagrams, UML diagrams, flowcharts, network diagrams, and many others.
49	Draco	Econometrics package with spreadsheets.
50	EASE	Bioinformatics software. Expression Analysis Systematic Explorer.
51	EasyReg Int.	Econometrics package.
52	EasySample	Simple sampling program
53	EDGAR	Experimental design programs.
54	ELKI	Environment for developing KDD-applications supported by index-structures. Compares data mining algorithms.
55	EMMIX	Fortran software to fit mixture models.
56	Epi Info	Designed for epidemiologists and other public health and medical professionals.
57	Epidata	Epidemiology package.
58	EqPlot	Equation graph plotter.
59	Erlang Calculator	Erlang B and C calculator. Used in call centers to predict service quality.
60	ErmineJ	Bioinformatics software. Performs analyses of gene sets in expression microarray data or other genome-wide data that results in rankings of genes. A typical goal is to determine whether particular biological pathways are "doing something interesting" in
61	ESA	Event Structure Analysis (ESA) is an on-line Java program that helps you analyze sequential events.
62	ESS	Emacs Speaks Statistics. Add-on package for emacs text editors.
63	EST	Chemometrics program. Excel add-in (xla) for the simulation of chemical equilibria in solution.
64	Esta+	General (and simple) stats package
65	Euler	Similar to Matlab.
66	EVE	Embedded Vector Editor. Vector graphics program.
67	Excellent Analytics	Excel plug-in that lets you import web analytics data from Google Analytics into a spreadsheet.
68	EzANOVA	Simple ANOVA program
69	Factor	Simple factor analysis program
70	Firebird	Relational database management system (RDBMS).
71	First Bayes	A teaching package for elementary Bayesian Statistics.
72	Fityk	Nonlinear least squares curve fitting.
73	FlexArray	Bioinformatics program. Statistical analysis and visualization of microarray expression data. Free to academic and government researchers.
74	FLOPC++	Formulation of Linear Optimization Problems in C++
75	Freemat	Similar to Matlab and IDL.
76	*GPower**	Power and sample size calculator.
77	G7	Econometrics package. Allows the building and usage of data banks.
78	GAP	Groups, Algorithms, Programming - A System for Computational Discrete Algebra.
79	GATE	Text mining program. Stands for: General Architecture for Text Engineering.
80	gCLUTO	Data mining program. Graphical clustering toolkit.
81	gdcalc	Gdcalc is a financial, statistics, scientific and programmers calculator.
82	Gemma	Bioinformatics software. Database and software system for the meta-analysis of gene expression data.
83	GEMS	Bioinformatics software (online forms). Gene Expression Module Sampler.
84	GenePattern	Bioinformatics software. Genomic analysis platform that provides access to more than 125 tools for gene expression analysis, proteomics, SNP analysis and common data processing tasks.
85	GeNIe	Graphical Network Interface. The graphical interface to SMILE.
86	GenMAPP	Bioinformatics software. Designed to visualize gene expression and other genomic data on maps representing biological pathways and groupings of genes.
87	GeoDa	Geospatial analysis and computation.
88	GeoGebra	Geometry software
89	Gephi	Interactive visualization and exploration platform for all kinds of networks and complex systems, dynamic and hierarchical graphs.
90	GGobi	Data visualisation.
91	Gist	Bioinformatics software. Contains software tools for support vector machine classification and for kernel principal components analysis.
92	Gitools	Bioinformatics Program. A framework for analysis and visualization of genomic data. Data and results are represented as browsable heatmaps.
93	GLPK	GNU Linear Programming Kit
94	GMM-GMR	GMM-GMR - Gaussian Mixture Model and Regression
95	Gnumeric	Spreadsheet program.
96	Gnuplot	Command driven graphing program.
97	GOBLIN	Graph optimization software.
98	Goldvarb	Designed to facilitate the statistical analysis of linguistic variables and social variables. Formerly known as Varbrul.
99	GoMiner	Bioinformatics software. A tool for biological interpretation of 'omic' data – including data from gene expression microarrays.
100	Google Plotter	A tool for adding statistical information to Google maps.
101	Gostat	Bioinformatics software (online forms). Find statistically overrepresented Gene-Ontology (GO) terms within a group of genes.
102	gretl	Stands for Gnu regression, econometrics and time-series library. Alternative to Eviews.
103	Grocer	An econometric toolbox for Scilab.
104	GSL	GNU Scientific Library. A numerical library for C and C++ programmers, providing a wide range of mathematical routines (over 1000 functions) such as random number generators, special functions and least-squares fitting.
105	GUIDE	Multi-purpose machine learning algorithm for constructing classification and regression trees.
106	HLM	Hierarchical linear and nonlinear modelling.
107	hMETIS	A set of programs for partitioning hypergraphs such as those corresponding to VLSI circuits.
108	Instats	General stats package.
109	Inverse Symbolic Calculator	An online form, where you enter a number (the answer), and the output is a list of possible sources (such as simple equations or well known constants) for that number.
110	IVEware	Imputation and Variance Estimation.
111	JAGS	Stands for just another Gibbs sampler. Bayesian hierarchical models using Markov Chain Monte Carlo (MCMC).
112	Java Treeview	Extensible Viewer for Microarray Data in the PCL or CDT format.
113	jHepWork	Mathematics and statistics package.
114	JMulTi	Time series and econometrics package.
115	KNIME	Data integration, processing, analysis, and exploration platform.
116	Kyplot	Graphs program.
117	LazStats	Menu driven general statistics program.
118	LimeSurvey	Survey application.
119	Lisp-stat	General stats package with emphasis on graphics.
120	LISREL (student version)	Originally limited to structural equation modelling (SEM) but now has other functions.
121	LiveGraph	A framework for real-time data visualisation, analysis and logging.
122	LiveInternet Site Statistics	A Google Chrome extension that lets you find at a glance the number of visitors and other data about the current webpage.
123	logQL	Flexible data analysis tool that employs a simple query language to generate the reports. Also known as Log Parser QL.
124	MacAnova	General stats package.
125	MacSHAPA	Macintosh-based software environment that supports observational data analysis, including the analysis of video.
126	MailWatcher Statistics	Program for Domino servers, which collects email data (and stores it into an SQL database), and can graph the data.
127	Matrix2png	Bioinformatics software. A simple but powerful program for making visualizations of microarray data and many other data types.
128	Matrixer	Econometrics package.
129	Matvec	Current capabilities range from matrix and vector manipulation to the analysis of linear and generalized linear mixed models.
130	Maxima	A computer algebra system. Mathematics program.
131	METIS	A set of serial programs for partitioning graphs, partitioning finite element meshes, and producing fill reducing orderings for sparse matrices.
132	MGridGen	Multilevel serial & parallel coarse grid construction library.
133	MicroConcord	Text concordance program.
134	MicrOsiris	Statistical and data management package.
135	Microsoft Excel Viewer	Open, view, and print Excel workbooks, even if you don't have Excel installed.
136	Ministep	Constructs Rasch measures from simple rectangular data sets.
137	MINUIT	Physics analysis tool for function minimization.
138	Mondrian	Data visualisation.
139	MONSTER	Minnesota prOteiN Sequence annotaTion sErveR (bioinformatics).
140	MorePower	Statistical power calculator for hypothesis tests.
141	Mx	Matrix algebra interpreter and numerical optimizer for structural equation modeling and other types of statistical modelling of data.
142	MYSTAT	Free SYSTAT version for students.
143	NLTK	Text mining program. Stands for: Natural Language Toolkit.
144	NORM	Multiple imputation of multivariate continuous data under a normal model.
145	Octave	Similar to Matlab.
146	Onto-Tools	Bioinformatics programs. Onto-Express, Onto-Compare, Onto-Design, Onto-Translate, Onto-Miner, Pathway-Express, Promoter-Express, nsSNPCounter, OE2GO, KUTE-BASE.
147	Open Code	A tool for coding qualitative data generated from text information such as interviews, observations or field notes.
148	OpenBUGS	Bayesian inference using Gibbs sampling.
149	OpenEpi	Epidemiology package.
150	OpenOffice	Open source alternative to Microsoft Office. Includes Calc, an alternative to Excel.
151	OpenSHAPA	Data visualisation program, suitable for video, sound, and electronic data logs.
152	Openstat	General stats package.
153	Orange	Data mining through visual programming or Python scripting. Extensions for bioinformatics and text mining.
154	PAFI	A set of programs (LPMiner, SLPMiner and FSG) that can be used to find frequent patterns in large and diverse databases.
155	PAMCOMP	Person-years And Mortality COMputation Program.
156	PARI/GP	Computer algebra system designed for fast computations in number theory (factorizations, algebraic number theory, elliptic curves), and contains functions for mathematical entities such as matrices, polynomials, power series, algebraic numbers etc.
157	ParMETIS	Extends the functionality provided by METIS and includes routines that are especially suited for parallel AMR computations and large scale numerical simulations.
158	PAST	Statistics program for palaeontology.
159	PAW	Physics Analysis Workstation. Designed for physics data analysis. Developed by CERN (European Organization for Nuclear Research).
160	PCP	Pattern Classification Program
161	PINT	Power analysis IN Two-level designs (for determination of standard errors and optimal sample sizes in multilevel designs with 2 levels).
162	Ploticus	Graphs program.
163	PlotLab VCL	Delphi/C++ Builder VCL components library for very fast Data Visualization
164	PopTools	Population modelling program.
165	POSDEM	A program to choose between sampling plans based on the population frame.
166	PostgreSQL	Relational database management system (RDBMS).
167	PQRS	Stands for probabilities, quantiles and random samples.
168	Program Mark	Ecological statistics program. Advanced mark-recapture modelling.
169	PS	Power and sample size calculator.
170	PSPASES	Parallel SPArse Symmetric dirEct Solver (intended for solving linear systems of equations).
171	PSPP	Alternative to SPSS.
172	QCA	Qualitative Comparative Analysis is a special-purpose program designed to analyze quantified data from multiple cases.
173	QPL	Questionnaire Programming Language.
174	QtOctave	A front-end for Octave.
175	Quail	Quantitative analysis in Lisp.
176	R	Widely used alternative to SPlus.
177	R AnalyticFlow	R front-end which allows data analysis by drawing flowcharts.
178	R Commander	GUI interface for R.
179	RapidMiner	Data mining.
180	Red R	Visual programming interface for R.
181	Regress+	Mathematical modeling for the Macintosh.
182	REG-X	Econometrics package which has functions for cointegration analysis, spectral analysis, dynamic modelling, Kalman Filtering and GARCH-M estimation.
183	RKWard	A transparent frontend to R.
184	RMAexpress	Bioinformatics software. A standalone GUI program for Windows (and Linux) to compute gene expression summary values for Affymetrix Genechip® data using the Robust Multichip Average expression summary and to carry out quality assessment using probe-level
185	ROOT	A system providing a set of OO C++ frameworks with the functions needed to handle and analyze large amounts of data. Developed by CERN (European Organization for Nuclear Research).
186	Rosetta	A rough set toolkit for data analysis / data mining.
187	RStudio	An integrated development environment (IDE) for R.
188	Ruby Stats	Port of the statistics libraries from PHPMath.
189	Rule Discovery System	Data mining program.
190	Rundom Pro	General stats package with emphasis on resampling procedures.
191	Rweb	An online form that processes R commands.
192	Sage	Mathematics software.
193	SalStat	General stats package designed for the analysis of science and psychology data.
194	SAM	Bioinformatics software. Sequence Alignment and Modeling System, a collection of flexible software tools for creating, refining, and using linear hidden Markov models for biological sequence analysis
195	SaTScan	Analyzes spatial, temporal and space-time data using the spatial, temporal, or space-time scan statistics.
196	SciDAVis	Scientific Data Analysis and Visualization.
197	SciGraphica	Data visualisation and analysis
198	Scilab	Numerical computation.
199	Shogun	A large scale machine learning toolbox.
200	SignalLab VCL	Delphi/C++ Builder VCL components library, for fast Digital Signal Processing (DSP) and for data visualization.
201	Simfit	Simulation, curve fitting, statistics, and plotting.
202	Simple Solver	Design and analysis of Boolean equations and state machines.
203	Simplex Method Tool	An online form that solves linear programming problems.
204	SimulME	A Java ME (J2ME) application with queuing calculator, stopwatch, Monte Carlo simulation, random number generator, probability distributions, Kolmogorov-Smirnov and Chi-square tests. Compatible with Sony-Ericsson, Samsung, Motorola and most Java ME (J2ME)
205	SIPINA	Statistical classification program.
206	SL Gallery	Statistical distribution graphs and calculations.
207	Smereka	Extensible personal freeform database and personal information manager.
208	SMILE	Structural Modeling, Inference, and Learning Engine. A fully portable library of C++ classes implementing graphical decision-theoretic methods, such as Bayesian net-works and influence diagrams, directly amenable to inclusion in intelligent systems.
209	SNNS	Stuttgart Neural Network Simulator
210	SOCR	Statistics Online Computational Resource
211	SOFA Statistics	Statistics, analysis and reporting program.
212	Software for Permutation Methods: A Distance Function Approach
213	SolverStat	Excel add-in which complements Excel's built-in Solver package by providing additional statistical analysis of nonlinear regressions.
214	Sonar	Survey administration and data collection system.
215	Sonson	Solution for digitizing vector data from images.
216	SpeQ Mathematics	Calculations and graphs.
217	SpreadCE	Spreadsheet program.
218	SQL Dumper	Allows dumping a SQL Server Database to a text file.
219	SSP (Smith's Statistical Package)	General stats package.
220	StatCalc	A Calculator that computes table values of 34 statistical distributions. It also computes moments, and many other statistics.
221	StatEasy	General stats package.
222	Statext	Statistics expressed in text (even the graphs).
223	Statibot	An interactive www-based expert system for basic statistical analysis.
224	Statist	General stats package.
225	Statistical Lab	R based general stats package.
226	Statistics 101	Simulations / resampling methods.
227	STATPerl	General stats package based on Perl.
228	Stattucino	Web spreadsheets with data analysis functions.
229	SUGGEST	Recommendation engine.
230	Swamy TVC Estimation	Econometrics package which performs the time varying coefficient estimation to decompose a set of parameters into their biased and bias free components.
231	TAMS	Text mining program. Text Analysis Mark-up System. Identifies themes in texts.
232	Tanagra	Data mining.
233	TerraView	GIS application that allows database queries and data analysis.
234	Tetrad	Search algorithms and statistical modelling.
235	TextStat	Text concordance program. A program for the analysis of texts.
236	TimesLab	Time series analysis program.
237	Tinn-R	R code text editor
238	TMeV	Bioinformatics software. Normalized and filtered expression files can be analyzed using TIGR Multiexperiment Viewer (MeV). MeV is a versatile microarray data analysis tool, incorporating sophisticated algorithms for clustering, visualization, classificat
239	TMS FlexAnalytics	An app to get statistics from google analytics and chart it automatically.
240	TSW	Time series and econometrics package.
241	TWOMOK	Non-parametric scale analysis for two-level data (ecometrics).
242	UIMA	Text mining program. Stands for: Unstructured Information Management Architecture.
243	VisiCube	Data visualisation tool.
244	ViSta	Visual statistics system with work maps, guide maps and interactive graphics.
245	VStar	Data visualisation and analysis tool for astronomy.
246	VSUstat	VSUstat - An online statistics program with functions for experimental design and has a layout similar to SPSS. Requires a Google Chrome or Mozilla Firefox browser, and an account log-in.
247	wCLUTO	Web-enabled data clustering application that is designed for the clustering and data-analysis requirements of gene-expression analysis (bioinformatics).
248	Weft QDA	A tool to assist in the analysis of textual data such as interview transcripts, written texts and fieldnotes.
249	Weka	Machine learning software written at the University of Waikato.
250	WinBUGS	Bayesian analysis using Markov chain Monte Carlo methods.
251	WinIDAMS	General stats package.
252	Winpepi	Epidemiology package.
253	WinSAAM	Windows version of CONSAAM. Simulation, analysis and modeling software for biology.
254	WordNet	A large lexical database of English.
255	X-12-ARIMA	Seasonal adjustment software produced, distributed, and maintained by the US Census Bureau.
256	XDAT	Visualization and analysis of multidimensional data.
257	XL Statistics	A set of Microsoft Excel (ver 97+) workbooks for statistical analysis of data.
258	XL Toolbox	Data analysis add-on for Excel.
259	Statistician Lite	Add-in for Microsoft Excel 2007 and Excel 2010. Formerly xlStatistician Lite.
260	Xnumbers	Multi precision floating point computing and numerical methods for Excel.
261	XploRe	Interactive statistical methods and data exploration.
262	Xpro	Exact parametric inference.
263	YACAS	Computer algebra system for symbolic manipulation of mathematical expressions.
264	YASSPP	A web-based server for predicting the secondary structure elements of a protein sequence (bioinformatics).
265	Yorick	Interpreted programming language for scientific simulations or calculations, postprocessing or steering large simulation codes, interactive scientific graphics, and reading, writing, or translating large files of numbers.
266	Zaitun Time Series	Time series software. Has the capability to deal with stock market data.
267	Zelig	Based on R. Automates various processes.
268	Zeptoscope	Java program for drawing plots.
269	Sigma Magic Software	Powerful and easy to use template based statistical analysis package built on top of Microsoft Excel
270	Salstat2	Derived from Salstat. It can be scripted by using python, it can make interactive charts, manipulate Microsoft Excel and also you can import and export data
271	Spatial calculator	An online calculator for spatial statistical functions
272	ILNumerics	Numerical Math Library for .NET / C#, with Statistical Toolboxes and rich Visualization Features.

Reading PDF files in R

2013-08-24T12:15:00.001+05:30

the basic process was to download
xpdf and add it to the windows path, thus making it accessable from
within R. Two methods follow:

Method One (easiest) - using the awesome ?system command:

(1) Download xpdf (whichever is the latest version):
ftp://ftp.foolabs.com/pub/xpdf/xpdf-3.02pl4-win32.zip
(2) Unzip it

# system(paste("[app]", "[pdf file]"), wait = FALSE)

> system(paste('"C:/Program Files/xpdf/pdftotext.exe"', '"C:/Documents and Settings/tony/Desktop/test/r-intro.pdf"'), wait=FALSE)

Method Two - if you want to use the tm package like I did last year,
?readPDF requires the following (not documented anywhere that I know
of, but this is what you do):

(1) Download xpdf (whichever is the latest version):
ftp://ftp.foolabs.com/pub/xpdf/xpdf-3.02pl4-win32.zip
(2) Unzip it
(3) Download the Redmond utility for adding files to your windows path
(free version button is in the top left of the page):
http://redmondlab.googlepages.com/path
(4) Unzip it
(5) Open the 'Redmond Path' application.
(6) Click on the green plus in the top left hand corner '+'.
(7) Naviagate to the folder which contains the files: C:/../
xpdf-3.02pl4-win32
(8) Add it and click Ok.

Then you can can do something like:

> library(tm)
> my.path <- and="" br="" esktop="" here="" in="" ocuments="" pdfs="" put="" settings="" tony="" your="">> Corpus(DirSource(my.path), readerControl = list(reader=readPDF))

There are some limitations to how well the conversions work depending
on the pdf file, but it was so long ago now that I'm afraid I don't
remember the details.

Reading and Text Mining a PDF File in R

2013-08-24T12:13:00.001+05:30


# download pdftotxt from 

# ftp://ftp.foolabs.com/pub/xpdf/xpdfbin-win-3.03.zip

# and extract to your program files folder



# here is a pdf for mining

url <- span=""> "http://www.noisyroom.net/blog/RomneySpeech072912.pdf"

dest <- span=""> tempfile(fileext = ".pdf")

download.file(url, dest, mode = "wb")



# set path to pdftotxt.exe and convert pdf to text

exe <- span=""> "C:\\Program Files\\xpdfbin-win-3.03\\bin32\\pdftotext.exe"

system(paste("\"", exe, "\" \"", dest, "\"", sep = ""), wait = F)



# get txt-file name and open it  

filetxt <- span=""> sub(".pdf", ".txt", dest)

shell.exec(filetxt); shell.exec(filetxt)    # strangely the first try always throws an error..





# do something with it, i.e. a simple word cloud 

library(tm)

library(wordcloud)

library(Rstem)



txt <- span=""> readLines(filetxt) # don't mind warning..



txt <- span=""> tolower(txt)

txt <- span=""> removeWords(txt, c("\\f", stopwords()))



corpus <- span=""> Corpus(VectorSource(txt))

corpus <- span=""> tm_map(corpus, removePunctuation)

tdm <- span=""> TermDocumentMatrix(corpus)

m <- span=""> as.matrix(tdm)

d <- span=""> data.frame(freq = sort(rowSums(m), decreasing = TRUE))



# Stem words

d$stem <- span=""> wordStem(row.names(d), language = "english")



# and put words to column, otherwise they would be lost when aggregating

d$word <- span=""> row.names(d)



# remove web address (very long string):

d <- span=""> d[nchar(row.names(d)) < 20, ]



# aggregate freqeuncy by word stem and

# keep first words..

agg_freq <- span=""> aggregate(freq ~ stem, data = d, sum)

agg_word <- span=""> aggregate(word ~ stem, data = d, function(x) x[1])



d <- span=""> cbind(freq = agg_freq[, 2], agg_word)



# sort by frequency

d <- span=""> d[order(d$freq, decreasing = T), ]



# print wordcloud:

wordcloud(d$word, d$freq)



# remove files

file.remove(dir(tempdir(), full.name=T)) # remove files

Building a Bubble Cloud

2013-08-22T19:23:00.002+05:30

For the 2012 Republican and Democratic national conventions, Mike Bostock , Shan Carter , and Matthew Ericson have created a series of visualizations highlighting the words being used in the speeches of both gatherings. These word-cloud-like word bubble clouds (what I’ll call bubble clouds, unless you can think of a better name) serve as a great interface for looking at the differences in the two conventions and for browsing through quotes from the talks.
Check them out here:

While there is a lot that could be discussed about all the little things that contribute to the quality and polish of these visualizations, in this tutorial we will look at some of the implementation details that make them tick.

Bubble Cloud Demo

I’ve created a basic bubble cloud visualization that tries to replicate some of the functionality of the NYT version. Click on the image below to see the demo .

Try clicking on the bubbles, then try dragging them around. Change the source text with the drop down in the upper-left corner.
The source code is available for download and use in your own projects.
This visualization uses D3’s force-directed layout, so if you aren’t familiar with that, you might check out my post on creating animated bubble charts , or the designing interactive network visualizations tutorial on Flowing Data.
While I won’t be going over the basics of the force layout, hopefully there is still enough in this implementation to keep things interesting. The main topics I want to cover here are:

The use of SVG and plain html components in the same visualization
Saving the state of the visualization using links
Creating a custom gravity effect
Creating a custom collision detection mechanism

Here is all the CoffeeScript that makes this visualization in case you want to follow along in the actual code. Sorry in advanced to the CoffeeScript haters.
Ready to get started? Let’s go!

You got Html in my SVG!

Typically when using D3 , the more advanced visualizations are made with SVG. In a previous tutorial, we dabbled in a SVG-less D3 world , but the trade-offs were steep.
These bubble cloud visualizations actually combine SVG and Html elements elegantly! Specifically, the bubble elements themselves are circle’s inside of a SVG element, but the text on top of them is actually maintained in regular div ’s. As a bonus, both sets of visual elements are backed by the data, so there is very little code duplication or overhead using this structure.
Why would you want this duality? My initial thought when I looked at the implementation was that this was for increased browser compatibility. I figured that in IE8 users might just see the html elements, but miss out on the bubble backgrounds.
This turned out not to be the case. Navigating to the site with an old browser just gets you a static picture. The interactive versions still require browsers with SVG support.
My current guess for why this implementation route was chosen is because of how SVG deals with text wrapping.
Surprisingly, the 1.1 version of the SVG specification does not deal with word wrapping . There are some proposed workarounds , and the 1.2 spec includes the textArea element, which is supposed to solve this deficiency. But none of the solutions look particularly clean, dynamic, or without speed costs.
So the words and phrases the visualizations that will be showing might be on multiple lines. SVG doesn’t support this natively, and it looks like there are some trade-offs with the SVG only workarounds. What should we do?
Why not just implement this part in HTML? This looks to be the path Mike et al. at the New York Times took.
And this combo plater actually works out pretty well. To get started, here’s how we setup the elements that will hold the nodes and labels:

      # node will be used to group the bubbles
      node = svgEnter.append("g").attr("id", "bubble-nodes")
        .attr("transform", "translate(#{margin.left},#{margin.top})")

      # label is the container div for all the labels that sit on top of 
      # the bubbles
      # - remember that we are keeping the labels in plain html and 
      #  the bubbles in svg
      label = d3.select(this).selectAll("#bubble-labels").data([data])
        .enter()
        .append("div")
        .attr("id", "bubble-labels")

Here this refers to the #vis div. For the labels, we add a container div called #bubble-labels. Adding the circles that make up the nodes for the visualization is stuff we have seen before, so lets focus on the labels. First we bind the same data we use to build up the nodes to .bubble-label elements in our label selection:

  # ---
  # updateLabels is more involved as we need to deal with getting the sizing
  # to work well with the font size
  # ---
  updateLabels = () ->
    # as in updateNodes, we use idValue to define what the unique id for each data 
    # point is
    label = label.selectAll(".bubble-label").data(data, (d) -> idValue(d))

So we can still use data bindings, just like in any other D3 visualization!
Now we build up our .bubble-label ’s by entering our (currently empty) selection and appending an anchor and some divs:

    # labels are anchors with div's inside them
    # labelEnter holds our enter selection so it 
    # is easier to append multiple elements to this selection
    labelEnter = label.enter().append("a")
      .attr("class", "bubble-label")
      .attr("href", (d) -> "##{encodeURIComponent(idValue(d))}")
      .call(force.drag)
      .call(connectEvents)

    labelEnter.append("div")
      .attr("class", "bubble-label-name")
      .text((d) -> textValue(d))

    labelEnter.append("div")
      .attr("class", "bubble-label-value")
      .text((d) -> rValue(d))

Here we can see that each label is an a with two div elements inside it. One to hold the word/phrase name, the other to hold the count. We will look at how these anchors work in the next section.
The last thing we need to do when creating these labels is determine the size of the elements holding this text. We want the text to be able to spill out a bit on either size of the bubble, but also to wrap around if it gets too long. Let’s look at one way this can be done. Here’s the code used to style and position the text:

    # first set the font-size and cap the maximum width to enable word wrapping
    label
      .style("font-size", (d) -> Math.max(8, rScale(rValue(d) / 2)) + "px")
      .style("width", (d) -> 2.5 * rScale(rValue(d)) + "px")

    # now we need to adjust if the width set is too big for the text
    # do this using a temporary span
    label.append("span")
      .text((d) -> textValue(d))
      .each((d) -> d.dx = Math.max(2.5 * rScale(rValue(d)), this.getBoundingClientRect().width))
      .remove()

    # reset the width of the label to the actual width
    label
      .style("width", (d) -> d.dx + "px")

    # compute 'dy' - the value to shift the text from the top
    label.each((d) -> d.dy = this.getBoundingClientRect().height)

Let’s analyze this in more detail. As the comments say, we set our font size based on the same scale we are using to scale the bubbles. Good. Bigger bubbles will have bigger text. We also get word wrapping for free by setting the width to an appropriate value.
But if the text is actually too small for the bubble, it will skew to the left. The width will pad it out and make it not line up correctly with the bubble. So we need to find out if the bubble or the word-wrapped text is bigger.
We can do this by adding the text to a temporary span element and grabbing its actual width using getBoundingClientRect . each is used here so we can also set the dx property of our labels, which will be used when positioning them.
With this corrected width known, we can reset our label’s width to account for smaller words.
Finally, we use getBoundingClientRect again to get the amount to shift the label down.
Its a bit confusing, but I think its a great solution to this issue – and now that its been worked out, we can all use it.
To close off this label issue, here is how the dx and dy properties are used in the tick function to position the labels:

    # As the labels are created in raw html and not svg, we need
    # to ensure we specify the 'px' for moving based on pixels
    label
      .style("left", (d) -> ((margin.left + d.x) - d.dx / 2) + "px")
      .style("top", (d) -> ((margin.top + d.y) - d.dy / 2) + "px")

And there you have it. A nice way to use SVG and Html together, keeping with the D3 paradigm of binding data and using the same data backend for both.

A Way to Save State

By saving state, I mean being able to return to or share a particular view of a visualization. The way to do this in web-based interactive visualizations is by modifying the url so as to encode the current state of the visualization in the url itself.
It’s an important detail, and one that is commonly missing from even professionally created visualizations. I liked Bryan Conner’s remark on this subject when discussing the Money Ball WSJ political donor visualization :

Having the ability to share a specific view as a link on Twitter or Facebook is definitely important and, as far as I’m concerned, should be standard for almost any contemporary visualization.

With this bubble cloud implementation, we really only have to keep track of two variables: the text being viewed and the current word selected. I’ll focus on the word selected tracking, which we try to manage in a simple and elegant manner.
Clicking a bubble causes an immediate change to the url. And it is this url change that serves as a signal to update the visualization. This keeps the code clean and neat . Let’s see how it works.
The functionality depends on modifying the url’s hash section which is the text that follows the hash symbol (#). Here is the click callback function, which is executed each time a bubble is clicked:

  # ---
  # changes clicked bubble by modifying url
  # ---
  click = (d) ->
    location.replace("#" + encodeURIComponent(idValue(d)))
    d3.event.preventDefault()

We use location.replace to kill the previous selection (if anything was selected before) and switch to the newly highlighted bubble. The use of the replace function also means the history won’t be polluted with a bunch of versions of the same visualization.
The hash component of the url is accessible via location.hash. What’s more, we can register to be notified when this component of the url changes by using the hashchange event. Combining these insights, the bubble cloud first hooks into hashchange in the visualizations initialization function:

      # automatically call hashchange when the url has changed
      d3.select(window)
        .on("hashchange", hashchange)

And then uses this event to trigger a change in the active bubble.

  # ---
  # called when url after the # changes
  # ---
  hashchange = () ->
    id = decodeURIComponent(location.hash.substring(1)).trim()
    updateActive(id)

  # ---
  # activates new node
  # ---
  updateActive = (id) ->
    node.classed("bubble-selected", (d) -> id == idValue(d))
    # if no node is selected, id will be empty
    if id.length > 0
      d3.select("#status").html("
The word #{id} is now active
")
    else
      d3.select("#status").html("
No word is active
")

This decouples the event handling implementation from the visualization and would allow us to add more callbacks for hashchange to modify other parts of the page if we wanted to. Here, we are just modifying a div, but think of the possibilities!
Of course you can use your own custom event to perform the same decoupling, but I think it is nice having state saving via url modification and user interaction changes wrapped up in this nice little package.

Spreading Gravity Thin

In D3, the gravity component of a force actually draws nodes towards the center of the force layout. It is useful to ensure your nodes don’t fly off the screen. But the default gravity implementation is symmetrical, meaning it pulls on a node’s vertical and horizontal positions equally.
But what if you want a non-symmetrical gravity? Say your visualization is wider then it is tall (like in this example) and you want nodes to spread out along the x axis a bit more? Well, then you can implement your own ‘gravity’ function to push your nodes around how you like!
To get your own forces working on the nodes, we will call them each iteration of the force simulation by binding the force’s tick event to a callback function:

  # The force variable is the force layout controlling the bubbles
  # here we disable gravity and charge as we implement custom versions
  # of gravity and collisions for this visualization
  force = d3.layout.force()
    .gravity(0)
    .charge(0)
    .size([width, height])
    .on("tick", tick)

Every time the tick event is triggered, it will call our function, which happens to be called tick. We also remove the built in gravity and charge forces by setting them to 0. We will implement both of these forces ourselves.
Here is the tick function:

  # ---
  # tick callback function will be executed for every
  # iteration of the force simulation
  # - moves force nodes towards their destinations
  # - deals with collisions of force nodes
  # - updates visual bubbles to reflect new force node locations
  # ---
  tick = (e) ->
    dampenedAlpha = e.alpha * 0.1
    
    # Most of the work is done by the gravity and collide
    # functions.
    node
      .each(gravity(dampenedAlpha))
      .each(collide(jitter))
      .attr("transform", (d) -> "translate(#{d.x},#{d.y})")

    # As the labels are created in raw html and not svg, we need
    # to ensure we specify the 'px' for moving based on pixels
    label
      .style("left", (d) -> ((margin.left + d.x) - d.dx / 2) + "px")
      .style("top", (d) -> ((margin.top + d.y) - d.dy / 2) + "px")

You can see for each node we are calling gravity and collide which perform all the work on our nodes. After which we update our node and label positions to their new x and y coordinates. collide will be in the next section, so let’s focus on gravity now:

  # ---
  # custom gravity to skew the bubble placement
  # ---
  gravity = (alpha) ->
    # start with the center of the display
    cx = width / 2
    cy = height / 2
    # use alpha to affect how much to push
    # towards the horizontal or vertical
    ax = alpha / 8
    ay = alpha

    # return a function that will modify the
    # node's x and y values
    (d) ->
      d.x += (cx - d.x) * ax
      d.y += (cy - d.y) * ay

You can see how the only thing we really tweak is how much alpha affects the x and y components of the gravitational pull.
The alpha parameter comes from D3’s force layout and is the cooling temperature for the layout simulation . alpha starts at 0.1 and decreases as the force simulation continues, getting to around 0.005 before stopping.
Here we are reducing the alpha value to be applied to the x movement by 8. When we multiple the movement towards the center by ax and ay, this makes the y movement stronger. Thus, while the nodes should remain centered vertically, they will be allowed to drift along the x axis.
Redistributing nodes like this allows more content below it to be ‘above the fold’ without running the risk of cutting off the bottom of the bubbles.

Don’t Burst Your Bubbles

As we saw in the bubble chart tutorial , a form of collision avoidance can be implemented in D3 by making the charge associated with each node a function of the size of the node. This provides a visually interesting, organic, experience, where bubbles push on one another in an particle-like way.
The bubbles in these word visualizations act in a subtly but significantly different manner. Here the nodes bounce off one another, maintaining a rigid parameter space around themselves. So how is this affect achieved in D3? By implementing a custom collision detection and avoidance algorithm!
Its actually not as complicated as it might sound. Here is the code:

  # ---
  # custom collision function to prevent
  # nodes from touching.
  # ---
  collide = (jitter) ->
    # return a function that modifies
    # the x and y of a node
    (d) ->
      data.forEach (d2) ->
        # check that we aren't comparing a node
        # with itself
        if d != d2
          # use distance formula to find distance
          # between two nodes
          x = d.x - d2.x
          y = d.y - d2.y
          distance = Math.sqrt(x * x + y * y)
          # find current minimum space between two nodes
          # using the forceR that was set to match the 
          # visible radius of the nodes
          minDistance = d.forceR + d2.forceR + collisionPadding

          # if the current distance is less then the minimum
          # allowed then we need to push both nodes away from one another
          if distance < minDistance
            # scale the distance based on the jitter variable
            distance = (distance - minDistance) / distance * jitter
            # move our two nodes
            moveX = x * distance
            moveY = y * distance
            d.x -= moveX
            d.y -= moveY
            d2.x += moveX
            d2.y += moveY

We compute the distance between each node pair and see if it is less then the minimum distance allowed by our visualization. If so, then we move the nodes away from one another. The fact that D3 stores the current x and y position of each node as part of the data associated with that node makes it easy to get the values required to perform these calculations.
Notice the input variable jitter. What’s it do? We use it to help scale the distance that will be used to move the colliding nodes. Higher values will make this distance larger, smaller values will reduce it.
I’ve connected the value of jitter to a range control underneath the visualization. This lets you explore what the collisions look like when the nodes really start pushing each other hard. The default value, 0.5 visually seems to make the interactions look natural.

The Quad Connection

The implementation as shown above is a brute-force strategy. All nodes are compared with all other nodes, regardless of how close they are in the visualization. While this works decently for a small number of nodes, it will start bogging down as the node count increases.
If we had some rough-grained insight as to which nodes were nearest one another, we might be able to speed things up significantly by only performing these calculations for these nearby nodes.
This is actually how the NYT version of the collision algorithm is implemented. How do they know which nodes are nearby? By using a Quadtree .
A Quadtree is a tree-like data structure in which each internal node has four children (hence the name). You can use a Quadtree to partition up a two-dimensional space such that areas with more nodes will have more partitions. You can then visit the nodes in an efficient manner and stopped traversing once you have made it out of the quadrants that could affect the node in question.
D3 has a sparsely documented Quadtree implementation which is used in the force layout in this manner to make it fast. While this Quadtree implementation is a great option to have, and should be considered for collisions between many nodes, I think the brute-force version provides the same basic idea, without more technical overhead.

Your Own Bubble Cloud in the Sky

That does it for this tutorial. Hopefully this provides a bit more insight into these great pieces from the New York Times (and hopefully they don’t mind me continuing to exploit their great pieces).
Again, the code is on github , so grab it and let your bubble clouds accumulate!
Taking this to the next level would involve splitting the bubbles based on some variable, like the dual convention version .
For that, check out Mike’s great demonstration of how they split the bubbles to get a sense of how to add this kind of visualization.
Enjoy!

Top 10 Sites for Creating Word Clouds

2013-08-22T19:15:00.002+05:30

Top 10 Sites for Creating Word Clouds

"A text cloud or word cloud is a visualization of word frequency in a given text as a weighted list (wikipedia)." This is a concept that is utilized a lot in education for: visualization, brain storming, and mind mapping." While this can be done using paid applications such as Inspiration, there are a number of free sites that accomplish the same.

Top 10 Sites for Creating Word Clouds

Wordle - Probably the oldest and most popular of the word cloud creators. Word clouds can then be edited w/ different fonts, colors, and layouts, and even shared w/ others online.
Tagxedo - A great way to make word clouds from sites, different social networks, and then change shape to meet your needs (ie, star).
Tagul - Similar to Tagxedo, a very easy to use site to create different shaped word clouds.
Worditout - A nice site for making word clouds by just copying/pasting text.
Word Sift - A excellent site for making word clouds and then editing them by altering the criteria in how it's created. Pictures are also generated to best represent your word cloud.
VocabGrabber - A wonderful site similar to Word Sift in making word clouds. However, the clouds generated are color coded for best displaying clouds.
Tag Galaxy - An innovative site that visually represents words in a galaxy style format. A must see!!!
Tag Crowd - A simple to use site for creating word clouds by pasting in text or a web page URL.
Quintura - Is a visual search engine where results are returned in a word cloud style format.
TwitterSheep - A fun site that generates it results by keywords found in bios from a user's followers.

Word cloud in R

2013-08-22T19:14:00.000+05:30

An update to the wordcloud package (2.2) has been released to CRAN. It includes a number of improvements to the basic wordcloud. Notably that you may now pass it text and Corpus objects directly. as in:

#install.packages(c("wordcloud","tm"),repos="http://cran.r-project.org")
library(wordcloud)
library(tm)

wordcloud("May our children and our children's children to a
thousand generations, continue to enjoy the benefits conferred
upon us by a united country, and have cause yet to rejoice under
those glorious institutions bequeathed us by Washington and his
compeers.",colors=brewer.pal(6,"Dark2"),random.order=FALSE)

data(SOTU)
SOTU <- ark2="" colors="brewer.pal(6," function="" random.order="FALSE)</pre" removewords="" stopwords="" tm_map="" tolower="" wordcloud="" x="">



This bigest improvement in this version though is a way to make your 
text plots more readable. A very common type of plot is a scatterplot, 
where instead of plotting points, case labels are plotted. This is 
accomplished with the text function in base R. Here is a simple 
artificial example:


states <- aine="" alifornia="" ansas="" aryland="" ashington="" assachusetts="" awaii="" c="" carolina="" daho="" dakota="" ebraska="" elaware="" ennessee="" ennsylvania="" entucky="" eorgia="" ermont="" est="" evada="" ew="" exas="" hampshire="" hio="" hode="" ichigan="" innesota="" irginia="" isconsin="" island="" ississippi="" issouri="" jersey="" klahoma="" labama="" laska="" llinois="" loc="" lorida="" matrix="" mexico="" ncol="2))" ndiana="" olorado="" onnecticut="" ontana="" orth="" ouisiana="" outh="" owa="" plot="" pre="" regon="" rizona="" rkansas="" rmvnorm="" states="" tah="" text="" type="n" virginia="" yoming="" york="">

Notice how many of the state names are unreadable due to 
overplotting, giving the scatter plot a cloudy appearance. The textplot 
function in wordcloud lets us plot the text without any of the words 
overlapping.


textplot(loc[,1],loc[,2],states)

A big improvement! The only thing still hurting the plot is the fact 
that some of the states are only partially visible in the plot. This can
 be fixed by setting x and y limits, whch will cause the layout 
algorithm to stay in bounds.


mx <- apply="" loc="" max="" min="" mn="" states="" textplot="" xlim="c(mn[1],mx[1]),ylim=c(mn[2],mx[2]))</pre">
Another
 great thing with this release is that the layout algorithm has been 
exposed so you can create your own beautiful custom plots. Just pass 
your desired coordinates (and word sizes) to wordlayout, and it will 
return bounding boxes close to the originals, but with no overlapping.


plot(loc[,1],loc[,2],type="n")
nc <- .5="" cex="50:1/20)</pre" loc="" nc="" states="" text="" wordlayout="">



okay, so this one wasn't very creative, but it begs for some further 
thought. Now we have word clouds where not only the size can mean 
something, but also the x/y position (roughly) and color. Done right, 
this could add whole new layer of statistical richness to the 
visually pleasing but statistically shallow standard wordcloud.







Filed under: R, wordcloud
7 Comments



25Jan/12Off

Beat me to it
So as a follow-up to the last post about the 2010-11 state of the union addresses, I was going to compare the 2011 to the new 2012 speech. It looks like Simply Statistics beat me to it, which is very gratifying. No need to do the work myself  




Filed under: wordcloud
No Comments



10Jan/12Off

Words in Politics: Some extensions of the word cloud
The word cloud is a commonly used plot to visualize a speech or set 
of documents in a succinct way. I really like them. They can 
be extremely visually pleasing, and you can spend a lot of time perusing
 over the words gaining new insights.


That said, they don't convey a great deal of information. From a 
statistical perspective, a word cloud is equivalent to a bar chart of 
univariate frequencies, but makes it more difficult for the viewer to 
estimate the relative frequency of two words. For example, here is a bar
 chart and word cloud of the state of the union address for 2010 and 
2011 combined.




Bar chart of the state of the union addresses for 2010-11



word cloud of the state of the union addresses for 2010-11

Notice that the bar chart contains more information, with the exact 
frequencies being obtainable by looking at the y axis. Also, in the word
 cloud the size of the word both represents the frequency, and the 
number of characters in the word (with longer words being bigger in the 
plot). This could lead to confusion for the viewer. We can therefore see
 that from a statistical perspective that the bar chart is superior.


... Except it isn't ....


The word cloud looks better. There is a reason why every infographic 
on the web uses word clouds. It's because they strike a balance of 
presenting the quantitative information, while keeping the reader 
interested with good design. Below I am going to present some extensions
 of the basic word cloud that help visualize the differences and 
commonalities between documents.



The Comparison Cloud
The previous plots both pooled the two speeches together. Using 
standard word clouds that is as far as we can go. What if we want to 
compare the speeches? Did they talk about different things? If so, are 
certain words associated with those subjects?


This is where the comparison cloud comes in.




Comparison plot

Word size is mapped to the difference between the rates that it 
occurs in each document. So we see that Obama was much more concerned 
with economic issues in 2010, and in 2011 focused more on education and 
the future. This can be generalized fairly naturally. The next figure 
shows a comparison cloud for the republican primary debate in new 
hampshire.



One
 thing that you can notice in this plot is that Paul, Perry and Huntsman
 have larger words than the top tier candidates, meaning that they 
deviate from them mean frequencies more. On the one hand this may be due
 to a single minded focus on a few differentiating issues (..couch.. Ron
 Paul), but it may also reflect that the top tier candidates were asked 
more questions and thus focused on a more diverse set of issues.

The Commonality Cloud
Where the comparison cloud highlights differences, the commonality 
cloud highlights words common to all documents/speakers. Here is one for
 the two state of the union addresses.







Commonality cloud for the 2010-11 SOTU

Here, word size is mapped to its minimum frequency across documents. 
So if a word is missing from any document it has size=0 (i.e. it is not 
shown). We can also do this on the primary debate data...




Republican primary commonality cloud

From this we can infer that what politicians like more than anything else is people  






The wordcloud package
Version 2.0 of wordcloud (just released to CRAN) implements these two types of graphs, and the code below reproduces them.


library(wordcloud)
library(tm)
data(SOTU)
corp <- 2010="" 2011="" as.matrix="" c="" collapse="\n" collected="" colnames="" commonality.cloud="" comparison.cloud="" content="" corp="" corpus="" ectorsource="" for="" function="" in="" library="" list="" max.words="Inf,random.order=FALSE)" n="" names="" nchar="" olorbrewer="" paste="" pre="" r2="" random.order="FALSE)" rcorp="" readlines="" removenumbers="" removepunctuation="" removewords="" repub="" repub_debate.txt="" rterms="" sotu="" sp="=speaker],collapse=" speaker="" splitat="" stopwords="" str_locate_all="" str_sub="" stringr="" stripwhitespace="" strsplit="" term.matrix="" termdocumentmatrix="" tm="" tm_map="" tmp="" tolower="" unlist="" wordcloud="" x="">



Link to republican debate transcript



















Filed under: R, wordcloud

Word cloud sites (Text mining)

2013-08-22T19:12:00.000+05:30

Wordle – (http://www.wordle.net/create)

free
no login required
No Save As option, but take a screen shot before you close to save it. Hit the PrtScn button and paste into another application, such as ppt or Word. If you paste into ppt, you can then save as "Other Formats" and choose jpg from the drop down in the "Save as type..."
keep word phrases together by using the tilde ~ (on the same key as the accent mark, to the left of number 1)
There may be inappropriate content in the Gallery, so have students go directly to http://www.wordle.net/create.

Tagxedo - (http://www.tagxedo.com/) “Tag Clouds with Style”

free
custom shapes and themes
font choices
save as jpg and/or png
can view history to choose the best one

ABC Ya - (http://www.abcya.com/word_clouds.htm)

free
no login required
very similar to Wordle
save as jpg

Tagul - (http://tagul.com/)

requires login with email
custom shapes
multiple fonts usage in one cloud (not available in Wordle)
tag words to act as pointers to url's if embedded on web page

Word It Out (http://worditout.com/)

login required, but don't need this to right-click and save as jpg or PrtScn
cutomizable by size, font, color scheme
allows user to ignore certain words to keep them out of the word cloud

Tag Crowd (http://tagcrowd.com/)

no login
free for nonprofit
shows frequency of words
can upload file to url
save as pdf

Wordsift - (http://www.wordsift.com/)

no login
no save as, so use PrtScn
from Stanford University
allows words to be listed by how common or rare they are
allows words to be listed in ABC order
Wordsift allows the user to click on words to view in an online visual thesaurus with dictionary, google images, and word sentence placement
view words by subject area and cross curricular areas by unique color coding and definitions that relate to specific disciplines
Watch this video lesson on using Wordsift with students on a lesson about Dr. Martin Luther King Junior.

Make Word Mosiac (http://www.imagechef.com/ic/word_mosaic/)

no login or email required
different shapes, colors, and fonts
may have limited use in the classroom
save as jpg (choose from the More button)
for personal use

VocabGrabber – (http://www.visualthesaurus.com/vocabgrabber/)

analyzes text and generates lists of the most useful vocabulary words then displays how those words are used in context
Copy text from a document and paste it into the box, and click Grab Vocabulary! VocabGrabber will automatically create a list of vocabulary from the text, which can be sorted, filtered, and saved. Click on any word in the cloud and a snapshot of the Visual Thesaurus map appears along with definitions for that word, and examples of the word in the text.

TagCloudGenerator – (http://www.tag-cloud.de/)

cannot paste text in this one
results are a moving flash file that can be downloaded

TagCloud - (http://www.tagcloud.com)

currently offline, getting an overhaul
check back to see if it's up