title: Problem with identifying outliers

text: I'm running a biological experiment with rodents, have two groups (each consists of 26 animals), where one is treated with a chemical, and one is control (saline). In one variable, there doesn't seem to be much of a difference between groups. However, when I identify outliers, the data becomes very significant ($p < .01$). I am wondering are my criteria to exclude those data points valid. I've used two methods to identify outliers. First of all, I was identifying outliers separately for control and treatment, as the data seems very different (check histogram attached for the treatment group)—is this wrong or not? To me, it sounds logical, as the groups are different (in the sense they are chemically treated differently). ![Histogram for treated group](http://i.stack.imgur.com/6vTZp.png) First method I mentioned—I checked for data that is 2.5 STDEV from the mean, and excluded those points. The second one is similar: http://www.wikihow.com/Calculate-Outliers—this one gives the same result, and identifies the outliers in the treatment group. However, as I would like to prepare a publication, I need a real "peer- accepted" test for outliers. I have found some tests in StatSoft Statistica, but the box whisker plot does not show outliers there. However, in the Grubbs' test, I get $p < .05$, so it is analytically detected as an outlier, right? So, my questions are: 1) Is it ok to calculate outliers differently for each group, or should I collapse the sample into one group ($N = 52$) and calculate it then, for both groups together? 2) Which method should I use to identify (and exclude) outliers in order for it to be accepted by a critical review? I attach the data points here: https://www.mediafire.com/?0qdsifib0hugd9u To me, it seems that any animal with the score higher than 54 is an outlier. Any help would be greatly appreciated!


title: Seasonal Data with Yearly "Step" Trend

text: Trying to plan for upcoming inventory needs for a clothing company. There is a strong seasonal component to the data. At the beginning of the year, sales are up; at the end, sales are down. Each year, at the beginning of the year, they acquire more sales people, and so they sell more each year. But! this upward "trend" only happens once a year. I tried doing a seasonal time series analysis with a trend, but this gave too high sales at the end of the year. How do I carry out this "step trend" analysis? Thanks!


title: Splitting a variable with nominal and numeric values

text: I have a variable that has both numeric and nominal components. The source has a documentation which helps in identifying which is which and for splitting into their proper components. I will do modeling on these data later on (R or WEKA). There are several of them and one of it is the `size` variable: the numeric components are 0 to 70 and the nominal values are `<10`, `10-20` up to `60-70`. The nominal part corresponds to unclear interviews with the patient, lost records (thus the inaccuracy) and other medical shortcomings. My plan is to create two new variables - `size-num` and `size-nom` that will have the aforementioned to sets respectively as values. Once I modify `size`, I face new problems:               +-------+---------+-------------+-------------+     | id    | size    | size-num    | size-nom    |     +-------+---------+-------------+-------------+     | 1     | >10     | ?           | >10         |     | 2     | 27      | 27          | ?           |     | 3     | 20-30   | ?           | 20-30       |     | 4     | 0       | 0           | 0           |     | 5     | <60     | ?           | <60         |     +-------+---------+-------------+-------------+        * In the numeric variable, how do I represent the removed nominal ones? If I just leave them out to be blanks then I will have to take into account missing values.   * In the nominal variable, how I do represent the numeric ones? I know that I can use them (numeric values) again but this will remove their continuous characteristic.


title: How to extract values of the forecasted times series of auto.arima

text: I'm working actually on the modeling of times series with auto.arima from package forecasts. But i enconter some problems regarding the extraction of values Lo 80, Hi 80, Lo 95, Hi 95 singularly. I have retrieved data to excel with function csv but it 's too long and there are some risks to loose data. I have retrieved only the mean or the point forecast with the function AA11a[["mean"]] as you can see on my codes in R               a<-ts(AA,start=1921,end=2009,frequency=1)     fit11a<-auto.arima(b,allowdrift=TRUE)     AA11a<-forecast(fit11a,50)     plot(AA11a)     plot(AA11a[["mean"]])      I want to ask how can i extract the Lo 80, Hi 80, Lo 95, Hi 95 singularly? I have tried this code AA11a$fcst$AA[,1] but it doesn't work and tell me : Null The other alternative is to send the output AA11a to excel and find singularly variables but before exploring that path i will ask whether somebody has an idea on how to extract different values i have cited above with codes in R. Thanks


title: What is my model statement for mixed models (nlme) in R?

text: I want to know if a covariate for each subject interacts with three types of trials, and the difficulty of those trials. My dependent measures are accuracy and response times (RT). For this question, I’d like to focus on RTs. Traditionally, people in my field have dichotomized the covariate of interest and used ANOVAs for analysis. I would like to treat the covariate as the continuous variable it is, and treat the subjects as random effects. I want to analyze this using mixed-models in `R` (`nlme`). The first 2 trial types can be either easy or hard and the third trial type is a combination of the first 2. These trials can be easy-easy, easy-hard, hard- easy, hard-hard. I expect people who have higher scores on the covariate to show a smaller difference between hard and easy RTs for at least 1 trial type. This is a repeated-measures design with each subject completing 3 blocks of 40 trials of each of the trial types (for trialtypes 1 & 2: 20 easy, 20 hard; for trialtype3, 10 easy-easy, 10 easy-hard, 10 hard-hard, 10 hard-easy). Stated differently, each subject completes 3 blocks of 120 trials with the various trialtypes randomly ordered. Only RTs for correct trials will be analyzed (resulting in an unbalanced design for RT data). Besides a counter-balancing of response keys, this is a completely within-subjects design. To summarize, what is the model (or models) that will allow me to test for interactions between trialtypes, difficulty, and the covariate using `nlme` in `R`?


title: What type of chart is this? (bubble charts with links)

text: Does anybody know what type of chart that is? (some sort of bubble chart but with links between bubbles) http://en.wikipedia.org/wiki/File:PageRanks-Example.svg ![enter image description here](http://i.stack.imgur.com/bhjcQ.png)


title: [Revised]Proving the expected \bold{density} of being the Nth order statistics is decreasing in sample size

text: (Sorry that I've previously formulated the question in a wrong way, which confused everyone including myself. This is a better version of the question. Thanks!) Here's another order statistics question that I wish to ask. **_Question_** Consider $n$ random variables $x_1, x_2,\cdots x_n\overset{iid}{\sim} D$. Where $D$ is some unimodal on 0, symmetric, continuous distribution with a finite variance (P.S. This condition might be overly restrictive, suggestions on loosening it would be greatly appreciated! ). The PDF for nth order statistics is $nF_{D}(x)^{n-1}f(x).$ I'm interested in the properties of the following "expected density" (I'm not sure there's a better way to put it) of the nth order PDF: $\displaystyle\int_{-\infty}^{+\infty}(n-1)F_{D}(x)^{n-2}f_{D}(x)\times f_{D}(x)\:dx$, which simplifies to $\displaystyle\int_{-\infty}^{+\infty}(n-1)F_{D}(x)^{n-2}f_{D}(x)^2\:dx$ What I want to show is that this expression decreases as $n$ increases. **_What I have gotten so far:_** By an integration by parts trick, we can show that the above can be expressed as: $\int_{-\infty}^{\infty}f_{D}(x)dF_{D}(x)^{n-1}$ $=f_{D}(x)F_{D}(x)^{n-1}|_{-\infty}^{+\infty}-\int_{-\infty}^{+\infty}F_{D}(x)^{n-1}df_{D}(x)$ $=-\int_{-\infty}^{+\infty}F_{D}(x)^{n-1}f^{\prime}_{D}(x)dx.$ Intuitively, I can see that as $-f^{\prime}_{D}(x)$ is positive on $x>0$ and negative on $x<0$. And as $n$ increases $F^n_{D}(.)$ shifts more mass to extreme values where $f^{\prime}_{D}(x)$ is very close to 0. So, eventually the whole integral will become smaller as $n$ increases. But I'm not sure how to proceed this argument formally. Any help will be greatly appreciated!


title: Histogram of uncertain data

text: I have a set of values, each of these values has its own mean and variance. I want to be able to account for this variance when I plot the histogram of the means. Something like an error bar on the bins. Is there a commonly used method to this?


title: Probability for completely unfair stats when rolling ability scores using 4d6 drop lowest for D&D

text: Players of a certain TRPG have characters with 6 ability scores, each ability score ranging from 3-18. One method of generating those is by rolling 4d6 drop lowest. That means four six-faced-dice are rolled, and the three highest results are added. What's the probability that, given 5 players, one player will have a highest ability score equal or lower than the lowest ability score of another player? The related question here shows how to get the distribution of 4d3 drop lowest, but how do I get from there to an answer of my question above? A good answer would explain the result in a way that a statistics novice can follow.


title: Statistical Blockade - Estimating Rare Event Statitistics

text: In the world of microelectronic design and simulation where the technology gets smaller and smaller (a few nanometers nowadays), the high density and replication of some of our circuits require us to estimate rare event statistics in our design flows. Traditional Monte Carlo cannot be applied as it would take years to simulate billions of samples. Therefore, We use some statistical techniques like Response Surface Modeling, Worst Case Distance (I think some call them Most Probable Point). Recently, There have been talks about a method called "Statistical Blockade". Searching this method would often point you to the same papers from Amith Singhee and Rob A. Runtenbar from IBM. Here is one of them: Statistical Blockade: A Novel Method for Very Fast Monte Carlo Simulation of Rare Circuit Events, and its Application **Question** : Without asking you to spending hours on this paper, I am wondering on whether similar extreme-value approaches are known/used in other disciplines.


title: Multi-variate normal distribution distance from vector sub-space

text: let X $\sim {\cal N}(\mu,C)$ be a random variable obeying multi-variate normal distribution in $\mathbb{R}^n$ and $U \subset \mathbb{R}^n$ be a vector space with $\dim(U)=n-1$. What is the probability of X to be at distance ($L_2$) d from U? Assume that d is small (I actually interested in $\lim\limits_{d \rightarrow 0} \frac{P(d)}{d}$ where P(d) is the above probability).


title: Benefits on kernel trick and one related question

text: My questions are:   1. What are the benefits of kernel trick? Can anyone summarize?   2. One thing I read about in one lecture note is , ' _we never need to explicitly represent feature vectors_ '. I do not quite understand what ' _explicitly represent_ ' means. Can any one explain? I saw it here, at page 16. PS: if it is appropriate, please add the tag _kernel trick_.


title: How to reduce residual error?

text: I have a time series as following. It has an upward trend and weak sesonality (found in ACF). I tried spectral analysis, but the residual error kept staying over 600, which is not small enough for a prediction. Any advice? Thx!               ts <- c(1530.040 3527.839 3327.613 3613.749 2971.990 3747.925 3352.686 5203.718 4550.774 3497.677 2719.688 2618.564 1378.669 3743.931 5094.959 3998.374 3649.445 3709.277 4078.814 6663.439 5774.322 4881.767     4233.501 4824.077 3207.395 6078.317 5398.443 3541.734 2852.284 3253.765 4716.344 5239.078 5267.733 4072.201 3796.414 4309.636 3751.775 5954.923 5627.322 3862.649 3673.138 4551.715 5904.327 5712.579 6108.412 4706.165 5230.390 3756.871 3079.763 5934.894 5172.300 3194.798 3167.828 4333.597 4512.727 7609.356 7728.995 6024.994 4947.637 4491.105)


title: How to get the correlation between two large (sparse) matrices?

text: I have two large (sparse) matrices (500000*500000) and would like to perform the mantel test to get a similarity measure.   * Is there a way to break-up these matrices into smaller ones and perform the test?    * Or can you suggest a way to do it efficiently without running into memory problems?


title: What are essential rules for designing and producing plots?

text: ### Background to Question: Previously on Cross Validated, we have had questions on:   * What is best practice when preparing plots?   * What are good tips available online for plotting two numeric variables? It was suggested by @david in the comments to this question that we should have a **community wiki** question with **one visualization rule per answer** that the community could vote on. ### Question   * What are the essential rules on designing and producing graphical representations of data? ### Rules   * One rule per answer   * Ideally, include a brief explanation of why you think it is a good idea   * Answers with examples (code and image) of good and bad practice preferred. _Feel free to edit the question if you think it could be refined._


title: Semi-supervised learning, active learning and deep learning for classification

text: Final edit with all resources updated: For a project, I am applying machine learning algorithms for classification. **Challenge:** Quite limited labeled data and much more unlabeled data. **Goals:**   1. Apply semi-supervised classification   2. Apply a somehow semi-supervised labeling process (known as active learning) I've found a lot of information from research papers, like applying EM, Transductive SVM or S3VM (Semi Supervised SVM), or somehow using LDA, etc. Even there are few books on this topic. **Question:** Where are the implementations and practical sources? * * * **Final update (based on helps provided by mpiktas, bayer, and Dikran Marsupial)** **Semi-supervised learning:**   * TSVM: in SVMligth and SVMlin.    * EM Naive Bayes in Python   * EM in LinePipe project **Active learning:**   * **Dualist**: an implementation of active learning with source code on text classification   * This webpage serves a wonderful overview of active learning.   * An experimental Design workshop: here. **Deep learning:**   * Introductory video at here.   * General site.   * Stanford Unsupervised Feature Learning and Deep Learning tutorial.


title: How to determine larger probability of drawing one element from a set when set size and number of draws with replacement vary?

text: Which is the better scenario in terms of probability with the following: **Scenario 1**   * I have a hat with 5 tickets in it (I have one).   * There is one single draw.   * I have a 20% chance that my ticket is drawn. **Scenario 2**   * I have a hat with 25 tickets in it (I have one).   * There are 5 draws.   * After each draw the ticket is returned before the next draw. i.e. potential chance of being drawn 5 times.   * I have a 4% chance each draw Which has the better probability of being drawn at least once? Ideally I'd like to know how to work this out (and why) so I can learn something :-)


title: How to calculate overall rankings/ratings from many pairwise comparisons?

text: I have a number of categories and a bunch of data about pairwise comparisons between those two categories; how do I learn information about overall rankings/ratings/etc. of all the categories? Let me provide more detail. I'm trying to learn what the general population is most worried about, by surveying them. I have about 20 categories of potential concerns, call them $A, B, ..., T$. I have picked a random sample of people. For each person in the sample I have done the following. I picked a pair of concerns, e.g., $E$ and $Q$. Then I asked the respondent: which are you more worried about, $E$ or $Q$? They could select one of three possible answers: $E$, $Q$, or "about the same for both". I have about 300 responses. Each respondent got a different, independently random selection of a pair of concerns. How should I analyze this data, to extract as much information as I can? For instance, can I learn anything about what concerns people are generally most worried about, or least worried about? Can I compute some approximate rankings or ratings or some comparisons? And, of course, I am concerned about computing statistical significance for any conclusions drawn. Any advice? P.S. Feel free to choose any reasonable model for how people will respond. For instance, you could assume transitivity for each individual (each person has an internal ranking of concerns, and answers the question based on that ranking; so if they would answer that they're more worried about $P$ than $Q$ if asked about $P$ vs $Q$, and would say they're more worried about $Q$ than $R$ if asked about $Q$ vs $R$, then if asked about $P$ vs $R$ they would say they're more worried about $P$ than $R$).


title: Marginal distribtuion with noninformative prior

text: For $i=1, \ldots, K$ and $j=1, \ldots,n$, assume the following model. \begin{align} X_{ij} \mid \mu_i & \sim N(\mu_i, \sigma^2) \nonumber \\\ P(\mu_i, \sigma^2) & \propto 1/\sigma^2 \end{align} We assume everything is independent. Then the posterior is as below \begin{align*} \dfrac{\mu-\bar X_i}{s/ \sqrt{n}} \sim T_{n-1} \end{align*} With this noninformative prior, is it possible to derive the marginal distribution of $\bar X_i$? Thank you very much!


title: Classifier weighted towards recall?

text: I have a classification problem where getting true positives is much more important than true negatives. To be clear, I know that roughly 10% of my population are actual positives, but I can assign some proportion (say 30%) of the population to be classified as positives without much cost, above all else I need to make sure that the actual positives are covered by this set. Unfortunately the classification tools I am using in WEKA seem to be balancing precision and recall, so that it is a) not assigning as many positives as it is allowed to, and b) getting quite a bad recall value. Is there a standard way of approaching this problem? My first guess would be to weight the cost function towards recall rather than F-score, but I don't see an easy way to do this in WEKA.


title: Best fitness metric for long term data

text: I'm trying to understand monthly data going back about 100 years. As you can imagine old data values can be very small, \$1 compared to recent values of \$100. What would be the best fitness metric to use as a forecast, based on the notion that being off \$1 100 years ago is a huge deal but \$1 now is not a big deal? Some of the options in my software are MAS, MSE, RMSE etc. Has someone done a visual comparison of fitness metrics online?


title: Looking for help identifing outliers in a pilot study to guide future hypothesis testing

text: I am looking to study a particular type of error on a cognitive test to evaluate for potential clinical implications. As there is no existing research on this variable, I would like to run a pilot evaluation on a relatively large pool of subject data from our clinic (feasibly this would be about 200 subjects) to see if there are any trends to guide hypothesis testing in future studies (e.g., are patients from a particular clinical population more likely to make this type of error than those from another population). The data is such that a patient could potentially make 0-15 of these errors, although it's likely that there will be a strong tendency toward the 0-5 range. Essentially, I am looking to identify subjects from within this pool who are making a relatively higher rate of these errors. Is there a recognized format for doing this type of pilot work? Barring that, what would your recommendations be? This is obviously very exploratory, so I feel like I would have some latitude in defining what constitutes an 'outlier' for the purposes of this study, but any references or suggestions would be much appreciated.


title: Multiplying pdf's: Statistics

text: If I have two different independent continuous variables, X and Y. Is it legitimate to multiply their pdf's to form a pdf for the ordered pairs (X,Y) as f(x,y). I've plotted the surface and the volume under it is 1, and the key points (peaks and troughs seem to be in the correct position). Is this a legitimate process? What issues do I need to be aware of?


title: Finite Population Correction Factor

text: 1. A simple random sample of 21 observations is drawn from a finite population of 41 individuals. If the population mean is 300 and the population variance is 7, then after applying the finite population correction factor the variance of the sample mean is: a. 0.1667 b. 0.3333 c. 0.5774 d. 7.0000 Using my formula: corrected variance=(²/n)((N-n)/(n-1)) I keep getting (b), but the key says the answer is (a).


title: Definition of Principal component

text: Suppose $u$ is the vector that maximises the variance of the projection of the data with design matrix $X$. Now, I have seen materials that refer $u$ as the (first) principal component of the data, which is also the eigenvector with largest eigenvalue. However, I have also seen that the principal component of the data is $X u$. Obviously, $u$ and and $Xu$ are different things. Can anyone help me out here and tell me what is the difference between these 2 definitions of principal components?


title: calculating conditional probability-bayes rule

text: I'd like to calculate the conditional probability in the following case: I was told that a box contains a BLUE ball. this is my evidence, my prior probability of a BLUE being drawn is 0.3. and this message can come from a friend(0.6) or an enemy (0.4) I also know the probabilities of: a friend observing BLUE : 0.16 a friend observing RED: 0.09 an enemy observing BLUE: 0.05 an enemy observing RED: 0.17 given that the ball is actually BLUE. I would like to know the conditional probability of the ball actually being BLUE given this information.


title: Looking for redacted text corpus

text: I am looking around for datasets mainly in the form of natural language text corpus, that has been redacted by experts. So far I found only the enron dataset, but it appears only a few specific things (like employee mail-id, etc.) were redacted, which will not be beneficial for feature extraction. Something like declassified government records or company data would be awesome; medical records would probably be comparatively easier to get, but are quite unrelated to the problem. Any idea if such dataset is available?


title: Estimating population average models in lmer or geepack

text: I would like to know how to estimate a population average model of a hierarchical logistic regression using `R` package `geepack`. The `Stata` code is:               xtlogit dep ind1 ind2 ind3, i(ind4) pa      I would like to reproduce this in `R` using `geepack` or any other method.


title: Sampling random numbers from a distribution with asymmetric confidence intervals generated by a bootstrapped estimate

text: So, I'm working on a problem where I have a sample distribution with a reported, mean and an asymmetric confidence interval from bootstrapped sampling of the initial population (this is from a meta-analysis - so, I also have the # of bootstrap replicates and the original sample size). Given this information, I would like to draw a random variable (I'm then going use it to construct my own bootstrapped sample using other similar variables - basically, a meta-analysis of meta-analyses). I'm sure this is simple, but I'm completely blanking on how one would do this - is there an appropriate distribution to draw from or way to back-engineer reported results to create one? R code would be great (if there is a package to do this), or just a see- this-reference, as I'm drawing a blank. N.B. I think my larger question here was, if you have the quantiles for a population, but no other distributional knowledge or data, is it possible to create a method to draw random numbers based on that information. But perhaps the answer there is no.


title: Random geometric graph connectivity

text: Given N random points uniformly distributed in the unit square, and a distance d, i can generate a matrix in the following format:                            E V1 V2     [1,] 0.5564821  1  2       [2,] 0.3373116  1  3     [3,] 0.3973278  1  4     [4,] 0.6066518  1  5     [5,] 0.9603731  1  6     [6,] 0.3612895  1  7     # more rows...      Where E is the edge linking vertex V1 and V2. I've just started learning graph-theory, so i ask: How can i determine if this random geometric graph is connected? Thanks!


title: Regression: what to do with categorical predictors?

text: I've got a few categorical predictors (like gender,...) and now I want to build regression models. So I've made the categorical predictors numeric by for example: "female" --> 1 and "male" --> 0. But when I do methods like nearest neighbors regression I have to standardize all the predictors (for example the weights). What to do here with the categorical variabels (that were made numeric)? Does this also have to be standardised? This seems so weird. Silke


title: Time Series forecasting with Incident Data: If I have a date field, do I need to recode into summary tables in R?

text: So I have a database of incidents, each identified by a unique number. Each individual incident is one row, and contains the date it occured. Data is roughly in the format:               #   Date   State   Pop.   Unemp. Rate   CategoricalVar1   Cat.Var2   Cat.Var3   Cat.Var4      So, I was going to do clustering for incidents, multinomial logistic regression for the categorical variables, based on quant and qual input variables with testing for IIA, C&RT similarly to the multinomial log (ignoring for now the domain expertise-based question of, "are the other factors predictive or causative of each categorical variable?"). There's no requirement for a time series component here. The goal abstractly is, given population/unemployment rate/all but "Y" categorical variables, can we predict the remaining categorical variable? Let's just accept the premise that this is how the logical law enforcement process would flow, and that the remaining variable is something of which awareness beforehand can be useful and actionable/preventable. Another different area of research would be by state and over time, trying to predict the number of a certain type of crime or location (as in, "office building/gas station/shopping mall/etc") of crime in future years. I have 21 years of data, 160k incidents. This obviously has a time series component. I'm working in R, and am messing around Stack concurrently, but it occurs to me the real issue here is how I'm conceptually tackling the TS framing with the data in the given format. The central question is: can I do TS as clumsily described above (or perhaps, in some other way that is obviously more logical/useful?) with incident data, where date is an X variable, or would it be easier to reshape the data into the following format, for each categorical variable:               CatVar1:     Date   State   Tot # of Incidents w/CatVar1 = 1   Tot w/CatVar1 = 2   etc.     1997     VA                15                             8            5      IE, does regrouping my ungrouped data pose any advantage? When would I want to regroup and when wouldn't I? I'm sorry, I know this is a beginner's question. I guess xts/zoo/forecast handles this? If we can concentrate here on whether my conceptual framing is wrong more than anything else, it would be most helpful to me. Thank you!


title: Mixing and dividing point processes

text: At the following figure at left side two realizations of point processes with different density (intensity) $\lambda_1$ and $\lambda_2$ is being mixed matching the center of the belonging areas to build a point process in the middle with intensity $\lambda$. Then randomly selected points as two sets extracted from it as shown at the right side.   Questions:   Is $\lambda=\lambda_1+\lambda_2$? and Is $\lambda=\lambda_3+\lambda_4$?   If two at left side were Poisson PP, Is the middle one a Poisson PP?   How about the two at right hand side? ![enter image description here](http://i.stack.imgur.com/aAtiG.jpg)


title: Determining the influential features for an outcome

text: I have a small table like this               Outcome     Feature_1_Value  Feature_2_Value   Feature_3_Value   .....        1               1              1.2               3.7        9.7             2.3            3.976             4.65        ....            ....           ....              ....      There is an outcome that depends on what values do a set of features have. The above table represents the training data. The outcomes and the features all are continuous values and are not discrete ( it is not 0/1/2 kind of values). Now, how do I analyze the training set of data to determine which of the features actually influence the output and how do I re-arrange the features in decreasing order of their influence? I would like to represent nicely to the end-user saying Feature_n_Value was twice more influencing the output than Feature_m_Value, but Feature_n_value was less than half influential for the output compared to Feature_p_Value, and so on. What is the mathematics behind such a calculation. I looked a bit at logistic regression. Would that be a right way to do? If so, then in logistic regression, I get a set of beta vectors (num of output values), each of the size of the dimension of the input (num of features). What do I do with the output beta vectors, though I know that you can predict an output given a fresh set of feature values and these beta values that are calculated using training data. Advice is highly appreciated. Thanks Abhishek S


title: Compute the variance of parameter estimates given limited number of samples

text: I'd like to infer the variance of estimated parameter $\hat\theta$ of the density function of $f(x;\theta)$ given only a limited number of samples $X_1,\cdots,X_n$. Bootstrapping doesn't perform well when the sample size is small. Is there any other methods to deal with the case?


title: Statistical test for multiple conditions and features with replicates

text: I have data from a biological experiment of the following form:                                          site1   site2   site3   site4     Condition A | replicate1 |  0.30 |  0.15 |  0.50 |  0.05     Condition A | replicate2 |  0.27 |  0.18 |  0.48 |  0.07     Condition B | replicate1 |  0.50 |  0.15 |  0.30 |  0.05     Condition B | replicate2 |  0.47 |  0.18 |  0.38 |  0.07      We have two conditions A and B, and two replicates for each condition, giving in total 4 samples. For each of these four samples we measure the usage (how many sequencing reads) at four sites. What I want to know is if the usage of these four sites is significantly changes between condition A and B. Any ideas for a clever way how to do this? Basically, the variation within a condition should be significantly lower than the variation between the two conditions. Yes, I just used the fraction of sequencing reads up there. The reason is that these are four separate experiments and cannot simply use the raw counts without any normalization. If I just look at the relative usage of each site in each experiment I can avoid the normalization step. I can scale those numbers to 100 and then I have integers if this helps. Is this an example of a classical ANOVA? I tried it in R, but isn't the sample size too low to get something meaningful?                  condition site1 site2 site3 site4     1         A    30    15    50     5     2         A    27    18    48     7     3         B    50    15    30     5     4         B    47    18    38     7          summary(aov(site1+site2+site3+site4~condition, data=data))


I recently performed an experiment to implement Gamification (a simple points system) into a web application to see if the points had an effect on the users in terms of usage. We have data on how much time users spend daily in the application, 1 month before Gamification was implemented, and 1 month after. The users selected for the experiment were sampled into 3 different groups:   1. Control group (users who did not use gamification) - 300 users   2. Points group (users who used the points system) - 300 users, and   3. Performance group (top 10% points earners) - 156 users **What statistical tools can I use to evaluate this data?** All I have done now is to look at how the total time for each group changed from the first month to the next, and got:   1. Control group: +1%   2. Points group: +28%   3. Performance group: +49% However, this is after I've removed what seemed to be outliers (new users appearing in the middle of the experiment etc.). Prior to this, it looked like this:   1. Control group: +41%   2. Points group: +105%   3. Performance group: +110% Is this enough proof that the points did, in fact, have an impact on usage?


title: Interpretation of interaction term in Cox PH model when centering LP on mean values of predictors

text: Let's say I have a Cox PH model for predicting the risk of dying, that in a simplified form looks something like this:               Absolute risk = 1 – baseline hazard^(exp[LP-LPmean]), where          LP (=linear predictor) = 0.54815*age - 0.06318*age [if on treatment]                               - 0.06351 [if using treatment]     baseline hazard = 0.98123     LPmean (=linear predictor for patient with mean values for covariables) = -0.03177       The model thus contains a coefficient for age and treatment and an interaction between age and treatment. The LP was centered at the mean values of the predictors, hence, the LP of the mean predictor values is subtracted from the LP in order to get predictions. I would like to get an idea about the effect of the interaction in the model and give a hazard ratio (HR) for treatment at a few different ages. If I do that as following:               for age 50: HR treatment = exp(-0.06318*(50) - 0.06351) = 0.04;     for age 60: HR treatment = exp(-0.06318*(60) - 0.06351) = 0.02;     etc..,      .. the effects are rather extreme (due to the centering of the LP). I can account for the centering like this:               for age 50: HR treatment = exp(-0.06318*(50-54.6) - 0.06351) = 1.25;     for age 60: HR treatment = exp(-0.06318*(60-54.6) - 0.06351) = 0.67;     etc..,      where 54.6 is the mean age in the study population. However, I am not sure whether this is correct and how I should interpret this. Many thanks!


title: Calibrating random forest regression

text: I'm trying to learn how to calibrate random forests for regression. I have found many informative sources on how to do this for classification but none for regression. The CORElearn R package has a random forest classifier with local weighting of basic models. But nothing for regression. I am basically trying to correct the bias random forests present when attacking a regression problem. A problem that I have seen worries many people but I have not seen explicit solutions to yet. Any new thoughts on this?


title: Euclidean distance score and similarity

text: I'm just working with the book Collective Intelligence (by Toby Segaran) and came across the Euclidean distance score. In the book the author shows how to calculate the similarity between two recommendation arrays (i.e. $\textrm{person} \times \textrm{movie} \mapsto \textrm{score})$ . He calculates the Euclidean distance for two persons $p_1$ and $p_2$ by $$d(p_1, p_2) = \sqrt{\sum_{i~\in~\textrm{item}} (s_{p_1} - s_{p_2})^2} $$ This makes completely sense to me. What I don't really understand is why he calculates at the end the following to get a "distance based similarity": $$ \frac{1}{1 + d(p_1, p_2)} $$ So, I somehow get that this must be the conversion from a distance to a similarity (right?). But why does the formular looks like this? Can someone explain that?


title: What are disadvantages of using the lasso for variable selection for regression?

text: From what I know, using lasso for variable selection handles the problem of correlated inputs. Also, since it is equivalent to Least Angle Regression, it is not slow computationally. However, many people (for example people I know doing bio-statistics) still seem to favour stepwise or stagewise variable selection. Are there any practical disadvantages of using the lasso that makes it unfavourable?


title: Logistic regression control variables vs. layer in crosstabs

text: I would like to know how SPSS controls for variables in logistic regression in order to report the correct "controlled" OR. The specific case which bothers my is as follows: Let's say I have one IV tumor size, which is either small tumor size or large tumor size. The other IV is risk classification with low risk and high risk. My dependent variable is complications: yes or no. Unfortunately people with a high risk often also have a large tumor. To assess the actual effect of the risk classification, I used the layer tumor size and did a crosstab, which gave me the OR for high vs. low risk patients when they have large tumors and high vs. low risk patients when they have small tumors. I assumed this was the controlling of variables. Subsequently I heard that, when one performs a logistic regression in SPSS under table "variables in the equation" the OR $(\exp(B))$ for each variable while controlling for the other variables is displayed. Unfortunately these ORs do not match with the ORs of the crosstab version. What is the difference in that case? Thanks a lot, I appreciate your help!


title: How difficult is it to train a gaussian mixture model compared to other models?

text: I have finally been able to wrap my head around the mechanics of how to initialize and train a multivariate Gaussian mixture model using expectation maximization algorithm. So I wonder how difficult this GMM and EM task is in comparison to all other **common** algorithms and models in machine learning. I appreciate any feedback. Thank you.


title: Gaussian kernel estimator as Nadaraya-Watson estimator?

text: I'm working on a problem from _"The Elements of Statistical Learning"_ (prob. 6.8): > Suppose that for continuous response $Y$ and predictor $X$, we model the > joint density of $X, Y$ using a multivariate Gaussian kernel estimator. Note > that the kernel in this case would be the product kernel $\phi_{\lambda}(X) > \phi_{\lambda}(Y)$. > > (a) Show that the conditional mean $E(Y|X)$ derived from this estimate is a > Nadaraya-Watson estimator. > > (b) Extend this result to classification by providing a suitable kernel for > the estimation of the joint distribution of a continuous $X$ and discrete > $Y$. I know that the Nadaraya-Watson estimator is just the weighted average (equation 2.41 and 6.2 in ESL): > $$\hat f (x_0) = \frac{\sum_{i=0}^N K_{\lambda}(x_0, x_i) y_i}{\sum_{i=0}^N > K_{\lambda}(x_0, x_i)}$$ Where $K$ in this case would be the multivariate Gaussian kernel function. I can think about how to extend this to a classification problem, but am not sure how to approach the first part of this question. Any pointers would be greatly appreciated!


title: Categorical or continuous scale for area chart?

text: I would like to plot the attached income distribution dataset (rendered as an image) as an area chart. As you can see, personal income is divided into 26 intervals of varying width. I also have the average and mean income in the intervals. To convey a truthful area graphic of this data, I wonder what my options really are?   Plotting the ordinal categorical data at hand would yield a big hump in the area chart for the 400-499 interval. But this is only because that interval is wider and the user can hence be misguided by the shape. Another issue with the categorical data is that the average of the "1000+ interval" is very far from 1000 (= 1644). An area graphic not taking this into account would do a bad job in showing the actual distribution. How would you go about and is there any way in which I can use the average/mean to "convert the categorical scale to a continuous scale"?


title: Coverage probabilities of the basic bootstrap confidence Interval

text: I have the following question for a course I'm working on: > Conduct a Monte Carlo study to estimate the coverage probabilities of the > standard normal bootstrap confidence interval and the basic bootstrap > confidence interval. Sample from a normal population and check the empirical > coverage rates for the sample mean. Coverage probabilities for the standard normal bootstrap CI are easy:               n = 1000;     alpha = c(0.025, 0.975);     x = rnorm(n, 0, 1);     mu = mean(x);     sqrt.n = sqrt(n);          LNorm = numeric(B);     UNorm = numeric(B);          for(j in 1:B)     {         smpl = x[sample(1:n, size = n, replace = TRUE)];         xbar = mean(smpl);         s = sd(smpl);              LNorm[j] = xbar + qnorm(alpha[1]) * (s / sqrt.n);         UNorm[j] = xbar + qnorm(alpha[2]) * (s / sqrt.n);     }          mean(LNorm < 0 & UNorm > 0); # Approximates to 0.95      From what I've been taught for this course, the basic bootstrap _confidence interval_ can be calculated like this:               # Using x from previous...     R = boot(data = x, R=1000, statistic = function(x, i){ mean(x[i]); });     result = 2 * mu - quantile(R$t, alpha, type=1);      That makes sense. What I don't understand is how to calculate _coverage probabilities_ for the basic bootstrap CI. I understand that the coverage probability would represent the number of times that the CI contains the true value (in this case `mu`). Do I simply run the `boot` function many times? How can I approach this question differently?


title: Why do Excel and WolframAlpha give different values for skewness

text: for the following 3 values 222,1122,45444 WolframAlpha gives 0.706 Excel, using `=SKEW(222,1122,45444)` gives 1.729 What explains the difference?


title: Calculate t without data set

text: How would you go about calculating test statistic t from a regression equation and the standard error of the slope?


title: How to compare similarity between two discrete distributions

text: I am trying to compare two temporal distributions. To give an example, consider two factories. For each factory we have a distribution of number of hours taken to produce a product. E.g. 40% of total products were produced in less than an hour, 20% in less than 2 hours but more than 1 and so on and so forth. What is a good way to measure similarity between the two distributions arising from two different factories. I considered using earth movers distance but its a distance metric and not a test. Thanks


title: How to best display results from a factorial design?

text: How can I create plots in R like this mockup made in Powerpoint? ![image of a grid](http://i276.photobucket.com/albums/kk30/baixiwei/grid.png) In other words, for each combination of x-axis and y-axis values, I want to plot two triangles which, together, form a square. I have two vectors of values, one of which determines the color of the lower triangle, and the other of which determines the color of the upper triangle. I can't figure out how to plot triangles which will fit together the way I want. Here's sample code:               # generate sample data     D = data.frame( expand.grid( x=seq(1,10,1), y=seq(1,10,1) ) )     D$color1 = sample( c('red','blue'), nrow(D), replace=TRUE )     D$color2 = sample( c('red','blue'), nrow(D), replace=TRUE )          # this will plot just one color in each cell, not two as I want:     plot( D$x, D$y, col=D$color1, pch=15, cex=4 )      In case it matters - the point of the plot is to show the results of simulated experiments for each combination of values of parameters x and y. Thus, each box represents the outcome of one simulation. It's a 2x2 experimental design with factors I'll call A and B, each with 2 levels, A1/A2 and B1/B2. So red color means A1>A2 on the outcome measure, and blue means A2>A1. Above the diagonal represents the result for B=B1 and below the diagonal represents the result for B=B2. Note, by the way, that I need to make several such plots, one for each value of a third parameter I'll call z. I'm open to suggestions if someone thinks there is a better way to visualize this data.


title: Modeling frequency over time

text: I have a two-dimensional data set that looks like $(t, x)$ where $t$ is a time in seconds when event $X$ happened. $X$ ranges from $[0, 200]$. I want to visualize the frequency of each $x$ at time $t$ over some time period. I guess this would be a bar graph with $x$-axis being event #, $y$-axis being frequency, and $z$-axis being time, $t$. Furthermore, I would like to group all events that happen within say a 5 second interval to count towards the same frequency bar on the $y$-axis. If there is a way to do this with R that would be even better. My goal is to get a sense how often some event occurs over the course of a day, and when certain events happen a lot or infrequently. If you know of a better way to understand this information, I am all ears.


Why is the design effect in most sample studies taken as 1.25? Who and how was it calculated first of all?


title: Pearson correlation coefficient on multiple parameters

text: First of all: I'm mathematically challenged. Question 1: I use the Pearson correlation coefficient to correlate people. As a real life example: I correlate those people by comparing the number of points they've given to movies. This works very well. Now I would like to correlate not 1 variable but multiple, e.g. correlate how much they like the individual actors. How should I do that? Just take the average of all n correlations? Question 2: if person $A$ and $B$ correlate for 0.9 than that is better then $A$ and $C$ who correlate for 0.8. But what if I had e.g. 25 movies while comparing $A$ and $B$ and 100 while comparing $A$ and $C$, should I then still consider $A$ and $B$? Or is there some extra weighting formula for this?


title: Bayesian logistic regression with mixed (two level) data

text: I have two datasets, the first on schools, and the second lists students in each school _who have failed in standardized test_ (emphasis intentional). The first dataset has the structure below (I'm sorry I could not figure out how to insert a proper data snippet):               School_ID   Total_White   Total_Black   Total_Asian   School_Revenue      The second looks like this:               Student_ID   School_ID   Race      I am trying to estimate the probability of failure, given student race and school revenue. If I run a multinomial discrete choice model on the second dataset, I shall clearly be estimating `P(Race | Fail=1)`. I obviously have to estimate the inverse, `P(Fail=1 | Race)`, by merging the two datasets by School_ID. Since all the pieces of information are available in the two datasets (`P(Fail), P(Race), Revenue`), I see no reason why this can't be done. But I am stumped as to actually how to implement in R. Any pointer would be much appreciated. Thanks.


title: Regression for rank order dependent variable

text: I am looking for a package to analyze rank-order data in the form of regression. My outcome (dependent) variable is a rank ordering in the sense that the respondents to a questionnaire were asked a single question: to assign the values 1 to 4 to four different objects (e.g. to rank four cars as best , second best etc). Their first preference was number 1, their second was number 2 etc. So, for each respondent we have four observations. My first reaction was to analyze the data as a repeated-measures design where I considered the four rankings as four observations for the same person (one response per option). This would be an easy case of ordinal regression (e.g. proportional odds model) using the excellent package of Frank Harrel (thank you Frank). However, these are not four independent observations per person, because if a respondent assigned the values to the first three objects, then the rank of the fourth object is pre-determined. Does anyone know if there is a package that can fit a more appropriate model such as the Sin‐Ho Jung and Zhiliang Ying (2003) "Rank‐based regression with repeated measurements data", Biometrica, 90 (3), 732-740.


title: What are some fellowships available for US students studying Statistics?

text: Hello Stats community, I wanted to ask if anyone has advice on finding graduate student fellowships for students in the US studying statistics, data mining, or data analysis? I know the NSF has a 3 year long fellowship, as well as the DoE


title: Bound on moment generating function

text: This question arises from the one asked here about a bound on moment generating functions (MGFs). Suppose $X$ is a bounded zero-mean random variable taking on values in $[-\sigma, \sigma]$ and let $G(t) = E[e^{tX}]$ be its MGF. From a bound used in a proof of Hoeffding's Inequality, we have that $$G(t) = E[e^{tX}] \leq e^{\sigma^2t^2/2}$$ where the right side is recognizable as the MGF of a zero-mean normal random variable with standard deviation $\sigma$. Now, the standard deviation of $X$ can be no larger than $\sigma$, with the maximum value occurring when $X$ is a discrete random variable such that $P\\{X = \sigma\\} = P\\{X = -\sigma\\} = \frac{1}{2}$. So, the bound referred to can be thought of as saying that the MGF of a zero-mean bounded random variable $X$ is bounded above by the MGF of a zero-mean normal random variable whose standard deviation equals the maximum possible standard deviation that $X$ can have. My question is: is this a well-known result of independent interest that is used in places other than in the proof of Hoeffding's Inequality, and if so, is it also known to extend to random variables with nonzero means? The result that prompts this question allows asymetric range $[a,b]$ for $X$ with $a < 0 < b$ but does insist on $E[X] = 0$. The bound is $$G(t) \leq e^{t^2(b-a)^2/8} = e^{t^2\sigma_{max}^2/2}$$ where $\sigma_{\max} = (b-a)/2$ is the maximum standard deviation possible for a random variable with values restricted to $[a,b]$, but this maximum is not attained by zero-mean random variables unless $b = -a$.


title: Composite variable from Likert items

text: I have six 5-point Likert items (1-Strongly Disagree, 2-Disagree, 3-Neither agree or disagree, 4-Agree, 5-Strongly Agree) which are related to my dependent variable. I sum the six variables which gives the frequency as follows:                   Frequency     6.00    1     9.00    1     12.00   6     13.00   3     14.00   3     15.00   18     16.00   11     17.00   4     18.00   17     19.00   11     20.00   14     21.00   22     22.00   15     23.00   9     24.00   11     25.00   2     26.00   3     30.00   2     Total   153      Now I want to recode the sum variable again into a 5-point scale (1-5) to get a composite score of my dependent variable. I will use this composite dependent variable in further analysis of ordinal logistic regression. My question is this: What is the best way to recode the sum again into 5-point scale? Or, please let me know if there exists any better way to create a composite scale from Likert items? just to clarify my quesion: I have six variables that are related to my DV, all are on 5 point likert scale. and I could not find the best way to reduce these six variables into one. The composite score may be considered as one of the option.


title: How to report general precision in estimating correlations within a context of justifying sample size?

text: I have a study with 130 participants where most of my analyses involve examining the relative size of correlations between various psychological variables. I'm trying to comply with recommendations that request authors to say something about the statistical power of their study. Statistical power focuses on whether the correlation is significantly different to zero, but the null hypothesis is not particularly interesting in my research. Everything is correlated with everything. It's the relative amount that matters. I'd like to say something in the method section of my report about the precision of estimating the true correlations.   * Is it a good approach to present the 95% confidence intervals on correlations when the correlation is zero, or perhaps some other focal value?   * Or is there a better approach to indicating the precision of estimating correlations in a given sample?


title: When controlling a binomial proportion, how to deal with proportions with low confidence?

text: I'm building a system to alert when the proportion of successes to attempts goes out of control. To visualize the data, I've produced some rolling control charts to get an idea of when things would have alerted. A sample of one such chart is here, roughly labeled in photoshop for some clarification: ![enter image description here](http://i.stack.imgur.com/qgOuX.jpg) Since each interval (day) has a widely varying amount of attempts, I'm concerned that I am going to have outliers occur due to sparse data (and low confidence) and that trigger an alarm. I'd like to only trigger an alarm when I have enough data to be confident in the generated proportions, but at the same time, leverage as much data as is available. Is it a good idea to incorporate confidence intervals in such cases, and only trigger if I am, say, 75% confident that the point is out of the lower control limit? It seems like if I have several points in a row that are out of the lower control limit, each with 50% confidence then it's a 75% chance that one of them is and I should alert. Am I off on a wild track here? Has this ever been suggested as a best practice? I am new to control charts and have just ordered Donald Wheeler's book. On a side note, I am beginning to question my decision to use geometric means and standard deviation as the chart doesn't ever appear to swing up as dramatically as it does down (geometrically speaking).


title: Boosting ranked lists

text: Suppose 10 people conduct a 10 studies which rank people's top five favorite colors. For instance, one survey may produce a list 1) red 2) blue 3) green 4) yellow 5) orange, where red is the favorite color on the list and orange is the least favorite on the list. The methodology of the 10 surveys are completely different, and each one has its own weakness. For instance, one survey may look at the color of people's clothes and assume that people like to wear their favorite color, while another survey may ask people to tell them their favorite color without asking any more questions about their second, third, forth, or fifth favorite color. Is there a known good "boosting" algorithm for combining the 10 surveys to get one very good ranking that combines the strengths of each survey?


title: Reporting $\chi^2$ test results in APA format

text: I'm trying to report the results of my chi-square test for independence in APA format but these results don't resemble the examples I've seen especially the 4 digits before the period in the chi-square test results i.e. 9355.19 and 9556.44. Am I doing this right? > "the χ2 tests showed significant results in the association of this > knowledge with years professional of experience χ2 (4, N = 302) = 9355.19, p > = 0.83 and with professional level χ2 (4, N = 302) = 9556.44, p = 0.30." **These are the values:**                        I agree I disagree I don't know |Total     Level 1      141         26           26 |  193     Level 2       29          5           12 |   46      Level 3       43         10           10 |   63      -----------------------------------------------     Total        213         41           48 |  302      Thank you very much for your help.


title: Using geometric mean / geometric standard deviation in statistical process control chart

text: I am trying to control a binomial proportion that hovers pretty close around 0.1% - 1.0%. Since it seems like dropping from 0.1% to 0.05% is equally as severe as increasing to 0.2% (halving the occurrences, doubling the occurrences), it makes sense to me to use geometric mean / standard deviation. Also, by using geometric standard deviation, it seems I have an added benefit of not having the lower limit become less that 0.0% (an impossible value for my data). Are there any examples that support using geometric calculations in statistical process control or is it a bad idea altogether?


title: Can individual odds ratios be added to get one pooled odds ratio to compare to a reference group?

text: I have the admission for 4 different hospitals (coded as binary) which are regressed on a list of independent variables. The multinomial logistic provides the odds ratio (O.R.) for independent variable for each hospital and then compare this to a reference hospital. Can I add each of the O.R. from each independent variable and get a combined/pooled O.R. for each hospital and then compare this pooled O.R. to a reference hospital and to the individual hospitals?


title: Detrending & cross-correlation function

text: I am looking for some help with my time-series data. What is the best method of detrending/transformation of these two variables, so I do not violate assumptions of stationarity when applying a cross correlation function (to find out if one series is leading the other)? ![enter image description here](http://i.stack.imgur.com/FnWsF.png) The cross correlation function shows that the raw series are correlated when HDA is shifted back two year, but applying this without transformation of data and error structure would be wrong. Because I am working with scaled variables from 0 to 1 and because my x variable contains zeroes, I am not sure what is the best way to go about this (I can correct for non-constant mean via OLS detrending, but don't know how to address non-constant variance and non- independence of errors). Could this type of question be answered with wavelets? The y variable is a time series of biomass (B) scaled to the maximum value of B. The x variable is a time series of the proportion of high density areas (HDA), scaled to the maximum value of HDA > y<-c(0.84685420,0.64096448, 0.61250603, 0.49262176, 0.34298023, 0.39499759, > 0.23326153, 0.27661148, 0.52848738, 0.66898569, 0.71739592, 0.74696673, > 0.78385469, 0.92071371, 0.96981193, 1.00000000, 0.95564700, 0.83825109, > 0.67275358, 0.59576917, 0.42636232, 0.33447999, 0.30739109, 0.14768687, > 0.14132776, 0.08885388, 0.07370519, 0.08349140, 0.05824787, 0.04986337, > 0.04616621, 0.04828163, 0.04666131, 0.02836843, 0.03283073, 0.09343192, > 0.06804694, 0.06146279, 0.12578685, 0.17464716, 0.30159781, 0.33469217) > > year<-c(1970:2011) > > x<-c(0.00000000, 0.27272727, 0.09297521, 0.53827751, 0.17786561, 0.09977827, > 0.10765550, 0.49889135, 0.29933481, 0.77922078, 0.85623679, 1.00000000, > 0.47568710, 0.97402597, 0.39430449, 0.56426332, 0.31774051, 0.47661077, > 0.27579162, 0.46487603, 0.38654259, 0.31168831, 0.15151515, 0.15293118, > 0.00000000, 0.22113022, 0.00000000, 0.00000000, 0.16201620, 0.15437393, > 0.00000000, 0.14229249, 0.00000000, 0.00000000, 0.00000000, 0.00000000, > 0.14354067, 0.08704062, 0.00000000, 0.23830538, 0.07718696, 0.05928854) > > ccf<-ccf(x,y,type="correlation", main="x=HDA/max(HDA), y=B/max(B)", > ylab="ccf") original biomass is > y<-c(131707, 99686, 95260, 76615, 53342, 61432, 36278, 43020, 82193, 104044, > 111573, 116172,121909, 143194, 150830, 155525, 148627, 130369, 104630, > 92657, 66310, 52020, 47807, > 22969,21980,13819,11463,12985,9059,7755,7180,7509,7257,4412,5106,14531,10583, > 9559,19563, 27162,46906, 52053) > > logy<-log(y) and I am looking at the scaled biomass so then that would be > logy.scaled<-logy/max(logy) ![enter image description here](http://i.stack.imgur.com/KWCbm.png) what would be the best way to apply a ccf? or other method to see the red series leads the blue in time?


title: Appropriate probability threshold for Jarque-Bera test

text: If a probability distribution has, say, 112 bins with around 29000 samples, with the maximum probability of a bin being less than 0.05, is the Jarque-Bera test an effective measure of conformance to a Normal distribution? I'm working with univariate data, whose probability distribution 'looks' normal, with some deviation. However, the Jarque-Bera test for normality indicates that it is _not_ normal. Also, since n ~ 29000, the JB statistic is inflated to a ridiculous value like 1474523. Am I doing something wrong? Any help would be most welcome. The Jarque Bera statistics and a histogram for these data are:               n = 28535.0     Skewness = -4.24685767755943     Kurtosis = 37.1765663807246     JB statistic = 1474523.40413686      ![Histogram](http://i.stack.imgur.com/NboDn.png)


title: Log-likelihood distance measure validity for clustering

text: I have calculated log-likelihood distances between 50 sequences according to the Formula (1): $$ D(X_i,X_j)= 1/2(\log p(X_i|Mod_j)+\log p(X_j|Mod_i)), $$ where $ p(X_i|Mod_j) $ is the likelihood of sequence $X_i$ being produced by model $Mod_j$, where $Mod_j$ is a corresponding Markov model of the given $Seq_j$, defined by its Transition Probability Matrix and Start Probabilities Vector. The measure is symmetrical as seen from the definition. To make the measure more "legible" and similar to the traditional measures, I compute distance$=(1-D)$ from formula (1). Thus, $D(X_i,X_i) = 0$ and the distance increases if the likelihood decreases. Now, I have a 50x50 Distance Matrix.I have run a "meaningfullness" check, and it seemed ok for me - i.e. more similar sequences had smaller distance and very different ones had very large distance. The distances seemed to satisfy the triangle inequality. However, I have noticed that: 1) the shorter sequences seem to be "closer" to all other sequences than longer ones. It seems that this distance measure is biased to favor short distances. 2) I have tried PAM-clustering with the distance matrix by converting my distance matrix to dist object in `R` by using as.dist(), and my results were very bad, even for 2 clusters or 49 ( max avg.silhouette width produced by `R` function pam was 0.28). With some numbers of clusters the avg.silhouette widths were even negative. I am coming to conclusion that my way of computing medoids is invalid/conceptually wrong. What could be the problem? Can log-likelihood distance matrix be used with medoids clustering at all? edit: I am including the heatmap of the distance matrix, where x- and y-axis represent sequences (1 through 50th). It looks strange to me but I cannot pinpoint what exactly doesn't feel right. ![heatmap](http://i.stack.imgur.com/RcSBc.png)


title: Conflicting results with GoF tests

text: I am a computer scientist with little statistics background and I am trying to find the best fitting distribution for some data set (using MATLAB). To assess the goodness of fit I use both Kolmogorov Smirnov (KS) and Anderson Darling (AD) tests, and here are the p-values for the same data set:               Distribution          AD Test     KS Test     Exp                   0.439       1.49e-7     Weibull               0.498       1.40e-6     Pareto                0.244       6.24e-14     Logn                  0.684       2.69e-4     Gamma                 0.595       2.16e-4      I use a significance level of 0.05, and as far as I know with a p-value < 0.05 the null hypothesis is rejected which is the case for the KS test results. Then what should I conclude based on this result? The KS test says that none of the distributions is a good fit while the AD test can't reject the null hypothesis. Edit: Here is an overview of what we do in our code:               fit_functions = { @wblfit, @expfit, @lognfit, @gpfit, @gamfit};     for i=1:length(fit_functions)         [varargout{1:x}] = fit (param, fit_functions {i});         [ad_result ks_result] = run_gof_tests(param, cdf_functions{i}, ... varargout{:});      and in run_gof_tests function we average 1,000 p-values to calculate the final p-values. Each p-value is computed by drawing fifty samples randomly from the data set. We have used this method due to reasons described in this tech report on pg 12: Modeling Machine Availability in Enterprise and Wide-area Distributed Computing Environments by Nurmi et al., UCSB Computer Science Technical Report Number CS2003-28. Thanks!


title: How to apply parametric and non-parametric bootstrap using R code?

text: I think its better to put like this. I want to use bootstrap for a test statistic. how can I apply the parametric bootstrap., and the non parametric bootstrap on the other hand. and what is the big difference of the two methods. it will be appriciated if you answered me with some examples. thanks.


title: Running transformed data with a Generalized Linear Model in SPSS

text: I am running a generalized linear model in SPSS on data highly skewed to the right with a bunch of zeroes (average # caterpillars per tree species branch sampled) so I decided to use the Tweedie distribution. In the model I ran both untransformed data (with `link=log`) as well as $\log(x+1)$ transformed data (with `link=identity`). The latter model had a much smaller (more negative) AICc value than the untransformed data with `link=log`. Is it valid to run the GLM with $\log(x+1)$ transformed data if `link=identity`? Or am I violating some kind of assumption about the model?


title: What is the variance-covariance matrix of the OLS residual vector?

text: For example, take the generalized regression model: y = X$\beta$ + $\varepsilon$, what would the variance-covariance matrix be for $\varepsilon$ hat?


title: how would you approach giving an introductory class about Bayesian statistics?

text: I need to give a lecture about Bayesian statistics, introducing it to people who have already basic knowledge of classic statistics (but not too much of it in general). I want to start with some example, showing the main differences between Bayesian statistics and classic statistics. Is anyone familiar with a very basic intro that I can derive slides from, or actual slides I can use? I think the best would be to start with an example... I would rather not have a coin toss example, but I can be persuaded to do that :-) I would also be happy to talk about some of the history.


title: has anyone implemented an autoencoder with random forests

text: I'm interested in exploring autoencoders which can be used to develop a compressed representation of data useful for machine learning. In my experience random forests are easier to work with and more flexible than linear models and so I'd like to try to use them to build an autoencoder. One might do this by using random forests to predict multiple outcomes then re-represent each data point as a binary sequence corresponding to the braches that it took. For example if a forest consisted of two trees with three branches then the code 011 101 would represent a datapoint that took the second and third branches of the first tree and took the first and third branch of the second tree. Is anyone familiar with work like this? I am interested in papers, implementations of multi-outcome random forests, and techniques that convert random forests into binary representations of data points. Edit: clarifying


title: combine two Normal Dists, with one the mean of the other

text: For an assignment, I need to find $E_{y|\theta} [- \frac{\partial^2}{\partial^2\theta} \ln[P(y|\theta)] ]$, where $y \sim \mathcal{N}(\theta, 1)$ and $\theta \sim \mathcal{N}(0, \delta^2)$. I can handle taking the partial derivatives, log, etc. Where I am stumbling is in finding $f_{y|\theta}$. I keep coming back to simply plugging in $\theta$ for $\mu$ in the regular normal dist. equation. But I can't believe that this is correct since otherwise, where would the effect of $\theta$ 's distribution be represented? Perhaps another way I might ask this is: How is $P(y|\theta)$ different from $P(y)$? Any assistance would be greatly appreciated. Thank you in advance


title: How to estimate the mode using non-parametric methods of a 4-variate random vector drawn from a continuous multivariate distribution?

text: I have a sample of size 10,000 of a 4-variate random vector coming from a (unknown) continuous multivariate distribution. **How can I estimate the mode of this density using nonparametric methods?** I am currently estimating the corresponding density using the sample and a kernel estimator. Then I take the value in the grid that maximises this estimator (there are several theoretical results that justifies this procedure). The packages `ks` and `np` are _very_ slow in this case, particularly in the calculation of the bandwidth matrix. Do you have any suggestion?


title: Assumptions needed for multiple linear regression

text: Are assumptions for multiple linear regression basically the same as simple linear regression? One has to check for linearity for each of the continuous predictors versus the outcome variable? If there are categorical variables, one has to check the linearity between each of the continuous predictors with the outcome variable for all levels of the categorical variables?


title: Computing Mallows' Cp with SAS

text: I have a question about SAS: How do you get the parameter estimate for Mallows' Cp procedure in SAS?


title: Significance testing or cross validation?

text: Two common approaches for selecting correlated variables are significance tests and cross validation. What problem does each try to solve and when would I prefer one over the other?


title: Comparison of two linear regression models with squared regressors

text: I need to compare two linear regression models that include the regressors in levels and in squares: $Y=a_1x^2+b_1x+c_1$ **and** $Y=a_2x^2+b_2x+c_2 $ Specifically, I need to test whether these regression models are significantly different (i.e., whether $a_1=a_2,b_1=b_2,c_1=c_2$). In the forum, I read that one should build respective interactions. Yet, I do not exactly understand what is meant by including “interactions.” What would that model look like? What should be compared with what, and what test statistic is best? (I do not know Stata, unfortunately.) Please advise. Both models are based on /using the exact same population.


title: Alternative to multi-class classification?

text: Problem Setting: I have a fixed number of class labels $|Y|$. In the test phase, I will be given a $x$ and a candidate $y \in Y$. I have to make a "yes/no" assertion whether $y$ is the correct class label for $x$. An easy way would be to use multi-class classification and output "no" for every pairing which is incorrect according to my classifier. For example, _if classifier predicts $y_3$ for $x$, then $\forall y \neq y_3$ output "no"._ But this seems like an overkill for this problem. Is there an alternative solution? I am thinking of using a maximum entropy classifier to obtain probability scores for each possible pairing of $(x,y_i)$.


title: How to choose q threshold for pFDR?

text: We are doing an activation map study and we want to control for false positives in multiple comparisons. Our map is around 5000 pixels. There's a lot of noise and variance so we see results that seem appropriate at p < 0.05 . About 10% of the map is significant activation. We have run the Storey pFDR procedure to get q values, generating the pi-zero with the QVALUE GUI program in R. We don't get one pixel with q < 0.05 . A few start coming in at .14-.15, with a lot of values in the .20-.26 range, .26 probably being the mode of the distrib rounding to .00, then a reduction. Now, Storey is clear ( in his PNAS article among others ) that the q value is a choice that doesn't have to match the p value for any particular reason. The PNAS supplmentary materials show some genetics values and the q values are often much much higher. I can't see a great way to choose my q threshold. Any suggestions?


title: Can time series data have both unit root and structural break?

text: My data rejects unit root, but shows structural break, is this possible? Thank you!


title: Pearson correlation is not significant although effect is .30 What could be the reason?

text: I have a correlation between two variables of 0.30, hence a medium effect. but it is not statistically significant and I want to report the reasons for that. looking at the test statistic I see that the sample size (as in all test statistics) has an influence. But what else could be the reason? Does the standard deviation or standard error have an influence? edit: Thank you for your helpful and fast answers!


title: Nested vs. conditional logit regression

text: I am trying to estimate a logit regression model with travel mode choice (categorical) being the dependent variable; explanatory variables include age (categorical), income (categorical), gender (binary), and travel cost (continuous). I also have information on origin town (6 digit code for each town, about 50 towns) of the travellers, which I want to control for in the logit model. I am using R. I saw some other posts on StackOverflow saying that similar problems could be solved using a nested logit (origin town as the nest). While some of my research suggested a conditional logit could be an option. Could you guys advise on what the key differences are between nested and conditional logit in solving my problem in this case?


title: Reporting of an interaction in a binary logistic regression

text: I've found some interesting results that I'm trying to write up appropriately, but I'm having a hard time finding any guidance into how to write up an interaction in a binary logistic regression (outcome is 0,1). The interaction was predicted, and this is not an issue. The issue is that I have two categorical predictors. One has 2 levels Var1(1 and 2) the other Var2 has 3 levels (creatively 1,2,3 in the coding). We have a significant Var1*Var2 interaction. Our output lists the beta, SE, Wald,df, Sig, and Exp(B) for Var 1, but due to the 3 levels of Var2, we get only the Wald, df and Sig columns for Var2 (overall, I assume this is a type of omnibus-like test) and the full set of results for Var1 at levels 1 and 2. A similar situation happens with the interaction term. On the overall interaction we have wald, df and sig, and then Var1 (level1) x Var 2 (level1) and Var1(level1) x Var2(level2) interactions with all 5 results.   Our premise is that we expect the interaction, and then do the appropriate follow up contrast test. This we can handle (I think). Our hangup is basically just how to report the interaction and the model results. The only examples I can find are either only 2x2 interactions (which don't have this complication) or no interaction binary regression write ups. Is there any reference for how to write up these findings that someone could point me towards? I guess secondly too, is there any good source for explanation of the difference between model fit statistics that SAS spews out? Likelihood Ratio, Score and Wald all seem to test whether the model as a whole with variables, is better than the intercept only model. From what I can see the LR is the preferred method (and the only one that SPSS reports by default), but is there any reference that could tell me why SAS lists the other two, or better when they would be preferrable measures? THanks, Ned


title: How to rebalance a dataset for a logit model in Stata?

text: I want to use logit model in Stata. In my dataset only 3% of my target variable observations are 1, and 97% are 0. How do I rebalance this data set, so I have more observations labelled 1, and fewer labelled 0?


title: 2d interpolation method for coarsely sampled image

text: I'm looking for a general method for 2d interpolation of a coarsely sampled image. I'll use an example, taken from the scipy.interpolate (Python) page. ![](http://bit.ly/10lHXma) Say, I have this image, but instead of interpolating in the sense of going from e.g. a 512x512 grid to a 1024x1024 grid, I want to use the red/blue spots as data and re-define the unobserved green areas accordingly through an interpolation (so, the area between two red spots would be more yellow than green, and the area between two blue spots would be blue, not green). Of course, I'll do this programmatically. Language isn't too important, just looking for a method. I've tried kriging (available in IDL/Python), but I don't understand it enough to constrain the parameters, and for large (above ~64x64) arrays, it takes a lot of CPU/memory capacity. Anyway, let me know if the question is too vague, and thanks...


title: How can I remove the z-order bias of a coloured scatter plot?

text: I'm currently working on a plot engine for my project. This engine should be robust for a wide range of inputs. In order to analyse the data, I'm plotting a series of graphs utilising _python/matplotlib_. Among them is the following: ![scatter plot](http://i.stack.imgur.com/6vhgg.png) I think this graph is not good because the data being plotted first (high pressures, red) have a lower z-order (i.e. they are overdrawn) than the blue bullets for low pressures. Thus introducing a bias when looking at the graph. The underlying reason for that is that the data is _bell-shaped_. First off, do you agree or disagree? I could leave it like it is because it is just one of many _views_ on the data. It could still be useful. However, if there is a way to make this graph better with some sort of trick, I'd be much happier. I already played with point size, _transparency/alpha_ and _edgecolor_. This only made it worse. A great way to remove the z-order bias in scatter plots is to bin the data and colour-code it accordingly (e.g. hexbin). But since I used the colour for the pressure information, I see no possibility to to something similar. Another idea would be to randomise the z-order, but I'm not sure how to do that and if the result would be better. Any other comments for improvements are appreciated.


title: How to model bounded target variable?

text: I have 5 variables and I'm trying to predict my target variable which must be within the range 0 to 70. How do I use this piece of information to model my target better? Thanks


title: Defining quantiles over a weighted sample

text: I have a weighted sample, for which I wish to calculate quantiles.1 Ideally, where the weights are equal (whether = 1 or otherwise), the results would be consistent with those of `scipy.stats.scoreatpercentile()` and R's `quantile(...,type=7)`. One simple approach would be to "multiply out" the sample using the weights given. That effectively gives a locally "flat" ecdf in the areas of weight > 1, which intuitively seems like the wrong approach when the sample is actually a sub-sampling. In particular, it means that a sample with weights all equal to 1 has different quantiles than one with weights all equal to 2, or 3. (Note, however, that the paper referenced in [1] does appear to use this approach.) http://en.wikipedia.org/wiki/Percentile#Weighted_percentile gives an alternative formulation for weighted percentile. Its not clear in this formulation whether adjacent samples with identical values should first be combined and their weights summed, and in any case its results do not appear to be consistent with R's default type 7 `quantile()` in the unweighted/equally weighted case. The wikipedia page on quantiles doesn't mention the weighted case at all. Is there a weighted generalisation of R's "type 7" quantile function? [using Python, but just looking for an algorithm, really, so any language will do] M [1] Weights are integers; the weights are those of the buffers which are combined in the "collapse" and "output" operations as described in http://infolab.stanford.edu/~manku/papers/98sigmod-quantiles.pdf. Essentially the weighted sample is a sub-sampling of the full unweighted sample, with each element x(i) in the sub-sample representing weight(i) elements in the full sample.


title: Uniform frequency from non-uniform (exponential) distribution?

text: I lack knowledge of statistical terminology, so I'll try to thoroughly explain my predicament and hope that it is understandable: I am programming a software for technical analysis of financial markets.   This software will receive a variable who's value represents current market conditions. In this particular case the variable represents a measurement, in standard deviations, of the financial vehicle's price in respect to its moving average. The exact nature of this variable is not relevant to my question - the important part is that it has non-uniform distribution. What I mean by this is that at any given time the variable has a much greater chance of being somewhere around 0 than of being at one of the extremes. Now, the point of my software is to find and evaluate _similar conditions_ in the market history. If I look for the exact same value (ex. 2.0021) I will not find it; so what I must do is create an interval around 2.0021 and evaluate whatever values fall within the interval. Then, when I receive a different input value, I will need to create an interval of a **different size,** such that approximately **the same amount of values will be found within this interval.** ## Here is a concrete example: > `CURRENT VALUE X`a `= 2.0021` Lets say I found that `100 values of X`a fall within the interval of `1.9992 < X < 2.0050`...   (The size of the interval (obviously) is `0.0058`) Next I receive another value for X: > `CURRENT VALUE X`b `= 4.5010` **Values around 4.5 are much more rare than values around 2.08**. Therefore, the interval surrounding `4.5010` must also be larger in order for my software to find the same number of results as (close to) Xa ; I am trying to create intervals in my data such that each interval contains the same number of samples within its bounds. In other words: The `CLOSER TO 0`, the more data points are present, the `SMALLER` the interval needs to be. The farther from 0, the fewer data points are present, the greater the interval needs to be in order to contain the same amount of results. And now to the heart of my question: How can I formulaically calculate the size of the interval around `X` based on the value of `X`? It is also crucial to mention that I have a ton (500,000) of real `X` values to work from. This is what I tried - I feel like I'm close but not quite there: I created a histogram out of my 1D data. Since the data is symmetric, and centered on 0, I created another histogram based on the data's absolute value. My Histograms I do not know where to go from here :( I tried doing an exponential regression to the data; it fits very well but I do not know how to use the formula to create the sort of bounds within my data that I want. I seriously regret not taking a Stats class in college... Can anyone help me out please?


title: Convergence in distribution from ratio of the two densities

text: I want to prove $[X|y] \rightarrow \delta_0(X)$, as $y$ goes to $\infty$, i.e. the distribution of $[X|y]$ converges to degenerate distribution at zero for large enough $y$. I know the density $f(x|y)$ upto a multiplicative constant and also the fact that $$ \frac{f(x|y)}{f(0|y)} \rightarrow 0 \quad \text{as } y \rightarrow \infty \quad \forall x \in \mathcal{D}(x)\setminus \\{0\\} $$ Using this fact, how do I prove $[X|y] \rightarrow \delta_0(X)$? (Or provide a counter-example).


title: How much variation is necessary for regression analysis

text: Let's say that I'm trying to find the factors which help predict whether I wear red shoes or not during a year long time-period. I could look at weather, day of the week, month, and other factors. However, how much variation does there need to be in the dependent variable for me to be certain that I am acctually measuring something that is meaningful. Do I need to have worn red shoes 2 times? 5? 50? 100? etc? Basically, how much variation in the dependent variable is necessary to develop a statistical model that generates somewhat meaningful results? Let's say that I have 150,000 rows and I am trying to explain variation that comprises 3.2% of that data. Is that enough?


title: Analyzing treatment effect with possibly flawed control data

text: I've got some rather messy data from a natural experiment. A number of subjects were measured (the measurements were hopefully-Poisson distributed counts and associated offsets), placed on a treatment at time T, and then measured again sometime afterward. However, there is some reason to suspect that the dependent variable has been changing over time, so we found a control dataset and measured those subjects before and after T as well. Unfortunately, I don't know which of the subjects in the control dataset were on the treatment - all I know is the ones who were on stayed on and the ones who were off stayed off. I can obtain this information, but it is rather time- consuming/expensive. There are two questions I would like to answer:   1. Does the treatment have an overall effect?    2. Does the treatment affect the majority of subjects? I don't really know how to answer 2 (some variant of McNemar's test?) so I would appreciate some advice there, but for 1 I've been setting up the problem something like this:               glm(counts ~ as.factor(subject.id) + before + offset(log(observation.time)),family=quasipoisson)      where before is coded 0 or 1. So I have two rows per subject. I've done that regression for both the control and test datasets, and the respective confidence intervals of before don't intersect, so I've been mildly optimistic. What is the best way to combine the two datasets in a single analysis, though? If I knew which of the control dataset had the treatment and which didn't, it seems fairly easy, but like I said, I don't know that.


title: Cochran-Armitage test for trend in a contingency table

text: I have a 2x4 contingency table and I'm testing for a monotone trend. How can the weights (scores) be applied in the following situation: The 1st category is expected to have no effect on the dichotomic variable. The 2nd and 3rd to have some effect. And the third to have by far the most. Until now, I applied weights 0, 1, 2, 3. The trend was significant. However, the scores were not optimally set for the test to be locally most powerful. Is it possible to set the weights such as: 0, 1, 1, 3 for example? I'm new to this test, so any advice as how to treat it would be helpful.


title: Formula for lmer()

text: **How does the formula for the`lmer` function work?** Some examples:               library('lmer')     lmer(Yield~(1|Batch),Dyestuff)     lmer(Yield~(Batch)*(1|rep(c(1,2,3),10)),Dyestuff)     lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy)      On the left of `~` sign is the dependent variable and on the right comes all the independent variables separated by the `+` sign. Is there any difference between `+` and `*`? What does `1|Batch` mean? Does it have to do with the intercept? With being a random or a fixed effect? What does `Day|Subjects` mean?


title: Correlating two questionnaires with grouped items

text: I need to correlate employee engagement (gathered data using the 9 item UWES questionnaire) and organizational commitment (gathered data using the 18 item Organizational Commitment Scale). The both of them can be divided into different components;   * UWES - `Vigour` (3-items in questionnaire), `Dedication` (3 items in questionnaire) & `Commitment` (3 items in questionnaire).   * OCS - `Affective Commitment` (6 items in questionnaire), `Continuance Commitment` (6 items in questionnaire) & `Normative Commitment` (6 items in questionnaire). I have an R-matrix set up correlating all the individual questions from both questionnaires, using Pearson's _r_. However I would like to correlate the components, and not individual items from the questionnaires. ![enter image description here](http://i.stack.imgur.com/ZjT4U.png) What would you recommend I use? Thanks!


title: Intuition behind (statistical) completeness

text: I was wondering if any of the members of this community would like to share his/her intuition about **completeness** in statistics. For the sake of "completeness", here's the definition, taken from Wikipedia: > Consider a random variable $X$ whose probability distribution belongs to a > parametric family of probability distributions $P_\theta$ parametrized by > $\theta$. > > The statistic $T$ is said to be complete for the distribution of $X$ if for > every measurable function $g$ (which must be independent of $\theta$) the > following implication holds: $$E(g(T(X))) = 0 \mbox{ for all }\theta \mbox{ > implies that }P_\theta(g(T(X)) = 0) = 1\mbox{ for all }\theta.$$ Now, to give you an idea of what kind of answer I am looking for, here's how I think about (minimal) sufficiency: I think of a statistic as a partition of the sample space. In that context, a statistic is **sufficient** for $\theta$ if this partition does not result in a loss of "information" about $\theta$; it is **minimal sufficient** if it is the coarsest partition which does not result in a loss of information (superlatives carry a uniqueness connotation, which I am ignoring here). **Note** : I am aware that a very similar question has been asked before, but I'm looking for a different answer.


title: Which machine learning method to use for geographic systems prediction?

text: I am trying to do experiments on geographic systems prediction. We're working on classifying the location where we sell product most. So, we need to analyze the hestorical data and predict the success of our product in specific areas, specific time. I don't have a lot of experience in statistics and modeling data beyond a high school statistics course so I'm kinda confused.               lat             long          city_block ..   event      date     41.98302392, -87.71849159,      23             1        17/04/2014       41.77351707, -87.59144826,      44             0       17/04/2014       41.77508317, -87.58899995,      24             0        17/04/2014       41.77511247, -87.58646695,      33             1        17/04/2014       41.77514645, -87.58515301,      21             0        17/04/2014       41.77538531, -87.58611272,      22             3        17/04/2014       41.71339537, -87.56963306,      39             0        16/04/2014      41.81685612, -87.59757281,      12             0        16/04/2014       41.81697313, -87.59910809,      02             3        16/04/2014      41.81695808, -87.60049861,      3              0        16/04/2014      41.75894604, -87.55560586,       4             0        16/04/2014       Event: represent number of items sold in location (lat/long) on (date). Can anyone give an example how to classify geographic systems, and predict values based on data we have. as far as I know, I can't use linear modeling (as I have two variables lat/long). Someone recommended to use density estimation Instead of using clustering. Using one-class SVM, is given points only from one class, and expected to learn a separation between members of that class and anything else. I read in another post that I have to use poisson, as it's counting model. Meanwhile, as I have a lot of negative points (90%), I will have problems finding good results. What do you suggest as an approach to tackle this problem ?


title: How to do multivariate machine learning? (predicting multiple dependent variables)

text: I am looking to predict groups of items that someone will purchase... i.e., I have multiple, colinear dependent variables. Rather than building 7 or so independent models to predict the probability of someone buying each of the 7 items, and then combining the results, what methods should I look into to have one model that accounts for the relationships between the 7 related, dependent variables (things they can purchase). I am using R as a programming language, so any R specific advice is appreciated.


title: How do we analyse likelihood in a dataset?

text: I am working to analyze poverty rate using census data. I have a huge dataset. I want to extract the likelihood from this dataset in order to create patterns for energy consumption. What is the best approach to tackle this problem?


title: Why is the decision function for probabilistic models a quotient (when we only consider two models)?

text: Take for example, that we want to find the probabilistic model for only **two** document types (doc can be + or -). I was trying to understand why the way that we classify a document model was with the following decision rule (under the assumption equal probability of the two type of models): $$log{\frac{P(D|\theta^+)}{P(D|\theta^-)}} = \begin{cases} \geq0, & \text{then +} \\\ <0, & \text{then -} \end{cases} $$ Why didn't we for example, choose to be the deciding expression: $$\frac{logP(D|\theta^+)}{logP(D|\theta^-)}$$ I feel it has something to do that what we originally had in mind to use: $$\frac{P(D|\theta^+)}{P(D|\theta^-)}$$ but probably using logs made is more mathematically convenient or something. But my point of my question is, why do we use the first decision rule that I posted? What is the intuition? Is the reason we are using it is for mathematical convenience or some deeper reason? What is wrong with the other two suggestion I have posted? * * * Notation: D = some document (just a bag of words for the document). +,- = two different document classes. $\theta^+$ = the parameter/model we have for the document class +. $\theta^-$ = the parameter/model we have for the document class -.


title: decision boundary of support vector machine when data is not linearly separable

text: Screenshot from this video: ![enter image description here](http://i.stack.imgur.com/xFSeO.png) This describes the decision boundary of support vector machine as a optimization problem with two constraints. But it seems to assume that the data points are linearly separable, what if they are not?


title: Can simple vector distance work as a SVM kernel?

text: I'm considering using the following simple function as an SVM kernel. It basically computes the distance between the 2 input vectors (norm): $K(\vec{x}_1, \vec{x}_2) = \left\| \vec{x}_1- \vec{x}_2 \right\|$, where $\vec{x}_1$ and $\vec{x}_2$ are $N$-dimensional vectors. I'm not familiar with documentation of such a kernel. Is it valid? Anyone experienced with such a kernel?


title: How to correctly interpret a parallel analysis in exploratory factor analysis?

text: Some scientific papers report results of parallel analysis of principal axis factor analysis in a way inconsistent with my understanding of the methodology. What am I missing? Am I wrong or are they. Example:   * **Data:** The performance of 200 individual humans has been observed on 10 tasks. For each individual and each task, one has a performance score. The question now is to determine how many factors are the cause for the performance on the 10 tasks.   * **Method:** parallel analysis to determine the number of factors to retain in a principal axis factor analysis.   * **Example for reported result:** “parallel analysis suggests that only factors with eigenvalue of 2.21 or more should be retained” That is nonsense, isn’t it? From the original paper by Horn (1965) and tutorials like Hayton et al. (2004) I understand that parallel analysis is an adaptation of the Kaiser criterion (eigenvalue > 1) based on random data. However, the adaptation is not to replace the cut-off 1 by another fixed number but an individual cut-off value for each factor (and dependent on the size of the data set, i.e. 200 times 10 scores). Looking at the examples by Horn (1965) and Hayton et al. (2004) and the output of R functions _fa.parallel_ in the _psych_ package and _parallel_ in the _nFactors_ package, I see that parallel analysis produces a downward sloping curve in the Scree plot to compare to the eigenvalues of the real data. More like “Retain the first factor if its eigenvalue is > 2.21; additionally retain the second if its eigenvalue is > 1.65; …”. **Is there any sensible setting, any school of thought, or any methodology that would render “parallel analysis suggests that only factors with eigenvalue of 2.21 or more should be retained” correct?** References: Hayton, J.C., Allen, D.G., Scarpello, V. (2004). Factor retention decisions in exploratory factor analysis: a tutorial on parallel analysis. Organizational Research Methods, 7(2):191-205. Horn, J.L. (1965). A rationale and test for the number of factors in factor analysis. Psychometrika, 30(2):179-185.


title: Time Series detrending with multiple polynomials

text: I'm currently working on energy demand forecasting using daily load data. I'm using data since 1990 and, in order to use ARIMAX models, I detrended the data (using a first-order polynomial) and then I removed seasonality. As we can see on this image ![monthly load data](http://i.stack.imgur.com/iQ5OG.jpg), in the last two years there is (unfortunately) a different trend with respect the previous 15 years. I was thinking to use two different first-order polynomials to perform detrend: one fitted on the data 1990-2008 and another one for 2008-2010. What is the best way to do this? I'm not sure if I must use a spline or not. Can you suggest me the best-practice method to perform this kind of detrending? Thank you, P.S. I'm a computer scientist and I'm using R (and MATLAB). **UPDATE** This is the loess smothing with span 0.9 ![LOESS span 0.9](http://i.stack.imgur.com/gRMMu.jpg) With span 0.5 ![LOESS span 0.5](http://i.stack.imgur.com/twiXz.jpg) and finally with span 0.1![LOESS span 0.1](http://i.stack.imgur.com/cgOWV.jpg) I will try to perform detrending with span 0.5. ** UPDATE 2 ** This is the comparison with the time-series with the original detrend (top, a single first-order polynomial) and loess with span 0.5 (bottom). There is an important difference in the last years. ![comparison between simple detrend and loess](http://i.stack.imgur.com/36Ryz.png)