title: How to detect a quasi separation problem for a data set?

text: Suppose that we have a two-column data set. One column consists of a hundred x=0 and a hundred x=1, whereas the other one consists of y's (1 or 0 response). Besides, suppose that the P(Y=1|X=0) = 0.001 and P(Y=1|X=1) = 0.05. When fitting a logistic model (glm(y ~ x)) to this data set, why are we having a quasi separation problem here? My doubts are as follows: Since P(Y=1|X=0) = 0.001, we have P(Y=0|X=0) = 0.999. On the other hand, P(Y=1|X=1) = 0.05, so P(Y=0|X=1) = 0.95. As a result, we would see Y=0 with a very high chance when x=0 or 1. Why does this situation imply a quasi separation problem? Next, why log[(0.05/0.95)/(0.001/0.999)] = 4 is the theoretical coefficient of x in the model? Thanks.


title: Comparing sales-person conversion rates

text: If I have two sales-persons for which I know their conversion rate (converted leads/total leads), and I want to find out whether the difference between the sales-persons is significant, what type of test could I use? Similarly if I have the conversion rate figure from last month for one sales- person and the conversion rate for this month, how can I discover whether this difference is significant?


title: Burns survival dataset - meaning of midpoint of set interval

text: I'm using the well-known USC Burns Survival dataset to explore logistic regression in R. The independent variable is burn area, and the outcome is binary survival (yes/no). In the documentation, burn areas are grouped as the 'midpoint of set intervals', and taken as a log ie.               log(area + 1)      for sample datapoints like (highest/lowest):               1.35    yes     2.35    no      Medically, burns are usually specifed in terms of surface area % of total body eg. 40%, or alternatively, as an estimate of surface area in square meters. My question is: how does this relate to the dataset independent variable ie. what does '1.35' actually mean in terms of % body surface area burnt? What would eg. 30% burns become in the dataset, using the 'midpoint of set interval'? Thanks guys Ref: http://statmaster.sdu.dk/courses/st111/data/index.html#burns


title: how to generate samples for model selection testing on adaboost

text: I want to assess some statistical methods for model selection on binary classification using adaBoost. For this, I have to generate artificial samples (input data) and create an oracle that has the form of an adaBoost function, i.e. a weighted sum of weak learners (in this case decision stumps). Now I'm wondering how I should generate the samples, and how I should choose the parameters for every stump, namely, feature, threshold and '>' or '<'. Will picking every option from a uniform distribution do? Any ideas?


title: What is a "contrast" in a Principal Component Analysis (PCA)?

text: I have been studying how to interpret principal components. I recently came across an example of a particular eigenvector: $$e_j^T = \left[ \frac{\sqrt{2} }{2}, \frac{-\sqrt{2} }{2}, 0, \dots,0 \right]$$ I am under the impression that, often in social sciences, "this is called a contrast." From what I can see, we have the difference between $X_1$ and $X_2$ (variables in our original dataset) playing a role in the $j$th eigenvector. What can this say about the relationship between $X_1$ and $X_2$, with respect to this eigenvector? There is clearly something going on with the equal, yet opposite, contribution of these independent variables, to our (possibly first) principal component. I guess the answer to this is the answer to the question, "what is a contrast?" Strangely enough, googling hasn't coming up with much of anything.


title: Sampling from a marginal when full density is given

text: I want to sample from a distribution with density $$ f(\mathbf x) = \int f(c) \prod_{i=1}^n f(x_i|c) dc $$ where $\mathbf x=(x_1,x_2,...,x_n)$. In my particular setup, is easy to sample according to the densities $f(c)$ and $f(x_i|c)$, but it is not obvious how to sample from the marginal given above. Would the following work?   1. Sample $c$ according to $f(c)$.   2. Sample $x_i$ according to $f(x_i|c)$, for $i=1,2,...,n$. Heuristic argument for why this should work:   * Obviously $f(\mathbf x,c) = f(\mathbf x|c)f(c)$, so the $(\mathbf x,c)$ pairs generated by my procedure above are distributed according to $f(\mathbf x,c)$. Marginalization amounts to ignoring the $c$'s in any pair $(\mathbf x,c)$. So the $\mathbf x$'s produced in step 2 of my procedure, taken by themselves, are distributed according to $f(\mathbf x)$, as desired. Is this sensible? Is there a short way to make it rigorous if so? Thanks!


title: Finding correlations from multiple irregular time series

text: I have data from different sources that is currently in multiple data frames. They generally look something like:               TimeStamp (MS Detail Converted with strptime), Field1, Field2, Field3       In all cases the Time Stamps are irregular, and the fields are either character data or int data. There are three types of correlations I would like to be able to do with this data (using the term loosely). In all cases these will sometimes be from the same source, and other times from different sources:   1. Numeric to Numeric   2. Numeric to Categorical (string). I am not sure of the proper term for this, but if field1 where to have values like A, B, and C, how might find if field2 has higher numbers than B is more likely to be the value of field1 The end goal is not only to test for correlations, but also to be able to discover them. I realize this is a pretty vague / big area. Mostly looking for basic examples, or point in the right direction of how to transform my data into format that will let me do these things.


title: Another 3 level hierarchical regression

text: I have a dependent variable: rtln (log transformed reaction time) and 3 predictors: choicenum (ranging from 1 to 6), ifrelevant (0 and 1), and condition (0 and 1), and I have a variable with subject IDs. I would like to estimate a model where choicenum is nested within subjects. I use the following notation:               lmer(rtln~choicenum + ifrelevant + condition       + (1|subject/choicenum),    data=myData)      for some reason this estimates 5 fixed effects for choicenum.   * Here's the output: http://dl.dropbox.com/u/22681355/output.tiff   * Here's the dataset: http://dl.dropbox.com/u/22681355/data.csv   * Here's the stata output of what I'd like to get: http://dl.dropbox.com/u/22681355/output.pdf


title: What is the distribution of the Binomial distribution parameter $p$ given a sample k and n?

text: If I have a $k$ successes in $n$ bernoulli trials, does the parameter $p$ of the binomial distribution follow a well-known distribution? There are some methods to calculate confidence intervals for $p$, I'm interested in the distribution for the exact method.


title: Multilayer perceptron 3 classes

text: I need to create a MLP for 3 classes. I don't know to to do the structure. I thought that having one hidden layer (don't know how many nodes) and and output layer with 2 nodes would be enough, like using dummies variables. The activation function could be sigmoidal, so I would have 00, 01, 10, and 00. But I would need to put for example that 01 and 10 belongs to the same class, or I would have 4 classes. Do you know what it's the best way to structure this? And what are the consequences of doing it the way I did it? Thanks!


title: Surprisingly similar hypothesis tests

text: I would like to ask you for pairs or even triples of statistical tests which are, surprisingly, similar. By _similar_ I mean they test essentially the same hypotheses and lead to approximately the same p value and by _surprisingly_ I mean they are hardly found in statistical texts. One for the start: **Jonckheere-Terpstra-Test and test for Kendall's rank correlation** The nonparametric 1-way ANOVA with ordered alternatives (known as "Jonckheere- Terpstra-Test") is similar to the standard test about Kendall's rank correlation, see e.g. this link. The latter of the two methods has the big advantage of belonging to a meaningful measure of association that can also be equipped by a confidence interval.


title: How to specify specific contrasts for repeated measures ANOVA using car?

text: I am trying to run a repeated measures Anova in R followed by some specific contrasts on that dataset. I think the correct approach would be to use `Anova()` from the car package. Lets illustrate my question with the example taken from `?Anova` using the `OBrienKaiser` data (Note: I ommited the gender factor from the example):   We have a design with one between subjects factor, treatment (3 levels: control, A, B), and 2 repeated-measures (within subjects) factors, phase (3 levels: pretest, posttest, followup) and hour (5 levels: 1 to 5). The standard ANOVA table is given by (in difference to example(Anova) I switched to Type 3 Sums of Squares, that is what my field wants):               require(car)     phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),     levels=c("pretest", "posttest", "followup"))     hour <- ordered(rep(1:5, 3))     idata <- data.frame(phase, hour)     mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)     av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)     summary(av.ok, multivariate=FALSE)      Now, imagine that the highest order interaction would have been significant (which is not the case) and we would like to explore it further with the following contrasts:   **Is there a difference between hours 1 &2 versus hours 3 (contrast 1) and between hours 1&2 versus hours 4&5 (contrast 2) in the treatment conditions (A&B together)?**   In other words, how do I specify these contrasts:   1. `((treatment %in% c("A", "B")) & (hour %in% 1:2))` versus `((treatment %in% c("A", "B")) & (hour %in% 3))`   2. `((treatment %in% c("A", "B")) & (hour %in% 1:2))` versus `((treatment %in% c("A", "B")) & (hour %in% 4:5))` My idea would be to run another ANOVA ommitting the non-needed treatment condition (control):               mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")     av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)     summary(av2, multivariate=FALSE)      However, I still have no idea how to set up the appropriate within-subject contrast matrix comparing hours 1&2 with 3 and 1&2 with 4&5\. And I am not sure if omitting the non-needed treatment group is indeed a good idea as it changes the overall error term. Before going for `Anova()` I was also thinking going for `lme`. However, there are small differences in F and p values between textbook ANOVA and what is returned from `anove(lme)` due to possible negative variances in standard ANOVA (which are not allowed in `lme`). Relatedly, somebody pointed me to `gls` which allows for fitting repeated measures ANOVA, however, it has no contrast argument. To clarify: I want an F or t test (using type III sums of squares) that answers whether or not the desired contrasts are significant or not. * * * **Update:** I already asked a very similar question on R-help, there was no answer. A similar questions was posed on R-help some time ago. However, the answers did also not solve the problem.


title: Equivalent of PCA for regression with lots of examples?

text: I'm writing code to approximate a function $y=f(\vec{x})$ where $y\in\mathbb{R}$ and $\vec{x}\in\mathbb{R}^N$ for medium-sized $N$ ($N$ between 20 and 50, approx.). I have a **ton** of examples, however, so my number of examples greatly exceeds $N$. But, the examples probably are somewhat redundant, and some of the $N$ variables probably aren't correlated much with $y$. While I can compute $\vec{x}$ in near zero time, I need to evaluate $f(\vec{x})$ often/quickly. So, I'd like to project $\vec{x}$ onto a lower-dimensional subspace before learning/evaluating $f(\vec{x})$. I know the standard technique for finding a lower-dimensional projection given a lot of examples is PCA. But in the regression context, I care more about the correlation of the projection of $\vec{x}$ onto each of the basis vectors with the value of $y$ then simply maximizing the variance of the projections. So my basic question is: **how can I compute such a(n) (orthogonal) basis?** I've looked into "partial least squares" and "correlated components analysis," but if I understand properly, they will only compute a single basis vector since $y\in\mathbb{R}^1$. I want fewer than $N$ vectors in my (approximate) basis but more than one. If it matters, I'll be using a kernelized regression technique, like $\nu$-SVR, LS-SVR, or kernel ridge regression. [Apologies if this is a repeat of another post -- I come from a different community and don't "speak stats" very well! Thanks in advance for your help!]


title: How to get the quantiles by simulations?

text: I would like to find the quantiles of a random variable $X$ by simulations. I plan to do it in the following way, however, I do not know which one should be more correct though the results look very very close. Let's say p=0.99 A)1.Simulate 100 random numbers of $X$, order them, and pick up the 99th one. 2.Do 1 for 1e4 times and take the average. B)Simulate 1e6(100*1e4) random numbers of $X$, and pick up the .99*1e6-th one from the ordered sequence. I guess A may look may "correct", and the standard error can be computed. But B is much faster. However, I do not know how to compute the standard error with B. The reuslts are close indeed. Thanks a lot for you helps


title: Can I pass an "at" parameter for the x-axis locations of bars, to an R barplot?

text: Hi in base graphics in R, I would like to overlay a bar plot on a chart which has irregularly spaced points on a curve. The bars should represent the change in these points since the previous sample. But I would like the bars to align to the points in the horizontal direction, so there should basically be a lot of whitespace, and wherever there's a point in the line chart, a bar should also appear at the same x-axis location, showing the change. The bottom of bars will start either on the x-axis (i.e. at y=0), or possibly underneath as a separate chart. How would I go about doing this in R? I know I can do this with bloxplots, as you can see in the chart in question below, where the boxes align to the various points (which are, fyi, south african government bonds). Now how do I add the bar plot? ![enter image description here](http://i.stack.imgur.com/rov8u.png)


title: Seasonal data forecasting issues

text: I have a seasonally decomposed data set. The data set has strong seasonality. Now I am trying to fit the 'seasonal part' of dataset into ARIMA model and tried to forecast (with SPSS). The problem is, I get exactly same values in forecasts as that of actual values. So, MAPE is coming to be 0.000 Is this obvious to happen or am I doing something wrong? P.S. The data set is here: http://mihirsathe.com/mihir/STI/STI/drugs/index.html


title: Sales Forecasting using Support Vector Machine

text: I have sales data for last three years 2011-2013. I want to use Support Vector Machine technique in R to do the predictions. I just wanted to know that the approach that I am using is correct or not? I am taking the sales data of 2011 and 2012 as an independent variables and sales data of 2013 as a dependent variable as the input demands independent variables and a dependent variable. My aim is to predict the sales for year 2014, where I will take sales of 2012 and 2013 data as in independent variable.FYI it a time series data. Am I doing right or what else I can do? I would appreciate any kind of help.


title: ratio of explanatory variables in linear regression

text: I wonder if anyone has any links or advice on specifying a ratio of two explanatory variables in a linear regression ie specifying the two independent variables plus their ratio. We have data where the ratio term seems to be significant. Is there anything on functional form or any other issues eg collinearity ? thanks, Glen I meant we are considering a model y = a *x1 + b * x2 + c * x1/x2 where y,x1,x2 are all continuous and wondered if this was reasonable - or if there were some good refs on assessing appropriate functional form of multiple covariates many thanks


title: Methodology of modelling sparse events

text: Let us say that we have a process that generates sequences of the following form: $$S=[0,0,0,0,0,0,0,0,0,0,0,0,0,\mathbf{4},0,0,0,0,0,0,0,0,0,0,0,0,0,0,\mathbf{1},0,0,0,0,0,0,0,0,\mathbf{3},0,0,0,\mathbf{5},0,0,0]$$ What is the best way to model something like this? Since $\mathrm{E}(x)$ is small and the number of samples is low one would be tempted to model this using a Poisson distribution. However, this seems quite limiting as give that a Poisson is fully described with the same parameter $\mathrm{E}(x)$, so the same Poisson could model generate another sequence with more zero events and higher values, which I consider been drawn for a different process. Would it make more sense to first model the non-zero occurrences with a Binomial, draw from it and them model the number corresponding to the non zero occurrences with a Poisson?


title: Pattern recognition techniques in spatial or spatio-temporal data?

text: I am working with weather forecasters and have access to historical climatology data. Given current weather conditions in an area of interest (i.e. the current "map"), we want to try to find the most similar "map" from the past data. The idea is to try and make a weather forecast by finding the best analog from past data. The data is represented as a regular X by Y grid (i.e. a matrix) of points, where X is the horizontal location and Y is the vertical location, and the (X,Y)th value in the matrix represents the response variable Z at that location. In addition, the grid points are evenly spaced out. As an example, Z can be a measure of surface temperature, which is measured at each of the grid points. We take care of seasonal effects by restricting the search in the past data to a window of +/- 15 days of the test date. For example, if we want to find the best analog for a map from 2013-06-19, we would only consider maps from 2012-06-19 +/- 15 days, 2011-06-19 +/- 15 days, etc. We also restrict the search to observations taken at the same time as the test date. For example, if the test data is an observation taken at noon, then we will only look at the past data taken from the same time. I have two questions. (1) Given two grids (or "maps" or matrices) of data, how can I best calculate the similarity between them? Are there methods that take into account the spatial nature of the data? For example, point (1,1) will be highly correlated with the nearby point (1,2), etc. I am currently using a very simple distance metric, where I just take the difference of the two maps and find the Frobenius norm. The map from the past that yields the smallest value is the 'closest' map to the test conditions. (2) I am new to spatial statistics and I am looking for literature that relates to what I am trying to do. What should I read to become familiar with working with grid data? What resources are there to learn about pattern recognition in spatial or spatio-temporal data? (I want to mention that I am working in R, so I would welcome package recommendations as well!)


title: Plot graph with more than one value on x and y-axis in R

text: I have a dataset containing two columns and a total of 90 rows. The data is from my experiment where in the first column I have an integer representing the quantity, in the second column I have a percentage. A small example:               Quantity   Percentage     1          53%     1          51%     1          67%     2          73%     2          69%     3          73%     ...        ...      As you can see in both columns the numbers can occur more than once. Now I wish to plot this in a graph (I was thinking a scatter plot) in R. I just am a real beginner in using R and statistics so I was hoping someone can help me out how to get a good graph. If someone has an other suggestion that would give a better representation, shoot! I just need to have a visual representation that shows the correlation between the two values.


title: Actual vs. Target

text: I’ve a data set contains two columns: actual value and target value. Is there a way where I can express the difference between these two columns as a single percentage value (that’s standardized, i.e., between 0-100%)? So I was wondering if I can calculate the standard deviation for the difference between values and then how can I standardize it to be a percentage between 0-100%? Any help will be appreciated.


title: Rates, underlying distributions, and measuring things

text: When comparing two metrics, say visitors and time spent on a website, by just taking the ratio (time on site)/(visitors) can this be misleading? And if so how can we describe this problem mathematically or with statistics. More context, the vanilla interpretation of such a ratio is that the "typical" visitor spends **t** minutes on a site. What if we look at the actual Time on Site by Visitor, and this distribution is hypothetically multimodal, or heavily skewed, basically something that would not lend itself to being interpreted as having a typical visitor? What I'm curious about is functionally describing what's going on here. I feel the answer is in the realm of these are actually random processes that we're taking a ratio of and hence the dynamics of each process should dictate the dynamics of the ratio of the processes...right? In physical problems this is probably not often an issue as the randomness is probably mostly measuring errors but in other systems that are inherently noisy maybe this should be a consideration.


title: Using non-stationary time series in basic regression models

text: I know that non-stationary time series data cannot be simply used in OLS regression. But if I take the first difference and the log of the variables, making them more stationary around zero, can I use these transformed variables in the OLS regression to generate some useful results? Does the estimated coefficients and R-squared become more accurate?


title: Signs on logistic regression betas flip relative to observed - expected counts from chi-squared test

text: I conduct a chi-squared analysis on some bins and conclude that an association between the bins and an event exists. I then calculate logistic regression coefficients to validate my hypothesis. Also, I always look at the observed vs expected count in each bin (from chi-squared) to guesstimate logistic regression results. For example, if `Bin 1` had an observed count = 152 but an expected count of 85, I would estimate that the over all chi-square result is likely to be significant and Bin 1 is likely to have a positive coefficient from the logistic regression. However, this `Bin 1` has a negative coefficient from logistic regression. Is my understanding wrong - `when observed >> expected, logistic regression coefficient should be positive?`(assuming here that results are significant).


title: Ordered logistic multilevel regression in imputed dataset

text: Is there any way to calculate an ordered logistic multilevel regression in an imputed dataset? I have tried Stata’s `GLLAMM` (Generalized Linear and Latent Mixed Models), but that won’t work with the `mi estimate, post:`command prefix. Should I “estimate anyway”, or is there a better way to do it, preferably R or Stata?               > mi estimate, post : gllamm q18 [...] , link(ologit) f(binom) i(country) eform          mi estimate: command not supported     gllamm is not officially supported by mi estimate; see mi estimation for a list of     Stata estimation commands that are supported by mi estimate.  You can use option     cmdok to allow estimation anyway.      Also, suppose my multilevel model has only two levels, i.e. workers nested in countries. IIRC, in the case of linear OLS regression, such multilevel model would be identical to a ordinary model that includes country dummies. If I use that workaround with ordered logistic regression, will my results be biased?


title: Correlation of two interval scaled variables

text: I have two variables A and B of interval scale type. A is the birth year of a person and B is the number of years spent at university. What method would you recommend to measure the correlation?


title: Correlation of a nominally scaled variable and another of interval scale type

text: I have two variables A and B. A is of nominal scale type and B is of interval scale type. A is the brand of the car a person drives. And B is the person's mean salary per year. What method would you recommend to measure the correlation?


title: Is is possible (or advisable) to do Change Point Analysis on sequence of groups with R?

text: I'm familiar with post-hoc testing with `ANOVA` for exploring differences between a sequence of groups, but recently I've been reading about Change Point Analysis (especially the `R` packages `bcp`, `changepoint` and `strucchange`). It looks like those packages only handle data where there is one data point per unit of time. I'm curious if they can be used with data where there are multiple data points per unit of time. Here's some example data representing the measurement of a single continuous variable on a number of specimens that have been dated to specific moments in time (no repeated measurements):               a<-data.frame(time="1000",x=rnorm(10,12,3))     b<-data.frame(time="2000",x=rnorm(50,13,4))     c<-data.frame(time="3500",x=rnorm(50,12,4))     d<-data.frame(time="5000",x=rnorm(7,14,5))     e<-data.frame(time="7000",x=rnorm(20,10,3))     f<-data.frame(time="7500",x=rnorm(15,11,3))     g<-data.frame(time="9000",x=rnorm(15,10,5))     h<-data.frame(time="9500",x=rnorm(35,30,2))     i<-data.frame(time="10000",x=rnorm(30,28,4))     a2i<-rbind(a,b,c,d,e,f,g,h,i)           library(ggplot2)     a2i$time<-as.numeric(levels(a2i$time))[a2i$time]      ggplot(a2i,aes(time,x))+stat_smooth()+geom_point()      ![enter image description here](http://i.stack.imgur.com/cASXS.png) Here's what I'd be most grateful for some advice on... Q1. Would it be valid to do the Change Point Analysis on a vector like the means or medians of the groups? That would allow me to start with a 'one data point per unit of time' input format which would suit the `R` packages, as I understand them. I've seen it done with environmental data like monthly gas concentrations (from daily observations), but I thought I'd check. Q2. Is there a kind of Change Point Analysis that I can do on the raw data in `a2g` that will give me some measures of the probabilities of changes across the sequence? For example, something that will detect the change from time=9000 to time=9500, using all the data points in the sample? I'm guessing that if it was possible, someone would already have implemented and I just need a pointer to the relevant function. Q3. In case Q2 can be answered 'yes', would the method change if the distribution of each group's values was non-normal (unlike my sample data)? Q4. If Change Point Analysis is completely the wrong approach here, please let me know. I'm basically just curious about methods other than `ANOVA` for these kinds of data. Any other suggestions would be most welcome.


title: R examples for Durbin & Koopman "Time Series Analysis by State Space Methods"

text: I was wondering if anyone has ported the examples from Durbin & Koopman "Time Series Analysis by State Space Methods" to R? You can find RATS code for the examples online and obviosly SsfPack/Ox but no signs of an R companion for this book...


title: How to interpret a three-way ANOVA?

text: I have used a three-way ANOVA to analyze the the effect of genes, transcription factors, and different conditions on the gene expression. Now I have 9 elements, SSa, SSb, SSc, SSab, SSbc, SSac, SSabc, SSe. How can I interpret these values? I appreciate if you introduce me some resources, or share your own experience with interpreting 3-way ANOVA. Best


title: Visualising relationship data

text: I have data about a certain workforce, with different worker types W1, W2, W3, W4, etc. I need a method to visualise the relationship between each worker type. In particular, I'm interested in visualising the proportion of each worker type (e.g., W1, W2, W3, etc.) that works with each other worker type. Essentially, assuming there are 10 worker types, I'm looking at a 10x10 matrix. I could visualise it like these: ![http://i.stack.imgur.com/ZLIoT.png](http://i.stack.imgur.com/sI7WM.png) ![http://i.stack.imgur.com/CvrLq.png](http://i.stack.imgur.com/BfmK0.png) ...but I'm open to any other ideas. Those two diagrams require programs that I find difficult to use. I'm looking for something more appealing than a matrix with coloured cells.


title: Clustering with shape prior

text: **EDIT STARTS** After seeing the comments and answers, I believe I started in a wrong direction. I have a set of rectangles, which I want to cluster as shown below. ![enter image description here](http://i.stack.imgur.com/V5vOs.png). The approach I took was to consider central points of each rectangle as a data point in $R^2$ and cluster them using Euclidean distance (K-means, K-mediods approach or any other method). I am **not** trying to discover "shapes" in the data, as I know the best shape would be rectangles. This is the reason I was trying to know if there is any way to include "shape priors" in the clustering methods. There might be algorithmic formulations of this problem, I would be glad if someone can point me to. I hope it is clear now what I want to do, I will be glad if someone can point me to right directions. **EDIT ENDS** Suppose, we have a set of points in $R^d$ which we want to cluster. We do have some idea about the shape of the clusters, i.e., we know the clusters should be spherical, ellipsoidal or rectangular. Is there a way to include this "shape prior" in the clustering process? (it would be more useful for me if this has been done in some well known algorithms like k-means) I did find some papers like this one and this one which talks about imposing some constraint into clustering, but the constraints are primarily in the form of "pairwise constraints" i.e. specifying whether two instances should be in same or different clusters.


title: Why is controlling FDR less stringent than controlling FWER?

text: I have read that controlling FDR is less stringent than controlling FWER, such as in Wikipedia: > FDR controlling procedures exert a less stringent control over false > discovery compared to familywise error rate (FWER) procedures (such as the > Bonferroni correction). This increases power at the cost of increasing the > rate of type I errors, i.e., rejecting the null hypothesis of no effect when > it should be accepted. But I was wondering how it is shown to be true mathematically? Is there some relation between FDR and FWER? Thanks and regards!


title: How can one perform a two-group binomial power analysis without using normal approximations?

text: I would like to do power analyses for hypothesis tests of (non-)equality of proportions in which the proportions are _very_ small. I would like to do so **without** using normal (or Poisson) approximations of the binomial distribution. There are several general types of power questions I'd like to be able to address.   1. _Post-hoc_ : Given $\Pr_1$ (probability of a success in group 1) and $\Pr_2$ and $N_1$ (sample size group 1) and $N_2$ to calculate the power of the design given $\alpha$.   2. _A priori_ solve for $N$ given $\alpha$, the ratio $N_1\over{N_2}$, $1 - \beta$ (power), $\alpha$, $\Pr_1$, and an expected $\Pr_2$   3. _A priori_ solve for $1 - \beta$ given $\alpha, N_1, N_2, \Pr_1$, and $\Pr_2$. An ideal response would include R code and point out any other givens that I forgot to point out. A simulation approach is not a suitable response due to the small proportions. With your solution, please also mention what kind of statistical test it is applicable to.


title: What are typical values to use for alpha and beta in Latent Dirichlet Allocation?

text: Specifically in the case where I don't know anything about the documents I'm working with. I'm looking a specific number or number range.


title: How to compute the confidence interval if I have 2 samples?

text: In the field of biology we often have $>1$ samples. For instance I can conduct two (but not many) experiments (replicates) each on a number of animals ($n_1, n_2$). The aim of the experiment is to estimate the proportion of the animals which prefer a certain habitat area. And I could estimate from each sample the proportion of those animals staying in that area: $p_1$ and $p_2$. Therefore, the sample mean of the proportion is: mean($p_1,p_2$). My question is: based on these two samples, how shall I present the confidence intervals of the sample proportions? Shall I present two CIs or is there anyway to compute one CI just based on these two samples? Animals have correlated intra-group behavior, so I can not simply combine these two samples into one and then use one of the standard method to compute CI. On the other hand, if I have a large number of samples, I could compute the CI directly from these samples. The method I could think of is to conduct bootstrapping 500 times from sample 1 and 500 times from sample 2, and then compute the CI from these 1000 random samples. Do you have any nice suggestions?


title: Inverse Mills ratio after OLS

text: **Short version of the question:** Is it possible to create a dependent variable in the first step of the Heckman Selection model such that it is possible to obtain the values for the calculation of the Inverse Mills Ratio for the second step by OLS? If so, are there any pitfalls regarding it? **Long version:** A two-step selection model is used, e.g., in the case that due to a selection process a significant part of the outcome variable $y$ is zero. To estimate unbiased coefficients of the independent variables, Heckman (1979) proposed a two-step estimator, where in the first step a Probit regression of a dummy indicating that $y>0$ on explanatory variables. Then the inverse mills ratio of every observation is calculated and included in a second-step OLS regression of the observations with $y>0$ of $y$ on explanatory variables. The inverse Mills ratio is the ratio of the probability density function and the cumulative density function of the normal distribution evaluated at the predicted outcomes $x*b_2$ devided by the standard error of the Probit estimation (normally the standard error is assumed to be one). The Selection model is the following (based on Wooldridge): $$ \begin{aligned} \text{Equation 1:}~~~~~~~~~~~~~~ y1 &= x_1b_1+u_1 \\\ \text{Equation 2:}~~~~~~~~~~~~~~ y2 &= 1[xb_2+v_2>0] \end{aligned} $$ where $x_1$ is contained in $x$. It is assumed that $E(u_1|v_2)=c_1*v_2$, i.e. we assume linearity. Moreover, $v_2\sim\mathcal N(0,1)$ and $(u_1,v_2)$ are independent of $x$ with zero mean. $(x,y_2)$ are always observed, $y_1>0$ only if $y_2=1$. $y_1>0$ only for $y_2=1$ Running OLS on the sample, where $y_1>0$ and thus $y_2=1$ yields biased estimates because of the correlation of the error terms. One has to include a correction for this: $$ E(y_1|x_1,y_2)=x_1*b_1+c_1*E(v_2|x,y_2)=x_1*b_1+c_1*h(x,y_2) $$ where $h(x,y_2)=E(v_2|x,y_2)$. In the censored sample where $y_2=1$, so we have: $$ E(y_1|x_1,1)=x_1*b_1+c_1*E(v_2|x,y_2=1)=x1*b_1+c_1*h(x,1) $$ As $E(v_2|x,y_2=1)=E(v_2|v_2>-b_2x_2)$=inverse Mills ratio (IMR) of $(b_2*x_2)$, we have: $$ E(y_1|x_1,1)=x_1*b_1+c_1*IMR(b_2*x_2) $$ as a corrected regression of the censored sample. Estimates for $b_2$ can be obtained by running a Probit regression of a dummy indicating that $y_1>0$ on independent variables $x_2$. As the Probit estimates are the inverse of the standard normal cumulative, $x_1*b_2$ of the Probit estimates have the desired feature that e.g. for the predicted probability 0.5 of $y_2=1$ half of $(xb_2+v_2)$ are below and half are above zero and thus the error term is positive for half of these observations. Why doesn't one just create a variable $y_2$ that is above zero for $y_1>0$ and below zero otherwise and runs an OLS regression on this variable as the first step regression and construct the Inverse Mills Ratio using the obtained estimates? If it is possible, what features does this variable need to have? An advantage of this procedure would be the possibility to estimate a fixed effects model in the first step with panel data by time-demeaning all observations. **Literature:**   * Heckman, James (1979): "Sample Selection as a Specification Error". Econometrica, 47 (1), 153–161   * Wooldridge: Econometric Analysis of Cross-Section and Panel Data


title: How to relate raw scores on five-point Likert to three conceptual categories?

text: My thesis is on employee productivity and I had a questionnaire which consists of a five-point Likert scale (5-strongly agree, 4-agree, 3-agree, 2-disagree, and 1-strongly disagree). On the other hand, I had a conceptual Framework of a 3 point scale (Highly Productive, Productive and Not Productive). **Is there any standard way to identify whether my respondents are highly productive, productive or not?**


title: I need a model that can predict based on multiple variables. How do I get started?

text: I have a problem where I have to predict a variable X that is dependent on several other variables a,b,c,d... I have the data containing the values of these variables a,b,c,d.. and also X up to a certain point in time T. I do not have an explicit relationship between a,b,c,d... and X but I have managed to fit the dependencies using a linear regression model. My goal is to predict X into the future beyond time point T. How do I approach this?. I had considered auto regression but then I would have to disregard the variables a,b,c,d.. as it would be based solely upon the past values of X. Would this cause a problem in accuracy? I had also thought of using auto regression separately on the variables a,b,c,d... and then using those projections and the original linear regression model to predict X at various points in time. Which would be a better model? Also is there any other way to implement this? And what are some good python libraries that I can use for these models and any you may suggest? Any help would be appreciated. Thanks!


title: What's a time series model for forecasting a percentage bound by (0,1)?

text: This must come up---the forecasting of things that are stuck between 0 and 1. In my series, I suspect an auto-regression component, and also a mean- reverting component, so I want something that I can interpret like an ARIMA--- but I don't want it to shoot off to 1000% in the future. Do you just use an ARIMA model as the parameter in a logistic regression to confine the result between 0 and 1? Or I learned here that Beta regressions are more appropriate for (0,1) data. How would I apply this to a time series? Are there good R packages or Matlab functions that make fitting and forecasting this easy?


title: How to overcome singularity problem in Linear Discriminant Analysis?

text: I've code for LDA which is failing as the matrices passed to calculate eigen values are not singular and hence lead to infinite eigen values? Can anyone recommend what can be done to fix this? I've tried is to compute rank of matrix and compare that with the minimum of the size of matrix to figure out singularity. What else can be done to overcome this problem?


title: Find the probability distribution of X?

text: Suppose Nokia store places 20 of its cell phones on a clearance sale, unknown to anyone 5 of these cell phones are defective. A customer selects 3 cell phones at random for inspection. Let X be the number of defective cell phones in the sample. Find the probability distribution of X? My attempt: $$ \begin{matrix} x & 0 & 1 & 2 & 3 \\\ P(X=x) & 0.015625 & 0.140625 & 0.421875 & 0.421875 \\\ \\\ \end{matrix} $$ I calculated the values using the binomial distribution: $$ P(X = x) = nCx \ p^{n-x} (1-p)^x$$ Where, $n$=5, $x$=0,1,2,3 and $p$=(.25) [from 5/20]


title: Time series modeling with high-frequency data

text: I'm looking for some forecasting advice when dealing with seasonal time series data that has a large number of observations. By "large" I only mean a few thousand --- I'm used to such sizes in Data Mining being considered pretty small, but it seems that in time series modeling that's pretty unwieldy for many of the tools I've tried. For example, here's a toy data set that records an observation once per minute, for five days:               set.seed(123)     t <- 1:(5*24*60)     x <- ts(15 + 0.001*t + 10*sin(2*pi*t/(length(t)/5)) + rnorm(length(t)), freq=length(t)/5)     plot(x, type='l')      ![time series plot](http://i.stack.imgur.com/xVSCN.png) (In my real operational data set, the values are observed at irregular intervals, but I've regularized them by doing something like `x <- approx(d$t, d$x, xout=1:(5*24*60))` first. Advice on whether that's advisable, or alternative approaches, is welcome too.) So the seasonality in this data set has a lag of 1,440 observations, which seems to be way outside the range that things like `auto.arima()` (in the `forecast` package) will find:               m1 <- auto.arima(x)     plot(forecast(m1))      ![prediction plot](http://i.stack.imgur.com/ccnGc.png) And I'm not quite sure how to interpret the `ets()` function here, but it doesn't seem to be able to handle this size data, and it didn't seem to pick up on the seasonality:               > m2 <- ets(x, 'MAZ')     > plot(forecast(m2))     Error in forecast.ets(m2) : Forecast horizon out of bounds     > m2$method     [1] "ETS(M,A,N)"      Where to go from here? Any suggestions? Thanks.


title: Experiment or simulation to undestand type I and type II errors

text: I am studying Type I and Type II errors and basic concepts of testing hypotheses. To better develop my intuition I would like to write a few simple simulations. What I'm hoping is to get a base simulation going, that I can then tweak to understand this area. I am not looking for code in a specific language. However, pseudo code would be great, or references to a well documented example. Thanks.


title: Is there an easy way to calculate significant difference between two largely overlapping correlations from same sample?

text: I am comparing three different measures comprised of 10-20 items each. Each of these measures is actually a short form of a longer 40-item measure and each of the short forms shares some similar items, and thus they are quite similar. I am interested in seeing how each of these short forms are correlated with the same dependent variable. I come up with _r_ = .54, .56, and .58. How can I assess whether these correlations are significantly different? It seems doing a simple _r_ -to- _z_ transformation would not be appropriate since the independent variables are quite similar and the dependent variable is in fact the same for all three correlations. As such, all of these *r*s overlap quite a bit. I found this one paper: Zou, G. (2007). Toward using confidence intervals to compare correlations. _Psychological Methods_ , 12(4), 399. But it seems way over my head and I am wondering whether there is an easier way to do it... Any guidance would be appreciated!


title: Estimating a 1-D Brownian motion process using noisy observations

text: This question is a follow-up on my previous question. Suppose I have a Brownian motion process that is defined as follows: at time $i=1$ random variable $X_1\sim\mathbf{N}(\mu,\sigma^2)$, and, for any subsequent random variables $X_i$ and $X_j$ such that $i<j$, $X_j- X_i\sim\mathcal{N}(0,(j-i)\sigma^2)$. I collect $n$ noisy observations $\\{Y_i\\}_{i=1}^n$ of $\\{X_i\\}_{i=1}^n$, where $Y_i\sim\mathcal{N}(X_i,N_0)$. I would like to construct $\\{\hat{X}_i\\}_{i=1}^n$ such that the squared error $\sum_{i=1}^n(X_i-\hat{X}_i)^2$ is minimized. Unlike in my previous question, the causality is not a constraint here (I have access to all $n$ observations and am interested in the "evolution" of the Brownian motion random process over the entire window of observation.) I am wondering if it is possible to obtain closed-form expressions for $\\{\hat{X}_i\\}_{i=1}^n$ from $\\{Y_i\\}_{i=1}^n$ as well as the squared error $\sum_{i=1}^n(X_i-\hat{X}_i)^2$ given $\sigma^2$, $N_0$, and $n$. Note that this problem is easily convertable into the continuous domain, which, as I have been told, admits solutions using Wiener filter. Since my problem (converted to that domain) requires a non-causal filter, the Wiener filter seems very simple. However, use of the Wiener filter requires that the signal is stationary and, in my case, Brownian motion is not stationary, as one can readily compute: $$\begin{array}{rcl}E[X_iX_j]&=&\int\int x_i x_j p(x_i)p(x_j|x_i)dx_jdx_i\\\ &=&\int\int x_i x_j \phi(x_i;\mu,i\sigma^2)\phi(x_j;x_i,(j-i)\sigma^2)dx_jdx_i\\\ &=&\int x_i^2\phi(x_i;\mu,i\sigma^2)dx_i\\\ &=&i\sigma^2 \end{array}$$ showing that autocorrelation actually depends on $i$ as opposed to $j-i$ (where $\phi(x;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$ is the pdf of the Gaussian distribution). Does anyone know where to go from here? Are there closed-form expressions for $\\{\hat{X}_i\\}_{i=1}^n$ from $\\{Y_i\\}_{i=1}^n$ as well as the squared error $\sum_{i=1}^n(X_i-\hat{X}_i)^2$ given $\sigma^2$, $N_0$, and $n$? This seems to me like a problem that someone thought about before. Unfortunately, I am not well-versed in filtering/signal-processing -- perhaps there is a version of a Wiener filter that could work. I would appreciate any guidance in figuring this out. **EDIT** : This may be soluble using the first principles from which the Wiener filter is derived, i.e. via the orthogonality principle and using the expectation operator as a dot-product. I'm fearful that this might get me stuck with a nasty integral equation, requiring the use of transforms (which I am not too familiar with to begin with) but I'll try it in the morning.


title: Backpropagation for convolutional network

text: I'm trying to develop the equations for back propagation for convolutional network. Let's say I have a network with one convolutional network (no pooling) with no activation function. The activation function of the layer is $f = x*k$ where $x$ is the input matrix and $k$ is my convolution kernel. If I derive $f$ by $k$: $\frac{df}{dk} = x * \frac{dk}{dk}$ If $k$ is a vector in size $N{\times}1$ then $\frac{dk}{dk}$ is a unit matrix $N{\times}N$. This derivative doesn't fit in size with what I expect and I can't propagate it further. I assume I'm doing something wrong... Thanks.


title: How to create an index with only categorical variables (dichotomous, nominal and ordinal) only?

text: I am currently conducting a study about the involvement of male in family planning. It is of interest to create an index for the involvement of male in FP. My variables were only categorical variables and mainly dichotomous variables. Is there anyone who is knowledgeable on creating index using logistic approach? or any suggestions for the creation of this index. Thanks so much


title: Finding the fewest needed samples for a regression method

text: I have a regression method that I have applied with success to tens of thousands of well-observed objects. I'm looking to estimate how well it works on poorly observed objects. I selected 50 objects with $R^2$ values exceeding $0.999$ that have each been observed at least 100 times and now I'm trying to sample small numbers of points from each of them and run the method to see how it fares. I want to know the distribution of $R^2$ as a function of sample size. For each object, for each sample size N=10..50 or so, I'm grabbing 100 random samples of size N, training the method on the sample, and then testing on the left out points, and computing the $R^2$ score from that. The result is that I have a table containing 50 objects, 40 sample sizes per object, and 100 $R^2$ scores per sample size per object, i.e.:               SampleSize_10 Object_1 R2_1     SampleSize_10 Object_1 R2_2     ... (skipping R2_3 through R2_99)     SampleSize_10 Object_1 R2_100     SampleSize_10 Object_2 R2_1     ... (skipping Object_2 R2_2 through Object_50 R2_99)     SampleSize_10 Object_50 R2_100     SampleSize_11 Object_1 R2_1     ... (skipping SampleSize_11 Object_1 R2_2 through SampleSize_50 Object_50 R2_99)     SampleSize_50 Object_50 R2_100      My question is: how should I collate this information? I was thinking that for each sample size, for each object I could take the mean of all of the $R^2$ scores for that object at that sample size, and then take the mean of all of the means from all of the objects at that sample size. But since $R^2$ is unbounded from below, I was thinking maybe I should use median instead of mean to control extreme outliers, although I'm afraid that might be disingenuous. Also, how can I get confidence intervals around the resulting number? I think the distribution would likely be normal due to CLT so I could probably get confidence intervals from STD but I suppose that's something I should probably check?


title: Analysis of sells (what to buy)

text: Statistical analysis is not my field, so forgive me if this is a very stupid question. I have a list of sold items by size. Shoes in this case               Size    Qty     35  2     36  1     37  4     38  4     39  32     40  17     41  23     42  57     43  95     44  90     45  98     46  33     47  16     48  4     total:  476      I have to tell the owner how much of every size to buy. The problem is, I can't say him. \- You should buy 95 shoes size 43 for every one of size 36... The usual practice is to buy the whole size curve and buy extras for the most selling sizes.   This is about a year worth of data. How should I present this information in an easy to understand way? What I want to present is a general rule. Something like "for every size curve, you should buy x additional shoes of size x".   The idea would be to later apply this approach to other clothing items.


title: SVM rbf kernel - heuristic method for estimating gamma

text: I read on this exchange a heuristic method of estimating gamma for the rbf kernel in SVMs. I was wondering if someone might be able to explain it to me in a little more detail? I believe you select 1000 (or a large number) of pairs of datapoints from the dataset then calculate the norm for the difference of each pair. Apparently, the inverse of the .1, .9 quantiles and median are good candidates for a suitable gamma for the rbf kernel. Thanks


title: how do I "build" a distribution from a list of values and then calculate p-values?

text: I have a list of values (A) that I consider to be the null distribution which I want to use to test the significance of another set of values (B) and get p-value for each item in B. I read this post Calculating p-value from an arbitrary distribution but I still don't know how to do it in R or other platform... can someone give me some links or pointers on how to do it? Edit: the question I'm trying to answer is given a list of genes, whose expression are most variable among a group of people. You can imagine I have a big matrix with each row representing a gene and each column representing a sample. So I have many genes and many samples. I want to find out those genes that show variation (I don't know the variation characteristics yet). I have now a list of genes who I defined as invariable (basically all of their values are above 0.8 for example) I thought I can calculate the variance for these genes and consider them as a null distribution. And then I can use this null distribution to test my genes of interest on their variance and get a p-value for each gene's variance. So in the end, I thought i could get a list of genes with significant p-value (or corrected p-value with some multiple testing correction).


title: Rejection sampling from a Gamma distribution using a Cauchy proposal

text: i'm trying to find the parameters $ \gamma,x_0$ of a standard Cauchy distribution : $$T(x)= \frac{1}{(\pi \gamma (1+(\frac{x-x_0}{\gamma})^2))} $$ To perform rejection sampling from a gamma distribution: $$ G(x) = \frac{1}{\Gamma(\alpha)}x^{\alpha-1}\lambda^\alpha \exp(-\lambda x) $$ And i have seen than in several textbooks as Bishop's pattern recognition it seems like if is straightforward to derive those parameters to obtain the values $$x_0 = \frac{\alpha -1}{\lambda}$$ and $$\gamma^2 = 2\alpha-1$$ But so far i just have obtained $x_0$ by making equal the derivatives of T and G and solving for x, but for the $\gamma$ parameter i just dont know how to find it, because all the attempts that i have made leads me to polynomial equations of degree >=3 and irrational roots... Thanks for your help!


title: Does the more complex model necessarily have a lower deviance?

text: I am aiming to compare the deviances of two models, an ordered probit and a multinomial probit, using the likelihood ratio test (obviously using the same data). However, I systematically obtain a higher deviance for the most complex model, i.e. the multinomial probit. I don't understand this result. I am using Matlab to perform the analysis, but as I try and find possible mistakes with the code I would also like to know from a statistical viewpoint if such a result if possible, and if it is, why?


title: How to estimate False Discovery Rate from p-value distribution?

text: I have learned many models and I calculated p-values for the cross-validation errors. I want to select significant models based on the false discovery rate (FDR). How can I estimate the FDR from p-value distribution?


title: Random sorting of an array of weighted scores

text: I'm trying to construct an algorithm that randomly sorts an array according to a weighted score for each item. In the most straightforward version of the algorithm, the output array can be constructed iteratively, where the next item added is chosen randomly according to the distribution of the weights of the remaining items. To illustrate with an example, if we have some items and weights like this: (code is PHP but that's not really important)               $items = array('a' => 5, 'b' => 4, 'c' => 1);     $output = weighted_random_sort($items);      Then in \$output, 'a' should be first with probability 0.5, 'b' should be first with probability 0.4 and 'c' should be first with probability 0.1. The overall distribution of $output should be:   * P('a', 'b', 'c') = (5/10)*(4/5)   * P('a', 'c', 'b') = (5/10)*(1/5)   * P('b', 'a', 'c') = (4/10)*(5/6)   * P('b', 'c', 'a') = (4/10)*(1/6)   * P('c', 'a', 'b') = (1/10)*(5/9)   * P('c', 'b', 'a') = (1/10)*(4/9) This algorithm is easy to implement with the iterative algorithm I described above. However I actually want to implement it in a different way, and here's where the statistics part of my question comes in. I want to implement it by perturbing each weight with a random function, and then just doing a plain reverse-sort of the perturbed scores. So basically, my question is, what kind of random function f(n) would allow me to do:               $items = array('a' => f(5), 'b' => f(4), 'c' => f(1));     $output = arsort($items)      and allow \$output to have same probability distribution as the original algorithm?


title: r :Does rpart require normalized data set for the input variable?

text: Does rpart function require normalized data when using the "anova" method? Im assuming so, but I do not like to assume. I have looked and no where has it said it needs to be, but the method makes me think it does. Thanks Scott


title: Why R function Chol() breaks for Hilbert matrices of order n > 12

text: Is there any theoretical reason that Chol() in R breaks for Hilbert matrix of order greater than 12? Thanks,


title: Calculate the tendency of a set of samples

text: I develop an application in which i constantly get samples of heart pulse. I defined an interval of t seconds. In each t seconds I have n samples. In every interval, I want to calculate the tendency of those n samples. For example, lets say I have n = 5, and I have samples with values {70, 88, 95, 103, 115}. I want to recognize I have a growth in the heart pulse, and I want to have some measure for the rate of grows/decrease/almost no change. I thought on two approaches for solving this problem.   1. I calculate the linear approximation for the n samples using linear regression by applying least squares implemented by the normal equations. (I treat each sample as two coordinates vector with x coordinate as time and y coordinate as heart pulse). I get from the normal equations a linear function of the form y = mx+b and my measure for the tendency is the slope, i.e. the m value.   2. The second approach is to calculate the correlation (pearson's correlation) between vector x and vector y when x is the time and y is the heart pulse. I'm asking which approach do you think is better for my problem (determining the tendency of heart pulses). Or, do you have a better algorithm for solving this problem?


title: Poisson Repeated Measures ANOVA

text: The following is a general background into my data analysis problem:   * I am working with count data that is non-normal.   * Data was collected at 2 time periods.   * I am interested in comparisons of species counts between treatments.   * There are 4 treatments. Data were initially analyzed using the Kruskal-Wallis rank sum test (a non- parametric ANOVA counterpart). A reviewer of my manuscript suggested using a Poisson distribution. How can I perform a repeated measures ANOVA with a Poisson distribution in R?


title: Is the restricted Boltzmann machine a type of graphical model (Bayesian network)?

text: In general, neural networks are not graphical models due to the choice of the cost function. But are restricted Boltzmann machines a special type of Bayesian networks?


title: Multiple resampling test/train dataset when choosing new models?

text: I have been reading several posts on testing multiple models on the same dataset, which can lead to problems controling type-1 errors. Mostly these posts have to do with data-mining on big datasets: How to draw valid conclusions from big-data and ensuring testing data doesn't influence training However, I know that this is done frequently. In fact, in my old work it was my job to find the best (logit) model given a dataset. In order to decide if a model is predictive or not, you have to validate against the test data. If the performance is poor, then you start from zero and create a new model. By the end you may have checked the test dataset dozens of times. I was recently asked if multiple re-sampling of the dataset would be a possible solution. I wanted to say 'no,' but I don't actually know _why_ this is bad. To give a specific example: Suppose I am looking for the best linear regression given a dataset of 1,000 observations. I split the data into training and testing sets. I formulate a model which is not satisfactory on the testing sample. So, I _redistribute_ the 1,000 observations into new training/test samples and attempt to find a new model. Each model will be trained and tested on their own specific training/test sample, which all come from the same 1,000 original observations. **My question: Why is this incorrect? What are the problems that are created with this methodology?**


title: Choice of K in K-Fold cross validation

text: I've been using the _K-Fold cross validation_ a few times now to evaluate performance of some learning algorithms, but I've always been puzzled as to how I should choose the value of K. I've often seen and used a value of `K = 10`, but this seems totally arbitrary to me, and I now just use 10 by habit instead of thinking it over. To me it seems that you're getting a better granularity as you improve the value of K, so ideally you should make your K very large, but there is also a risk to be biased. I'd like to know on what the value of K should depend, and how I should be thinking about this when I evaluate my algorithm. Does it change something if I use the _stratified_ version of the cross validation or not?


title: How to test hypothesis on a fraction of two coefficients?

text: Say you have an equation: $$y = b_0 + b_1X_1 +b_2X_2$$ How would you test the hypothesis that $\frac{b1}{b2} = q$?


title: Vector autoregression - number of variables to use

text: When building a vector autoregression model is there some theory that would guide me in chosing the number of variable to include? For example, I have about 3000 data points and I would like to get an idea of how many lags of explanatory variables to ues.


title: PP-plots vs. QQ-plots

text: What is the difference between probability plots, PP-plots and QQ-plots when trying to analyse a fitted distribution to data?


title: When should one consider using GMM?

text: One of the things which makes econometrics unique is the use of the Generalized Method of Moments technique. What types of problems make GMM more appropriate than other estimation techniques? What does using GMM buy you in terms of efficiency or reduced bias or more specific parameter estimation? Conversely, what do you lose by using GMM over MLE, etc.?


title: Multivariate linear mixed model in R

text: I have run into a problem with respect to an application of linear mixed effects model using lme4 package and I wondered if I could seek your help. This is my model in a multivariate setup where $Y_{ik}$ is the gene expression value for a give tissue type $k$ and $G$ is the genotype with $G_i \in \left(0,1,2\right)$ for an individual $i$. \begin{equation*} Y_{ik} = \beta_{0} + \beta_{k} + \beta_{1}G_{i} + b_{i} + b_{k}G_{i} + \epsilon_{ik} \end{equation*} where $\beta_{0}$, $\beta_{k}$ and $\beta_{1}$ are the fixed effects and \begin{equation*} b_{i} \sim N\left(0,\tau^2\right) \qquad b_{k} \sim N\left(0,\Psi^2\right) \qquad \epsilon_{ik} \sim N\left(0, \delta^2\right) \end{equation*} are the random effects. I am testing for the variance component $\Psi^2$, a 3 df test in my case since I am simulating 4 tissue types: \begin{equation*} H_0 : \Psi^2 = 0 \qquad H_A : \Psi^2 > 0 \end{equation*} In other words, are there any tissue specific gene expression differences influenced by the genotypes/SNPs? I am generating null simulations $\left(\Psi^2 = 0 \right)$ under the aforementioned conditions and here's my code:               require(lme4)     nobs = 150; k = 4; ngene = 1; maf = 0.3; sdB = 4; sdW = 2     alpha = 0.05     Geno = rbinom(nobs,2,maf)     IND = as.factor(rep(1:nobs,each=k));     Tissue = as.factor(rep(1:k,nobs)); data = data.frame(IND,Tissue)     data$Geno = Geno[data$IND]     beta0 = 8; beta1 = 5; betak = c(1,1,1);     beta = c(beta0, beta1, betak)     cMat = model.matrix(~ Geno + Tissue,data)     means = cMat %*% beta     bi = rnorm(nobs, mean = 0, sd = sdB)     sdBt = NULL     if(!is.null(sdBt)){             bk = c(0,rnorm(k-1, mean = 0, sd = sdBt))             bkG = Geno %*% bk             Y = t(apply(means, 1, function(m) rnorm(1, mean=m, sd=sdW))) + rep(bi,each=4) + as.numeric(unlist(split(bkG,1:NROW(bkG))))         }else{                                   Y = t(apply(means, 1, function(m) rnorm(1, mean=m, sd=sdW))) + rep(bi,each=4)         }     data = data.frame("IND"=IND, "Gene" = round(as.vector(Y),3), "Tissue" = Tissue, "Geno" = rep(Geno,each=4))     Fit = lmer(Gene ~ Geno + Tissue + (1|IND) + (0+Geno|Tissue),REML=T,data)     Fit0 = lmer(Gene ~ Geno + Tissue + (1|IND),REML=T,data)      lrt = anova(Fit,Fit0)     lrt[2,7] < alpha      I ran 1000 such null simulations generating one SNP and one gene at a time under varied samples sizes with 1000 being the max sample size and checked the type I error rate, which is unfortunately, never close to 0.05. Is there anything wrong with the code -- w.r.t. how I am generating the datasets and running the mixed effects model? I really appreciate any help. Thanks. **EDIT:** I tweaked my code a little bit and the ran null simulations only to get conservative estimates of the type I error (between 0.02 and 0.03). Here is the reason why -- true parameter being on the boundary of the parameter space It has been an incredible learning experience!


title: arima function in R, data= and subset= not available? How to specify training set?

text: I would like to fit my arima model on some subset of my full data, say, based on observations within a certain time range. Normally, say in fitting an `lm()` object, I would use the `subset=` argument with some logical vector, but this does not exist for the `arima` function (nor does it exist in other packages like `forecast` or `TSA`). Any suggestions on how to restrict my model to be fit on observations within a certain time range? I can think of two possibilities, but neither are particularly good. When I specify my response variable and my `xreg=` variables, I know I could wrap each variable I want to use within a `window()` statement. Alternatively, I could do a window statement on my whole dataframe to make a new dataframe, and then use               arima(ShortDataframe$var1,           xreg=model.matrix(~ ShortDataframe$var2 + ShortDataframe$var3 + ShortDataframe$var4)      but I am hoping for a solution that is user-friendly.


title: How to measure pre- and post- intervention effect with single metric capturing percent and absolute magnitude of change?

text: My colleagues and I are attempting to best quantify the benefits of an intervention program on the number of emergency department (ED) visits per patient during pre-intervention and post-intervention periods (Year 1 and Year 2). We're working with a health care dataset, but in principle my question could apply to a number of areas. Simply creating a variable that measures the % change in annual ED visits from Year 1 to Year 2 ignores the absolute difference in the two ED rates. For example, a reduction of 100 ED visits during Year 1 to 50 visits during Year 2, and a reduction of 2 ED visits to 1, are both coded as a 50% reduction. Meanwhile, measuring the absolute change in ED visits captures the magnitude of change but ignores the change relative to the patient's baseline # ED visits in Year 1. A drop of 1 ED visit from 100 to 99, or 2 to 1, are both coded as -1. We want to capture _both_ the % change and the absolute difference in a single quantity measuring the effectiveness of the intervention. Realizing this is a somewhat subjective question - are there recommendations how best to construct the response variable (which will be modeled with a multivariate regression)? Thanks!


title: Multiple regression in directional / circular statistics?

text: I'm trying to develop a predictive model for an angular dependent variable (on $[0,2\pi])$ using several independent measurements – also angular variables, on $[0,2\pi]$ – as predictors. Each predictor is significantly but not extremely strongly correlated with the dependent variable. How can I combine the predictors to determine a predictive model for the dependent variable that is optimal in some sense? And how can I rigorously identify the strongest predictor(s)? For variables on Euclidean space(s), I'd employ multiple regression (or similar) and principal components analysis. But the periodicity of all variables mucks with these approaches, e.g., 0.02 must be highly correlated with 6.26, but not with 3.14. How are "the usual" procedures generalized to directional / circular statistics? Any insights or cites to useful references would be useful. (I'm already aware of the texts by N. Fisher and Mardia & Jupp, but don't have handy access to these.)


title: How to fit this neuron firing model with R?

text: I originally posted this as an answer elsewhere but in retrospect it seems more like a question: What is the sample-size range for which the median should be preferred to the mean as a measure of central tendency? And why? I have R code that simulated neuron activity based on the model described at these links: http://www.izhikevich.org/publications/spikes.htm http://www.izhikevich.org/publications/spikes.pdf I know how to simulate one but not fit the parameters to existing data. Take this pattern which was generated with the code below. Each neuron could have different a, b, c, d, k1, k2, k3 parameters and also the I parameter (Input Current) could vary over time or be constant for each neuron. The model is described in the pdf above. Anyone know how to fit this model using R? What I would like to be able to do is fit a set of parameters to actual firing data (which would be voltage over time, same as what is simulated). ![enter image description here](http://i.stack.imgur.com/AGp2b.png) **R code to simulate Neuron:**               a=.1; b=.2; c=-65; d=2     k1=.04; k2=5; k3=140     I=4          Ntime=1000     v<-matrix(nrow=Ntime)     u<-matrix(nrow=Ntime)          v[1]=-65     u[1]=0     for(i in 1:Ntime){       #Fire Neuron if Voltage>30         if(v[i]>=30){           v[i]<-c           u[i]<-u[i]+d         }              #Calculate Voltage Change         dv= k1*v[i]^2 + k2*v[i] + k3 - u[i] + I         du= a*(b*v[i] -u[i])              #Update Voltage         v[i+1]<-v[i]+dv         u[i+1]<-u[i]+du          }          plot(v, type="l", xlab="Time (ms)", ylab="Voltage")


title: Imputation before or after splitting into train and test?

text: I have a data set with N ~ 5000 and about 1/2 missing on at least one important variable. The main analytic method will be Cox proportional hazards. I plan to use multiple imputation. I will also be splitting into a train and test set. Should I split the data and then impute separately, or impute and then split? If it matters, I will be using `PROC MI` in `SAS`.


title: Is there a systematic study of the effects of choosing the wrong distribution to model data?

text: I would like to know if there is any study showing what are the effects of modelling data using the wrong probability distribution. Any reference to papers or books explaining it (possibly for a novice/intermediate audience) is fine.


title: Finding optimal parameter values using a Bayesian model

text: I have a problem with the following setup. I've been reading "Doing Bayesian Data Analysis: A Tutorial with R and BUGS" and it seems like the Bayesian approach is a good one, but I'm not entirely sure how to model it. I've got roughly 20 events that I'd like to analyze. Each event should be similar, and I can apply the same analysis to each. I'll assume that the events are independent (they might not be, but I'm not worried about that). For each event, I have a function of three parameters: threshold, s, and t. Threshold is a continuous variable between 0 and .5, and I'm guessing that the useful values are around .1. s is an integer between 15 and 2, and t is an integer between s and 1. I can run the function over all values of s and t and a sample of values for threshold, probably .05, .06, ..., .2 (or whatever, this is flexible). The function gives me a raw score, where higher is better, and negative values are, well, negative; they represent an actual loss. It will be straightforward enough to normalize the raw score to a value between 1 and -1. What I want to know is what are the best values for the parameters (or, technically the distributions of the parameters) that will lead to the best(ish) score(s), after looking at all the events. What I'm thinking I'll do is the following: Model threshold as a normal-like variable - using either a normal distribution, or beta, or gamma. I'm leaning towards a beta or gamma distribution because it is fairly easy to model my high uncertainty about the prior, and I don't expect the distribution to be perfectly bell-shaped, or even symmetric. Similarly, model s and t as normal-like distributions (again, leaning towards beta or gamma, for the same reasons). Even though the values are integers, it's no big deal at the end of the day if the mean of the posterior distribution is, say 8.6 - I can work with that. If there were distributions that handled discrete values, that would be cool, but, again, not necessary. I realize the binomial distribution is discrete, but it's more about modelling binary outcomes, if my understanding is correct, so I'm going to steer away from it in this particular instance. My question is, how do I utilize the score, either in raw form (guessing not that useful) or the normalized form (hopefully useful)? It seems like the most important piece of information, but I don't know how to incorporate it. I'm not trying to build a regression model; I'm looking more for the optimal results. And I don't want to do a logistic regression; again, it's not a prediction, but also because it's not a simple yes or no - best score or not - I want to weight parameters that give a normalized score of say .98 nearly as high as the ones that produce a 1. If I were doing this without Bayesian analysis, I think what I would do is something like, for each data point, for each parameter, add a value equal to normalized score for that data point. Then, at the end, look at, say, the top 3 or 5 values for each parameter, and try to make a decision based on that. But I'd like to understand how to model this in a Bayesian way. Any advice appreciated, and thanks in advance!


title: Specification of a mixed model with nesting with lmer

text: I've been told my model wasn't statistically correct, but I'm not sure what I should do. What I have is basically 2 **species** : A and B. For each species, I have 3 **replicates** : 1, 2, and 3 (=6 populations in total). I'm applying a **treatment** with 3 levels to each population (with different individuals being used for different treatment levels). And I'm testing the effect of the **species** , the **treatment** and of their **interaction** on the response variable. I therefore need to nest my replicated populations within 'species', as I am interested in the 'species' effect. **The model I wrote was:**               model <- lmer(response variable ~ treatment*species + (1|population), data=data)      with 'population' being the combination between 'species' and 'replicate', and therefore each species has 3 distinct levels of population (A1, A2, A3 and B1, B2, B3). **This should thus be equivalent to:**               model <- lmer(response variable ~ treatment*species + (1|species:replicate), data=data)      So I thought the **nesting** was OK like this. But I got the complaint that, because I have a split-splot design, I should include the treatment*replicate(species) interaction in the model, so that populations are treated as the appropriate unit of replication. Am I missing something, or is my model already correct with this lmer specification? As I have already nested replicate within species, I thought it was taken care of. Thanks!


title: Quantifying variable importance for GLMM using hierarchical partitioning (in R)

text: I am interested in quantifying variable importance for a binomial logistic mixed-model regression. My model has 5 fixed effects, and 3 random effects (2 nested). I am doing model selection and multimodel inference via the `MuMIn` package in R. I originally planned to use Akaike weights to calculate variable importance using the `MuMIn` package, but then I came across Murray & Connor 2009 (full text here) which shows that calculating the independent effect through hierarchical partitioning performs much better for their simulated models. However, their simulated models do not contain random effects, and the `hier.part` package which can perform which analysis, also does not have an option to include random effects. My question is: do the random effects need to be taken into account in calculating independent effects using hierarchical partitioning?


title: Performance comparison of “patternnet” and “newff” for binary classification in MATLAB R2014a

text: I have a binary classification problem for financial ratios and variables. When I use `newff` (with `trainlm` and `mse` and threshold of 0.5 for output) I have a high classification accuracy (5-fold cross validation – near 89-92%) but when I use `patternnet` (`trainscg` with `crossentropy`) my accuracy is 10% lower than `newff`. (I normalized data before insert it to network - `mapminmax` or `mapstd`) When I use these models for out-sample data (for current year- created models designed based one previous year(s) data sets) I have better classification accuracies in `patternnet` with better sensitivity and specificity. For example I have these results in my problem: # Newff: Accuracy: 92.8% sensitivity: 94.08% specificity: 91.62% Out sample results: accuracy: 60% sensitivity: 48% and specificity: 65.57% # Patternnet: Accuracy: 73.31% sensitivity: 69.85% specificity: 76.77% Out sample results: accuracy: 70% sensitivity: 62.79% and specificity: 73.77% Why we have these differences between `newff` and `patternent`. Which model should I use? Thanks.


title: Cross-validation prediction standard error

text: I am using leave-one-out cross-validation to evaluate a linear regression model. In subsequent analysis, I need three specific values for each observation: observed value, predicted value, prediction standard error. Prediction standard error values can be retrieved from function `predict.lm` setting argument `se.fit = TRUE`. The following code (adapted from here) can be used to do what I currently need:               library(faraway)     gala[1:3, ]     c1 <- c(1:30)     gala2 <- cbind(gala, c1)     gala2[1:3, ]     obs  <- numeric(30)     pred <- numeric(30)     se   <- numeric(30)     for (i in 1:30) {          model1  <- lm(Species ~ Endemics + Area + Elevation,                        subset = (c1 != i), data = gala2)          specpr  <- predict(model1, gala2[i, ], se.fit = TRUE)          obs[i]  <- gala2[i, 1]          pred[i] <- specpr$fit          se[i]   <- specpr$se.fit     }     res <- data.frame(obs, pred, se)     head(res)       obs       pred       se     1  58  70.185063 5.524249     2  31  72.942732 6.509655     3   3  -8.303608 7.055163     4  25  20.948932 6.998093     5   2 -15.953141 7.403062     6  18  27.274440 6.220029      I searched through the documentation of some of the packages that offer functions for cross-validation, but did not find any that saves prediction standard errors. Is there any package that already offers such functionality?


title: assumptions of REML

text: In Patterson and Thompson 1971, $L$ is the log-likelihood of $y$ $$ L = const - \frac{1}{2}log |H| - \frac{1}{2}n log (\sigma^2) - \frac{1}{2\sigma^2} (y - X\alpha)'H^{-1}(y-Xa) $$ They consider the transformations $Sy$ and $Qy$ with respective log- likelihoods $L'$ and $L''$ and the assumptions:   1. $S$ is of rank $n-t$ and $Q$ is of rank $t$   2. $Sy$ and $Qy$ are statistically independent From these assumptions it follows that $L=L' + L''$. Can someone give a formal proof? Of course 2. explains the additive part on the left, but the equality is not clear to me. Thanks


title: Visualizing Mutual Information Against TF-IDF for Text Corpus Data

text: I'm working on a data visualization project for the semester and have decided to work with a corpus of discussion forum data focused around debate over political issues (available here). I'm visualizing different frequency counts of unigrams, bigrams, and trigrams, but I was also curious about using different metrics for visualization such as TF-IDF and Mutual Information. One of my ideas for a visualization was to plot n-grams in a scatterplot showing their TF-IDF scores against their Mutual Information scores (Mutual Information being MI between an n-gram and a topic of debate within the corpus, say, "abortion"). My thought was to have each of these metrics as an axis and have each data point represent an n-gram within the plot. My questions are (a) is this informative at all, or is it completely statistically unsound or redundant to plot n-grams for these two metrics this way, and (b) if something like this would make for a decent visualization, is Mutual Information the best correlation metric to use for one of the axes? Would something else like Chi-square be more appropriate? Thanks.


title: Joint distribution of dependent exponential variables

text: I am trying to model a joint probability constraint of a union of correlated / uncorrelated exponential functions. Formally, assume we have random variables $A_i \sim exp(\theta_i), \forall i \in {1…n}$. The rate of each exponential distribution is a random variable itself (suppose we have a general distribution, G) thus $\theta_i \sim G_i(.)$. The objective is to get a closed form expression for the following chance constraint: $Prob( \bigcup_{i=1}^{i=m}(A_i > s_i )) \leq \epsilon$, $\forall s_i$ constants   1. If we assume that random variables $A_i$ are independent then the joint constraint becomes :$\sum_{i=1}^{i=m} Prob(A_i \geq s_i) \leq \epsilon $. In that case each $ Prob(A_i \geq s_i)$ can be modeled as a mixture of exponential distribution. If $G_i(.)$ is a discrete probability distribution, in that case we could easily obtain a closed form expression for the constraint . (Is this correct?)    2. If variables $\theta_i$’s are correlated, I am stuck on what could be possible. If I don’t model dependence and assume independence in that case I guess I could get an upper bound on the joint constraint. Also, I was wondering if assuming linear dependence between $\theta_i$’s would make sense or not. Hence I could pose $\theta_i = a_i*\lambda + b_i$ with $a_i$ and $b_i$ different linear coefficients, and $\lambda$ a common random variable. I am again confused whether this would make sense or not. I have never used copula’s but I guess they could be useful in that case (Any comments appreciated) I would appreciate your comments.


title: Need a clear and simple auto-regressive model example

text: This may be hard to find, but I'd like to read a well-explained auto- regressive model example that:   * uses minimal math   * extends the discussion beyond building a model into using that model to forecast specific cases   * uses graphics as well as numerical results to characterize the fit between forecasted and actual values Ideally I would like to understand how to apply the AR model with a simple example. As in - given an array of a car's velocity over time:               [(23,1), (28,2), (19, 3), (25, 4), (30, 5)]      Where first number is the speed and second number is time - how can I apply AR to this data set and forecast the next data point? If you can't provide an example - could you point me to literature where I can find it?


title: What to report and how to report Cox-regression survival analysis in a research paper

text: I am a medical graduate trying to analyze large data and write a paper about it. Can any one explain what should I included in the article about the Cox- regression model. HJP


title: Cox-regression model: question about acceptability, meaning, & reporting

text: Apologize for insufficient info in my previous version of this question. So, I am comparing patient survival in a retrospective study between two different groups, using SPSS. Now the difference came significant by KM, so I tried to adjust multiple covariates for confounding effect in a Cox model. I did a backward Wald test by putting main group in block 1 and co-variates in block 2. The output window gave me this:               Case Processing Summary                      N   Percent     Cases available in analysis Event   6065    10.3%         Censored    24860   42.1%         Total   30925   52.3%     Cases dropped   Cases with missing values   28189   47.7%         Cases with negative time    0   0.0%         Censored cases before the earliest event in a stratum   0   0.0%         Total   28189   47.7%     Total       59114   100.0%     Dependent Variable: patient survival time (days)                           Omnibus Tests of Model Coefficients     -2 Log Likelihood     115511.442 Omnibus Tests of Model Coefficientsi                                          Step  -2 LL  Overall (score)   Change From Prev Step   Change From Prev Block                    Chi-square  df  Sig.    Chi-square  df  Sig.    Chi-square  df  Sig.     1a  113504.047  1998.887    23  .000    2006.433    22  .000    2006.433    22  .000     2b  113504.119  1998.833    21  .000    .072    2   .965    2006.361    20  .000     3c  113504.185  1998.423    19  .000    .066    2   .968    2006.295    18  .000     4d  113504.302  1998.059    17  .000    .118    2   .943    2006.177    16  .000     5e  113506.039  1996.799    15  .000    1.737   2   .420    2004.441    14  .000     6f  113507.146  1986.667    14  .000    1.108   1   .293    2003.333    13  .000     7g  113516.047  1978.462    7   .000    8.900   7   .260    1994.433    6   .000     8h  113518.520  1972.010    6   .000    2.473   1   .116    1991.960    5   .000     aStep 8 TX_MONTHJ   0.011   0.026   0.183   1   0.669   1.011   0.961   1.063         AGE 0.043   0.001   1483.638    1   0   1.044   1.042   1.047         AGE_GROUP   0.701   0.119   34.677  1   0   2.015   1.596   2.544         AGE_DON 0.006   0.001   56.755  1   0   1.006   1.004   1.007         NDIAL2_TRR(1)   0.129   0.041   9.873   1   0.002   1.138   1.05    1.233      I interpret this as saying that no difference was observed, but I'm confused about what to report. Could you also verify if this test looks good with the likelihood ratio so big?


title: Link Anomaly Detection in Temporal Network

text: I came across this paper that uses link anomaly detection to predict trending topics, and I found it incredibly intriguing: The paper is _"Discovering Emerging Topics in Social Streams via Link Anomaly Detection"_. I would love to replicate it on a different data set, but I'm not familiar enough with the methods to know how to use them. Let's say I have a series of snapshots of network of nodes across a period of six months. The nodes have a long-tailed degree distribution, with most having only a few connections, but some having a great many. New nodes appear within this time period. How could I implement sequentially discounted normalized maximum likelihood calculations used in the paper to detect anomalous links that I think might be precursors to a burst? Are there other methods that would be more appropriate? I ask both theoretically and practically. If someone could point me to a way to implement this in python or R, that would be very helpful. Anyone? I know you smart folks out there have some starting thoughts for an answer,


title: Finding ways to bid for items, which has a normal distributed price

text: This is an interview questions, I am not quite sure how to solve. The question is stated as this: Suppose that you want to buy a specific type of car, but you don't know anything about cars (you can't find the underlying quality or the real value). But dealers know about cars and they know exactly how much a car is worth (real value). Now you are given a list of prices that different dealers bought at, but the list doesn't show which dealer bought which car. The price is normal distributed. Assume there are infinite number of dealers and the price dealers bought at is considered as the true value of the car. If you can only bid once for each dealer, and that dealer either accept or decline your offer. Dealer will only accept your offer if and only if your offer is above their buying price. The question is, how to bid or what price to bid, so that you minimize the loss (the price you offer minus the real value). Assume that you only need to buy one car so you can bid as many times as you want until your offer get accepted. (Interviewer mentioned mean and standard deviation, but I still couldn't figure out.)


title: Meaning of "standard deviation units" in a test for statistical significance (basic)

text: For my chemistry class, I have to do some statistics on a set of data. One of the things I have to do is eliminate outliers via the 4D rule. Unfortunately my teacher is also absent, so I can't get this question answered. The 4D rule goes like this on my stat sheet:   1. Calculate mean and standard deviation.   2. Establish a range of acceptable results by adding and subtracting **two standard deviation units** to/from the mean.   3. Calculate new mean and standard deviation after excluding the results that are out of range.   4. Multiply new standard deviation by four and calculate the difference between each of the original results (without exclusions) and the new mean.   5. If the differences calculated is greater than 4 * standard deviation (calculated in step 4), it is valid to exclude the trial. And my data points are: 1.1, 4.5, 9.4, 9.7, 9.9, 10.2, 22.2. (they are rounded short for this post) Mean: 9.6, standard deviation: 5.3 Is it correct to then say that the acceptable range is -1 to 20.2. So I can only exclude the last data point and not the first two? Would the results be difference if I changed the measurements from micrometers to nanometers then?


title: How should tiny p-values be reported? (and why does R put a minimum on 2.22e-16?)

text: For some tests in R, there is a lower limit on the p-value calculations of 2.22e-16. I'm not sure why it's this number, if there is a good reason for it or if it's just arbitrary. A lot of other stats packages just go to .0001, so this is a much higher level of precision. But I haven't seen too many papers reporting p < 2.22e-16 or p = 2.22e-16. Is it a common/best practice to report this computed value or is it more typical to report something else? (like p < 0.000000000000001)


title: Three Parameter Pearson 5 Distribution

text: I need some information about the following - Is "Pearson Type III Distribution" as given in lmomco package (CRAN R) same as Three Parameter Pearson 5 Distribution? If not, how do I estimate the parameters of Three Parameter Pearson 5 Distribution? Regards Katherine


title: Confusion matrices with percentages rather than number of instances?

text: Most of the confusion matrices I've seen contain the number of instances in each cell. Isn't a confusion matrix with the percentage of instances in each cell easier to read? Is this approach wrong or does it go against some unwritten rule with regards to confusion matrices? Such an confusion matrix will look like this, where each of the 10 class labels makes up 10 percent of the dataset and the total is 100 percent. 9.06 percent of the dataset belonged to class 1 and was assigned to class 1. Therefore 90.60 percent of class 1 instances are classified correctly. ![example](http://i.stack.imgur.com/pReOK.png)


title: 'Stationarity requirement' why?

text: I am working on measuring variability of geotechnical data. I see in the literature it is mentioned that, non-stationary data should be first converted into stationary data (for example by trend removal), so that further statistical analysis can be done on the detrended stationary data. Also, it is mentioned that most of the statistical tools assume that the data is stationary, that is why it is important to make the non-stationary data stationary. My question is exactly where is this stationarity assumption is required, i.e., if my data is non-stationary and I still use statistical analysis which result is going to be in error? For example I'm finding the correlation coefficient of the data. If I go on using the non-stationary data will I get wrong value? Why?


title: Various models not improving basic rpart result

text: I have a data set with 10,000 or so samples in it and 100 or so features. I've created a training set and test set and am trying to predict a numeric value. I've used rpart to determine the most important feature by having it generate just two nodes on the training set. I take that and apply it to the test set:               1) root 3288 159847905.80 34.59281022         2) FeatureOne< 0.455 1946  87599096.20 19.00446043 *       3) FeatureOne>=0.455 1342  71090239.25 57.19707899 *      So I can see that only using a single FeatureOne improves the model fairly significantly. My next steps are to use different models to try to improve upon this result. I would think that, in theory, training algorithms that take into account more variables should be able to improve upon the original, basic rpart result of only looking at a single variable. I include FeatureOne as well as many other features that I know are predictive (albeit slightly less than FeatureOne). The problem I'm seeing is that I've tried several modeling algorithms and most of them underperform the basic rpart result that I saw originally. In fact, most are much, much worse than just looking at dataset[dataset$FeatureOne>=0.455,] from the original rpart result. Here's a result using a support vector machine from the e1071 library using the svm command:               1) root 3243 156686211.90 32.84847980         2) predFromSvm< -76.292973 1738  79285639.65 18.49298044 *       3) predFromSvm>=-76.292973 1505  76628786.36 49.42645847 *          rmse = 245.51     cor  = 0.1068      Here's a result using randomForest:               1) root 3243 156686211.90 32.84847980         2) predFromRf< 42.16955833 1631  71891832.78 17.54488657 *       3) predFromRf>=42.16955833 1612  84025916.57 48.33245037 *          rmse = 220.78     cor  = 0.039      Here's a result using glm:               1) root 3243 156686211.90 32.84847980         2) predFromGlm< 32.33149937 1826  79956524.01 15.16580504 *       3) predFromGlm>=32.33149937 1417  75422994.20 55.63504587 *          rmse = 292.29     cor  = 0.0490      Why is it that these algorithms are _worse_ than simply looking at the most important variable as determined by a basic rpart command? I would think these algorithms would be significantly better than this since they are looking at many variables including the highly predictive FeatureOne. Any suggestions on what to try next or ways that I might be using the algorithms incorrectly?


title: Determine an unknown number of real world locations from GPS-based reports

text: I'm working on some software which should determine **real world locations** (f.e. speed cams) from several **GPS-based reports**. An user will be driving when reporting a location, thus the reports a very inaccurate. To solve that problem I have to cluster reports about the same location and calculate an average. My question is about **how to cluster those reports**. I read about Expectation-maximation algorithms and k-means clustering, but as I understood I would need to determine the number of real locations in advance. Are there any other algorithms, which don't need the exact number of real locations, but instead use some edge conditions (f.e. minimal distance) ? A report contains _longitude_ , _latitude_ and _accuracy_ (in meters). There is no name or anything else which could be used to identify duplicates. Another obstacle could be that it will be common, that there is only one report for a real world location. That makes it difficult to distinguish outliers from good data.


title: Can I justify performing a two-way anova where data is normally distributed but has heterogeneous variances?

text: I'm an undergraduate psychology student currently finishing up my dissertation, but a few days ago, while going through my SPSS output, I realized my Levene's test is significant, but my data is normally distributed. Is there any way I can justify that it is ok to still perform a two-way between-subjects ANOVA?


title: What kinds of features are worth taking from objects on images?

text: I managed to label and get set of pixels for every single object of interest in the image. Now I'm looking into clustering and am not sure what features are worth using. I currently have two in mind $x_1 $= Areas of objects $x_2 $= Mean Grayscale Pixel Intensity Now i'm not sure if I should get more features in order for clustering to work better. I was thinking maybe feature 3 be a ratio of feature 1 and 2, or get the histogram of each object, and find a parameter to describe the histogram. I am able to take any of these features of reigionprops as well, but I'm not sure how much useful would half of them be.


title: Exponential Family: Observed vs. Expected Sufficient Statistics

text: My question arises from reading reading Minka's "Estimating a Dirichlet Distribution", which states the following without proof in the context of deriving a maximum-likelihood estimator for a Dirichlet distribution based on observations of random vectors: > As always with the exponential family, when the gradient is zero, the > expected sufficient statistics are equal to the observed sufficient > statistics. I haven't seen maximum likelihood estimation in the exponential family presented in this way, nor have I found any suitable explanations in my search. Can someone offer insight into the relationship between observed and expected sufficient statistics, and perhaps help to understand maximum likelihood estimation as minimizing their difference?


title: Regression to solve system of a system of linear equations?

text: The question I am going to ask probably has a straight-forward answer, I just don't know how to frame my question in a way that would be 'normal.' I have a set of a set of linear equations that I would like to do a regression on. The main equation will always look like the standard matrix linear equation system: $Ax = b$ where A is a 3x3 matrix, x is 3x1 and b is 3x1. However, I can gather data to make 6 equations of this form and A should be the same for each one. $Ax_1 = b_1 \\\ Ax_2 = b_2 \\\ Ax_3 = b_3 \\\ Ax_4 = b_4 \\\ Ax_5 = b_5 \\\ Ax_6 = b_6$ Please note that $x_1$, $x_2$, $b_1$, $b_2$, etc. are **not** single values, but are 3x1 vectors. So for a given matrix A, I can produce 6 sets of x's ad b's. Because I am gathering data from sensors that will have noise, I don't want to just solve $Ax_1 = b_1$ and throw away the data I can get from the other 5 systems. Is there a way to do a linear regression not on a system of linear equations, but on a system of a system of linear equations?