title: What could "directional mean" be in this context?

text: I have some reasonable idea of what directional means are in the context of spatial statistics, but I am stumped by this use of the term in the methods section of this computational biology paper (emphasis mine): > The raw ES values were normalized to account for variable numbers of shRNAs > across different genes by dividing the raw ES by the **directional mean** of > a size-matched null distribution generated by 100,000 random permutations of > a hairpin set of the same size. To sum up the relevant part: one is supposed to compute the directional mean of a distribution obtained from a set of scores (no vectors involved). How would one go about getting a directional mean from a set of values (or associated density function)? **Edit:** one piece of information that _could_ be relevant: the score values whose distribution is being plotted and whose "directional mean" is supposed to be extracted, are obtained through a goodness-of-fit test applied to a subset of values ("the hairpin set") out of the whole ranked list (the microarray data, converted to a list of differential expression values). I also doubt the spatial aspect of the microarray data could play any role here (it's long been eliminated through pre-treatment).


title: Can posterior distribution for a continuous variable be greater than one?

text: I already asked this question here, but I am not sure where would be better to ask it? This might sound a dumb question but I am really confused about it. According to Bayes' rule we do have the following: $$p(\theta|X)=\frac{p(\theta)p(X|\theta)}{\int{p(\theta)p(X|\theta)d\theta}}$$I know that probability density function can be greater than one in general but it seems to me because there exist discrete summations of denominator integral which would is greater than the nominator, therefore posterior probability density function cannot be greater than one in any point. Is this correct?! To explain reasoning a bit in more detail: Suppose we are interested in $p(\theta=\theta_0|X)$ and we know that: $$\int{p(\theta)p(X|\theta)}d\theta\approx\sum\limits_{n}{p(\theta_i)p(X|\theta_i)}$$ but now only consider the summations which include $p(\theta_0)p(\theta_0|X)$. Then the denominator will obviously be larger than the nominator.


title: Readable introduction to measure theory

text: I'm interested in learning more about nonparametric Bayesian (and related) techniques. My background is in computer science and though I have never taken a course on measure theory or probability theory, I have had a limited amount of formal training in probability and statistics. Can any one recommend a readable introduction to these concepts to get me started?


title: Suitable Distribution/Methodology for Estimating Data Points

text: I am trying to model future sales data for products which have very low sales volume. I am a programmer with a smattering in statistics, so I apologise in advance if this qustion is naive! My question is what distribution is most suitable for my sales profile and is it possible to verify the distribution. The exact framing of the problem is that we might have a product that has say 6 units of inventory and we may sell 8 units a year with a "standard deviation" of 5 units (i.e. the sales are lumpy so we calculate a standard deviation, but its not really a normal distribution)...we want to say with a certain probability how many days inventory we have left. For high volume products we can assume the normal distribution and its pretty easy to back out the days inventory left. However for low volume products we can't assume normal distribution (as obviously it is bounded by the fact that we can't have less than 0 sales). I have looked at Poisson, but I am not sure if that is the most suitable. Could someone point me to some resources that would help me identify the right model/technique to use? Thanks, Mike.


title: Weibull Distribution v/s Gamma Distribution

text: What is the difference between the intuition behind Gamma and Weibull distributions? Is there any relationship between the two densities ? Kindly help.


title: External model validation using new data for prediction: How large of a drop in $R^2$ is significant?

text: I need to validate a model using `external model validation' and I have a question relating to deciding when a drop in $R^2$ when compared to $R_{prediction}^2$ is significant. **DISCLOSURE** : This is not a homework question. I am seeking clarification on something that is in the textbook that I am using for an applied regression class (textbook: A second course in statistics: regression analysis). I have model whose goal is to predict $y$ from a set of covariates. The model has already been fit and I would like to validate the model using new data. My textbook states the following: **FROM MY TEXTBOOK** One simple technique is to calculate the percentage of variability in the new data explained by the the model, denoted $R_{prediction}^2$, and compare it to the coefficient of determination $R^2$ for the least squares fit of the final model. Let $y_1, y_2, ..., y_n$ represent the $n$ observations used to build and fit the final regression model and $y_{n+1}, y_{n+2}, ..., y_{n+m}$ represent the $m$ observations in the new data set. Then $$R_{prediction}^2= 1 - \left(\frac{\sum\limits_{i=n+1}^{n+m}\left(y_i - \hat{y_i}\right)^2}{\sum\limits_{i=n+1}^{n+m} \left(y_i - \bar{y} \right)^2}\right)$$ where $\hat{y_i}$ is the predicted value for the $i$th observation using the $\beta$ estimates from the fitted model and $\bar{y}$ is the sample mean of the original data. If $R_{prediction}^2$ compares favorably to $R^2$ from the least squares fit, we will have increased confidence in the usefulness of the model. However if a significant drop in $R^2$ is observed, we should be cautious about using the model for prediction in practice. (Page 316) **MY QUESTION** Unfortunately, this is all the textbook says about using this technique and I have no idea how I would go about determining if a significant drop occurred. For my model I have an $R^2$ of 0.8283 and using the technique described in the text, I have a calculated $R_{prediction}^2$ of 0.7444. Could anyone point me into the direction of other sources that would be able to describe the method of determining whether or not this is a significant drop ---I performed several google searches but the results have not turned up anything useful. Thank you for taking the time to read my question.


title: GA or ANN, and suggestions

text: I am working to approximate, to extrapolate rather, values of a function from a few known values. Here is what I have:   1. Less than 100 known I/O pairs.   2. A monotonously positive correlation.en   3. Something that is not exactly Gaussian, Poisson, or inverse binomial. Right now I can choose between a genetic algorithm and an artificial neural network. If I use a genetic algorithm, then I will use it to approximate the relation with a polynomial or a rational function. If I use a neural network, I will simply find the first Google result for a C++ NN library (unless I am advised otherwise). Here are my questions. Suppose I have access to a long time (1-2 week max) as well great computational might (university clusters).   1. Would I achieve a closer approximation with a GA or an ANN?   2. Which would be faster?   3. For those ANN specialists out there: I am unacquainted with the world of ANN technology. If I am looking for portability, programmability, speed, and adaptability for such extrapolation, which C++ toolkit would you recommend? Thank you.


title: Hidden Markov Model: Average steps needed for the filtering density to reach a certain value

text: I need a hint for solving the following problem. The movement of a certain robot is modeled with a HMM. The robot moves on a circular path and at every time step, it either stays in the same location with $\epsilon$ probability or moves to the next location, with $1-\epsilon$ probability (CCW). Each location of the corridor is labeled from $1$ to $N$ and the robot has a sensor telling whether the current location has a prime label or not. The sensor works correctly with $\delta$ probability. It is known that the robot starts at a non-prime location. I am asked to find the average number of steps for the filtering density $p(x_i|y_{1:i})$ to reach a maximum value (mode) equal or larger than %90, given the sensor correct reading probability $\delta=0.99$. I am stuck with this problem but a possible analytical way of solving it which comes to my mind is to treat the average number of steps $i$ where $p(x_i|y_{1:i})$ attains a mode of $0.9$ or greater as a random variable. Then calculating the expected value of $i$, namely $E[i]$ should give me the answer. The problem is I am not sure whether this can be solved by this way and I am suspicious that there is a much easier way to solve this which I fail to see. I am looking for any hints or comments about whether my way of thinking is correct or not. This is a homework question and therefore I am only looking for hints. Thanks in advance.


title: Optimization of multiple objective functions with constraints

text: I have a number of objective functions like:               y1 = a11* x11 + a12*x11*x11+ a13*x12+.......     y2 = a21* x21 + a22*x21*x21+ a23*x22+.......     :::::     :::::     :::::      These are multiple objective functions. However, the constraints of the objective functions have dependency on each other. Something like,               x11+ x21 + x22 < const1     x12 + x21 > const2          :::::     :::::      What is the way to optimize such a system of equations? I would ideally like to use R to do the same? y1 = 0.32 x1 + 0.21 x1*x1 + 0.49 x2... y2... y3 . . . The equations that i have is a non-linear function. These are non-linear regression equations or non-linear Market Mix Models. The x's are TV spend, Digital Spend etc. I want to 'include' all these models, use some constraints on them and optimize the spends in all the models with respect to constraints like...x1 (TV spend) < 100, TV + Digital spend < 500. I want to be able to say that of an amount of 100, i should spend, 30 on model 1, 20 on model 2( equation 2) etc


title: Model selection for unbalanced data

text: How to do model selection for unbalanced data? how many data points from the whole data set should be selected for model selection? how many for training and testing?


title: Increase the sample size of a Latin Hypercube study?

text: I want to create a climate model ensemble, testing 5 parameters (real, uniformly distributed between two values), using a latin hypercube approach. The problem is that I'm not sure how many replications I want to do. Is it viable to do one latin hypercube of 20 samples, and then another of 10? How would I make sure the 30 samples are relatively evenly distributed? Or would it be possible to do 10-sized LHCs, and do multiples? ie. do 3 or 4 LHC samples of size 10, making sure that each LHC was independent of the others?


title: Are non-square latin hypercubes viable?

text: At https://github.com/OpenMDAO/OpenMDAO-Framework/issues/599 it is stated that non-square Latin Hypercube experimental design is not well defined (I assume that for higher dimensions that means hypercube must have the same length in every dimension). I'm wondering why that is. Would it not be possible to, for example, have a LHC scheme where each row in each dimension had _one or more_ samples, and the sample space was relatively evenly sampled, with no correlation between variables? I guess this is actually making a hypercube out of a hyperrectangle, but is there any reason that wouldn't work? For example, I have a bunch of parameters that I want to change in a physical model. Some have 4 states, some only two. Could I make a hypercube where the Parameters with only two states just have those states duplicated (so that there are two model runs with each of those two states)? how about doing the same with 4-, 3- and 2-state parameters? Could I make an LHC scheme where each dimension has length equal to the smallest common multiple of the states count for each parameter (ie. 12)?


title: How to measure smoothness of a time series in R?

text: Is there a good way to measure smoothness of a time series in R? For example,               -1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0      is much smoother than               -1, 0.8, -0.6, 0.4, -0.2, 0, 0.2, -0.4, 0.6, -0.8, 1.0      although they have same mean and standard deviation. It would be cool if there is a function to give me a smooth score over a time series.


title: Martingale Residuals and Schoenfeld Residuals

text: If the Martingale Residuals indicate that a predictor $X$ violates the proportional hazards assumption but the Schoenfeld Residuals indicate that it does not....which one should you trust?


title: Moore-Penrose generalized determinant

text: Is there a function in R to calculate the generalized determinant of a singular matrix? (similar to the `ginv()` used to compute the generalized inverse)


title: Random Forest: more trees gave worse results

text: I'm new to ML and decided to start learning by having a go at number recognition using the MNIST subset on kaggle http://www.kaggle.com/c/digit- recognizer The images are 28x28 pixel greyscale of the digits from 0 to 9, with ~34000 in the training set and 28000 in the test set. Random Forest using the raw image and 2000 trees gives a score of 0.96829. I did some preprocessing of the images to extract more features and trained / tested on a equal sized subsets of the training data. Compared to normal random forest for subsets around 1000 - 3000 I was getting about a 3-4% improvement (e.g. 94% vs 90% for RF) After training with the entire training set for 200 trees, I scored 0.97557 (a ~0.7% improvement on RF). However, increasing to 2000 trees I scored 0.97457 (0.1% less).   * Does this mean the features are 'saturated'? i.e. any minor performance gain / loss is just random fluctuations?   * Is there anything I can try (apart from more features) to improve the result?   * Is there any way to weight some features more than others so they are more likely to be evaluated?


title: Dealing with class imbalance in multi-label classification

text: I have a set of around 300k text examples. As mentioned in the title, each example has at least one label, and there are only 100 possible unique labels. I've reduced this problem down to binary classification for Vowpal Wabbit by taking advantage of namespaces, e.g. From:               healthy fruit | bananas oranges jack fruit     evil monkey | bipedal organism family guy     ...        To:               1 |healthy bananas oranges jack fruit     1 |fruit bananas oranges jack fruit     0 |evil bananas oranges jack fruit     0 |monkey bananas oranges jack fruit     0 |healthy bipedal organism family guy     0 |fruit bipedal organism family guy 1 |evil bipedal organism family guy     1 |monkey bipedal organism family guy     ...        I'm using the default options provided by VW (which I think is online SGD, with the squared loss function). I'm using the squared loss because it closely resembles the Hamming Loss. After training, when testing on the same training set, I've noticed that all examples are predicted with the '0' label... which is one way of minimizing loss, I guess. At this point, I'm not sure what to do. I was thinking of using cost-sensitive one-against-all classification to try to balance the classes, but reducing multi-label to multi-class is unfeasible since there exists 2^100 label combinations. I'm wondering if anyone else have any suggestions.


title: Predictive model from ordinal panel data

text: Here's a hypothetical problem analogous to one that I am having trouble with: Prostate cancer biopsies are scored according to Gleason pattern (0-10). The scale is ordinal. Over 7 warrants surgery. Patients at risk for prostate cancer were followed with serial biopsies at irregular intervals. Our data consists of a series of biopsies (date and result) for many patients (i.e. 1000 results for 100 patients), along with patient characteristics that may influence risk like sex, birthdate (i.e. age), race etc. What models could be estimated with this data to predict the probability of Gleason > 7 at next biopsy? Or if simpler, what model could predict Gleason score at next biopsy? Intuitively it should incorporate the last biopsy results, time since last biopsy, and the other patient characteristics. I generally use R. Relevant papers, packages, etc. would be appreciated. Thanks in advance.


title: What is the Comparison p used in a sample size calculation for a 1-proportion test in Minitab?

text: In order to calculate a sample size for a 1-proportion test in Minitab, one has to input the power and the comparison proportion. Could someone provide a definition for the comparison proportion please?


title: Relationship between LASSO T and LARS number of steps k

text: We can see on the figure (cf _Least Angle Regression_ p30, Efron, Hastie, Johnstone, Tibshirani - link: Least Angle Regression) that there is a direct relationship between:   * LASSO T absolute norm of $\beta$: $T(\beta) = \sum_j\vert\beta_j\vert$   * and the number of steps $k$ computed by LARS. I am trying to find a mathematical or at least a direct relationship between both LASSO T and LARS k, like $A = B$ or $A = x \Rightarrow B = f(x)$.


title: Should I use missing value using imputation or listwise or pairwise deletion methods?

text: I have 60,000 data and around 45% of them is missing and the missing values are random. Can I simply use listwise or pairwise deletion or do I have to use imputation? If imputation is recommended which imputation is the best one?


title: Analyze and visualize participants response towards particular condition

text: **Method:** I presented 15 participants with audiovisual clips. There were six different clips I presented, and for each clip there was emotional incongruence between visual and auditory information. The combinations were:   * visual negative with auditory positive    * visual negative with auditory neutral    * visual positive with auditory negative   * visual positive with auditory neutral   * visual neutral with auditory negative   * visual neutral with auditory positive I asked participants to make an emotional judgements, whether the clip was happy, angry or neutral. There were five repetitions. **Research question** : Are participants driven by visual or auditory information in their emotional choices, when they watch displays with incongruent emotional information? **Basic analysis:** Because there were no “correct” responses per se for the incongruent stimuli, I calculated the tendency to respond correctly when emotion was presented auditorily or visually. The tendency was estimated by subtracting the proportion of “auditory correct” responses from the proportion of “visual correct” responses for the six incongruent conditions. For example if the incongruent display combined visual negative with auditory positive information, and participants responded that it was a negative interaction, then it was counted as visual response. The estimated indices varied from 1 (participants always responded correctly to the visual information) to −1 (participants always responded correctly to the auditory information). Below is a figure produced in R that shows averaged tendency scores, with standard error bars. The data is available here - columns represent conditions, rows - tendency scores for each participant. ![enter image description here](http://i.stack.imgur.com/9cAsP.png) Now, getting to the point I have two questions: > **Question 1:** What would be the best way to check for significant > differences between those different clips? I know I could simply run a series of 30 paired t test's but it doesn't seem as a pretty solution. I don't think I can use ANOVA here - I have no clue how I could define levels for all those incongruent conditions. Maybe I could organize this data in some different way, but I really not sure what it would be. > **Question 2:** Is there a better way you could use to visualize this? I don't love the above plot - in fact I really don't like it. The oversized legend is very heavy to read, and so are the patterns differentiating each condition. I could use color instead, but I am a bit color blind to be honest, so I would prefer to avoid it. **EDIT** _Question 1_ I decided to just compare the 6 pairs which had the same type of incongreunt information using paired t test. Works kind-of-ok to get the differences. _Question 2_ has been answered well by @AndyW and @gung - I decided to use pure SE bars with vertical orientation of x axis. Done in R using `segplot()` - I need to tweak the details, but it's roughly the idea. ![enter image description here](http://i.stack.imgur.com/9wvAN.png)


title: Integral of survival function

text: If `X` is a nonnegative random variable representing the life of a component having distribution function `F`, and S is the survival function and `T` be a nonnegative random variable representing the time until some specied event. What is the result of this integral? $$ \int_0^x \int_0^t S(u) du dt + \int_x^\infty \int_t^\infty S(u) du dt $$


title: Chances of meeting people who's initials contain a J

text: I am not a statistics expert but was interested in the validity of the reasoning below. Somebody argued that many people had the letter 'J' in their initials. Can someone comment on the validity of the following approach and how it could be improved, notably what probability distribution would a list of firs tnames sorted by popularity follows? * * * According to wikipedia the most common first names for males in the US according to the 1990 census are, in that order: > > 'James','John','Robert','Michael','William','David','Richard','Charles','Joseph', > and 'Thomas' Let's use an approximate probability distribution (any better way?) that gives the probabilities to these names according to their order i (starting at 0) following: > P(name_i) = 0.15 / (i + 2) This is a probability distribution for at least 1000 names (sums to 1), and starts like this: > 0.0749, 0.0499, 0.0374, 0.0299, 0.0249, ... meaning that someone has 7.49% chances of being called James. The sum of the probabilities of having a name starting with J from these list is 14% (P(James)+P(John)+P(Joseph)), lets call it _Pj_. Of course Pj is in fact higher because more names than these 8 in the list of 1000 may start with J. Now let's assume a child is given three first names, the probability that at least of them starts with J given pj is **36%**!! (from (1-(1-pj)^3)). So every time you meet a new person, you have one third chances to see him having a J on his business card. Let's call this probability pIj, probability of at least one initial starting with a J. However, when you say overwhelming majority, let's assume that you mean 8 out of 10. The probability that out of 10 people you meet in a row 8 have a name starting with a J is, using binomial distribution, only **0.5%** I'm afraid, by comb(10,8) _pIj^8_ (1-pIj)^2 . However, if you settle with 5 our of 10, your chances rise to **16%**. * * * FYI the original question and post can be found here.


title: - 2 Log L vs. Deviance

text: I have a question. > Is it better to use deviance or - 2 Log L for the likelihood ratio test? Or > are they basically the same thing? SAS returns the "deviance" sometimes. But I just use the -2 Log L value for the likelihood ratio test (e.g. subtract the -2 Log L for two models to compare them). Is is that okay?


title: Sum of uniform distributions

text: Let’s say $X_1$ ~ $uni(0,1)$ $X_2$ ~ $uni(0,1)$ $X_3$ ~ $uni(0,1)$ And $Y=0.1X_1+0.3X_2+0.6X_3$ What’s the $F(Y)$ (i.e., CDF)?


Does the non-normality matter in using regression for prediction? Hi all, In the Q-Q plot of the residuals after linear regression, the residuals turned out to be highly non-Gaussian. Most of the points (95%) are below a 45º straight line. And those below the straight line are all on the lower side. The shape of the curve looks like $f(x)=x^{1/5}$ for $x$ on $[0, 1]$. (The 5% points in the middle of this curve are above the 45º straight line). It turns out that the dependent variable $y$ has data that are highly non- Gaussian. They are all in $[0, 1]$, but mostly clustered around $1$. So I tried various ways of transforming $y$. The latest one I've found was to do $y_{new} = y^7$. Using this transformation, the $y_{new}$ data become much more symmetrical than before, but they still in no way look Gaussian in a histogram. Being very disappointed, I don't have any more weapons in my bag...Instead, I began to wonder, if my end goal is just to get the $\hat{y}$, i.e. the prediction of the data on a wider data-set, does the non-normality even matter, in terms of accuracy for prediction? Please shed some light for me. Thank you! * * * The $y$ data looks very similar to the data generated using the following script:               mydata=rt(10000, df=5)     mydata=mydata[mydata<0.8]     mydata=(mydata-min(mydata))/(max(mydata)-min(mydata))          hist(mydata, 100)


title: Standardized residuals in R's lm output

text: I have a quick question: if I plot the diagnostic plots to an R regression, a couple of them have "Standardized Residuals" as their y-axis such as in this plot: ![enter image description here](http://i.stack.imgur.com/kZwaM.png) My question is this: **what are the residuals standardized over?** That is, let us assume that in my model, there are 100 predicted values hence 100 residuals.   1. Standardized residual $e_i$ is defined as $(e_i - \bar e)/s_e$(realized residual - mean of all 100 realized residuals)/(standard deviation of all 100 realized residuals)?   2. Since each residual $e_i$ is itself a realized value out of a distribution of possible realizations for this single residual $e_i$, is this residual $e_i$ normalized by its own mean $\bar e_i$ and variance $\text{Var}(e_i)$ (as opposed to the mean and variance from all other values 1 to 100 as described above)? I tried finding documentation clarifying this distinction but could not find any that was beyond doubt. Any help answering this question, and providing me tips on how to find such answers through documentation in the future would be immensely appreciated.


title: Control charts of seasonal data

text: I have been asked to develop control charts for my company - or some other means of quickly identifying when ticket sales are "out of control" and not just a random blip. The data we are looking at are seasonal with a positive trend (sales revenue). It has been suggested that I use an exponentially-weighted moving average (EWMA) chart but my (limited) understanding is that these are only appropriate for use on data with a constant mean. My gut instinct is to apply an ARIMA model and to assess whether observations fall within confidence limits derived from that... but I'm not confident enough to be sure this is the right approach. My question therefore is:   1. Is it appropriate to use control charts to monitor revenue sales   2. What would be the best model to use Please bear in mind that it needs to be administered by somebody with undergraduate maths, and explained to true lay people.


title: Online reference for review of introductory statistics material

text: Can anyone suggest a good **review** reference (necessarily free online PDF) to brush up college-level statistics? I'm looking at something which typically covers a first course in statistics ideally covering:   * Absolute basics (mean, variance etc.)   * Large Samples   * Confidence testing   * Hypothesis testing (normal, t, chi, F)   * Linear Regression/Correlation I don't need a book reference. I need something like a quick review which goes over the important aspects. Anything in free form will do : Presentations, Lecture Notes or otherwise. I just need it to be succinct; depth would be much appreciated.


title: Test of independence vs test of homogeneity

text: I am teaching a basic statistics course and today I will cover the chi-squared test of independence for two categories and the test for homogeneity. These two scenarios are conceptually different, but can use the same test statistic and distribution. In a test of homogeneity, marginal totals for one of the categories are assumed to be part of the design itself -- they represent the number of subjects selected for each experimental group. But since the chi- squared test revolves around conditioning on all marginal totals, there are no mathematical consequences to distinguishing between tests of homogeneity and tests of independence with categorical data -- at least none when this test is used. My question is the following: is there any school of statistical thought or statistical approach that would yield different analyses, depending on whether we are testing for independence (where all marginals are random variables) or a test of homogeneity (where one set of marginals are set by the design)? In the continuous case, say where we observe $(X,Y)$ on the same subject, and test for independence, or observe $(X_1, X_2)$ in different populations and test if they come from the same distribution, the method is different (correlation analysis vs t-test). What if the categorical data came from discretized continuous variables? Should the tests of independence and homogeneity be indistinguishable?


title: Testing for variable independence with count data and assumption violations

text: I used a survey to collect information about injuries suffered while practicing and competing in gymnastics. I have count data on the injuries that looks something like this:                          Floor   Vault   Rings ...     Head         33      45      7     Neck         43      19      17     Shoulder     21      39      25     ...         ...     ...     ...     Toes         16      9       3      I wanted to see if there was a relationship between the body part injured and the equipment being used. I thought that using a chi square test of independence would work but I've hit a snag, I think. When the participants filled out the survey, they were able to check multiple responses. For example, if someone was injured while performing the vault, they could have injured both their head and neck and both would show up in the injury count. Thus one injury contributed to two cells. I've considered running separate analyses for each body part but I don't think that exactly answers my question. Is there another way to tackle this problem or a different statistical test that might work?


title: Choosing significance level alpha and max_lag for Granger causality

text: I am really new to Granger Causality and am just rabidly experimenting with this MATLAB function:               granger_cause(x,y,alpha,max_lag)       http://www.mathworks.com/matlabcentral/fileexchange/25467-granger-causality- test I have two series x and y which I need to check for causality. I need to check if y at a lag of 180 samples has a causal relationship with x or not. Based on my experimentation, I noticed that if I vary the values of alpha and the lag, then I am able to find that almost anything seems to correlate with x (such as a random set of data). Can someone tell me just how should I should vary the alpha and lag values? I found that if I put alpha as 0.8 then I get a good answer for x and y but not for the random data; however considering that other resources tell me to have alpha at 0.01 or so, I'm pretty suspicious about this. Can someone please help me out here?


title: Prediction errors in orthogonal regression

text: I am reading up on orthogonal regression for possible implementation in one of my datasets. My question relates to computation of prediction errors as "mean prediction errors -MPE" or "mean absolute prediction errors-MAPE" in orthogonal regression where my dependent and independent variables (both measured in the same units) may not have the same errors as they may come from different sources. In normal least squares where the dependent variable is regressed on the independent variable, the MPE or MAPE can be taken from mean of the predicted - observed. How does one compute the MPE or MAPE in the scenario when both the variables have error in them. This may be simple but was not obvious to me. Hope the question is clear enough. VJ


title: What is an alternative to the book "Probability and Statistics" and covers the same material?

text: Can someone recommend a textbook(s) that covers the same topics as Probability and Statistics by by Michael J. Evans and Jeffrey S. Rosenthal? I find this particular textbook to be very poorly written and terse.


title: Confusion related to expectation notation

text: I have a confusion regarding notation of the expectation. Consider the following: $$E_x(x^2) = \sum_x x^2p(x)$$ $$E_x(x|y) = \sum_x xp(x|y)$$ So the same notation here is giving different meaning right isn't it? Up to now I was assuming $E_x$ taking $p(x)$ and calculating the expectation. However after looking at $E_x(x|y) = \sum_x xp(x|y)$, I am a bit confused


title: Detect current activity of a time series

text: Let's say I have a time series that indicates the number of visitors to a website, so each value is non-negative. For example: [100, 120, 40, 200, 1000, 1400, 1300, 1500] In this case it's clear that the website is busier over the last days. Similarly, [2000, 1000, 100, 100, 50, 70, 20, 30, 20, 30] would indicate that the traffic on the website is decreasing. What would be a good way to calculate a value in the range [0, 1] to indicate this? Moreover, assuming these two series belong to two different sites, how would I be able to compare one with another? (I assume some kind of normalization should be made over the results, is that correct?) Thanks in advance.


title: fitting a model of each subject and perform a t-test over the estimated parameters

text: In my field, psychophysics/experimental psychology, sometimes a lot of measures are collected for each subject. Then, a model is fitted for each subject and with the estimated parameters sometimes a global test like a t-test or ANOVA is performed to reach conclusions. It seems that this approach sometimes produces wrong conclusions and I wanted to know references that already made this point. Below a toy example: Let's suppose that I want to know whether the color of an object (red or green) influences its visibility. To measure visibility I presented a tiny object in a computer screen and ask the subject to report whether can or cannot see the object. After she responds, I display another tiny object and so on. I use four different sizes and two colors and presented them randomly. Each size and color is presented 50 times. This is a possible outcome of the experiment: ![Results toy example](http://i.stack.imgur.com/3ou12.jpg) Each panel shows the data for one subject. The y-axis is the frequency of reports 'yes, I saw it'. The curves are cumulative normal distributions fitted using generalized linear models and the vertival lines are the estimated means of the underlying distributions. The confidence intervals of the fits indicate that the color influences the visibility for each subject. For subjects 1,2 and 3 green is easier to perceive. For subjects 4,5 and 6 red is easier. Given that the color does not affect equally each subject, it is possible that if I use the mean of the fitted distributions to perform a t-test it turns out not significant. So sometimes in my field, based on the results of the t-test, it is concluded that color does not affect visibility. The individual analysis for each subject is not performed and discussed. So I am searching for references (articles/books) that point this problem out. Below the code of the toy example in R:                library('plyr')      library('ggplot2')           set.seed(100)      n<-50      size<-1:4      prob<-c(0.1,.3,.5,.7)           variables1<-expand.grid(subject=1:3,color=c('red','green'))      df1<-ddply(variables1,.(subject,color),function(d){          if (d$color=='red') nyes<-rbinom(4,n,prob)              if (d$color=='green') nyes<-rbinom(4,n,prob+.2)          data.frame(size,nyes,nno=n-nyes,p=nyes/n)       })           variables2<-expand.grid(subject=4:6,color=c('red','green'))      df2<-ddply(variables2,.(subject,color),function(d){         if (d$color=='green') nyes<-rbinom(4,n,prob)             if (d$color=='red') nyes<-rbinom(4,n,prob+.2)         data.frame(size,nyes,nno=n-nyes,p=nyes/n)       })           dat<-rbind(df1,df2)           curves<-ddply(dat,.(subject,color),function(d){        model<-glm(cbind(nyes,nno)~size,binomial(probit),d)        x<-seq(0,5,.1)        y<-predict(model,data.frame(size=x),type='response',se.fit=T)        data.frame(x,y=y$fit,ymin=y$fit-y$se.fit,ymax=y$fit+y$se.fit)      })           means<-ddply(dat,.(subject,color),function(d){        model<-glm(cbind(nyes,nno)~size,binomial(probit),d)        coe<-coef(model)        m<- -coe[[1]]/coe[[2]]       data.frame(m)      })           q<-ggplot()+        facet_wrap(~subject)+        geom_point(data=dat,aes(x=size,y=p,color=color))+        geom_line(data=curves,aes(x=x,y=y,color=color))+        geom_ribbon(data=curves,aes(x=x,ymin=ymin,ymax=ymax,fill=color),alpha=.3)+        geom_linerange(data=means,aes(x=m,color=color),ymin=0,ymax=.5)      q           t.test(m~color,data=means,paired=T)


title: how to derive Kolmogorov distribution $K_\alpha$ for big n?

text: The Kolmogorov–Smirnov test rejects the null hypothesis at level $\alpha$ if $$\sqrt{n}D_n>K_\alpha$$ where $$D_n=\sup_x |{F}_n(x)-F(x)| $$, and $$\operatorname{Pr}(K\leq K_\alpha)=1-\alpha$$ , where $K$ follows Kolmogorov distribution. There are tables for $K_\alpha$, but this website gives formula for $K_\alpha$ for big $n$. For example, it says, if $n$ is over 50, $$K_{\alpha=0.01} = 1.63 $$ $$K_{\alpha=0.05} = 1.63$$ $$K_{\alpha=0.01} = 1.36$$ $$K_{\alpha=0.1} = 1.22 $$ $$K_{\alpha=0.15} = 1.14 $$ $$K_{\alpha=0.2} = 1.07 $$ How could such number be reached?


title: What is the difference between the "coef" and "(exp)coef" output of coxph in R?

text: I have been trying to discern what exactly the "coef" and "(exp)coef" output of coxph signify. It seems that the "(exp)coef" are comparisons of the first variable in the model according to the group assigned in the command. How does the coxph function arrive at the values for "coef" and "(exp)coef"? Additionally, how does coxph determine these values when there is censoring involved?


title: What is the "Value" output of Survreg in R?

text: I am exploring the use of the `survreg` function in R to analyze my current experiment. Does anyone know what the "Value" column in the output of the function stands for?


title: Determining odd time series

text: I have a number of time series which are derived of same underlying data. However, the data in each comes from a different source, so they may be slightly lagged or differently enriched but essentially representing the same underlying data. Is it possible to find out if there is one of these series that is odd and severely delayed or out of line? As an example, imagine you subscribing to market data for the same stock from 5 different providers and trying to find out which is running slow in a dynamic way. However, stock prices are very deterministic so that makes the problem slightly simplified. In my case I am working on data that doesn't have a single source to check.


title: How do you check the linearity of a multiple regression

text: With Univariate regression it's very easy to see if the correlation is linear or otherwise, but in a multiple regression things can look very complicated, polynomial or other correlations can be hidden by the interaction of of the variables. I understand that you can use power terms in the variables but the relationship to the dependent variable should still be linear. How do I know if multiple regression is appropriate?


title: Tweaking axis properties in a timeseries graph

text: The following code:               require(zoo)     data <- read.csv(file="summary.csv",sep=",",head=TRUE)     cum  = zoo(data$dcomp, as.Date(data$date))     data = zoo(data$compressed, as.Date(data$date))     data <- aggregate(data, identity, tail, 1)     cum  <- aggregate(cum, identity, sum, 1)     days = seq(start(data), end(data), "day")     data2 = na.locf(merge(data, zoo(,days)))          pdf(file='timeseries.pdf',width=9,height=5)     par(bty = 'n')     plot(data2,xlab='',ylab='entropy (bytes)',axes=FALSE)     axis(side = 2, at=c(991, 20000, 40000, 53048))     axis(side = 2, at=c(10000, 30000),labels=FALSE)     lines(cum,type="h",lwd=0.3,col=rgb(0.64,0.08,0.00))     axis.Date(side = 1, days, at=c("2007-07-25", "2008-01-01", "2009-01-01", "2010-01-01", "2010-06-21"))      Yields the following Graph: ![alt text](http://i.stack.imgur.com/k3Yei.png) The questions (all related to the bottom axis) are:   1. How can I make it display '2007/07' and '2010/06' in the beginning and end labels?   2. How can I make it to automatically add minor ticks to the months (without labels)?


title: Analysing regression results

text: This question is a bit embarrassing, but I do not have much experience of statistical analysis. The problem I have is that the analysis I do of the results in my paper is so meager. My results from a fixed effects regression basically consist in the effect of one or two central explanatory variables. I can comment on their sign, significance and magnitude. I can perhaps discuss whether there is a problem of omitted variable bias in the model. Is there something more I can do? As of now, the analysis won’t even be a page long. My knowledge is limited and I would appreciate all kind of advice or input here. The research question is very narrow: are there economies of scale in health care at the state/county level? I have a cost measure as dependent variable and the number of medical care events and physician visits as the explanatory variables capturing the size of the health care. Then I try to control for geography, demographics and socioeconomic status. A limited specification could look like this: $$ \text{Cost}_{it} = \beta_0 + \beta_1 \text{number_of_visits}_{it} + \\\beta_2 \text{number_of_events}_{it} + \beta_3 \text{percentage_age65}_{it} + \\\beta_4 \text{area}_i + \beta_5 \text{taxable_income}_{it} + \\\beta_6 \text{percentage_private_healthcare}_{it} + \epsilon_{it} $$


title: Connection between feature selection and hypothesis testing

text: I have a dataset for cell-phone accounts and I am trying to predict whether or not an account will cancel given some input features. One such feature is the number of devices an account owns. I am trying to figure out if this is a predictive feature. I was thinking that one way to do this is to calculate the mean number of devices for all positive examples (accounts that do not cancel) and negative example (accounts that do cancel). Then I would do a hypothesis test for means to determine if they are significantly different. So I'm just wondering if this is a valid approach for feature selection.


title: Overfitting in Genetic Programming

text: I've recently started experimenting with Genetic Programming as an optimization tool. I'm still a little confused as to how to reduce overfitting in this framework. A couple of techniques I've read about   1. Limiting the number of generations   2. Limiting the tree (genome) size   3. Adding a complexity penalty to fitness function   4. Incorporate cross validation into fitness function Of the methods above, 4 is the most intriguing/confusing to me. I've read this post on SO How to avoid overfitting when using crossvalidation within Genetic Algorithms, but it didn't give me too much insight on the actual implementation. Let's say I'm trying to implement cross validation in my fitness function. I can "train" my GP on some time period (let's say 2000-2010) and evaluate it's fitness on another time period (let's say 2010-2014). However, won't this just introduce overfitting on the time period from 2010-2014. The GP will just gravitate towards finding good programs on the testing time period. Alternatively, you could say that you want a program that performs well on 2000-2010 (performance P1) and 2010-2014 (performance P2). Say your fitness function requires that the performance of the GP on those two time periods has to be high and similar (i.e. Fitness = (P1 + P2) * 1/(P1 + P2)), but then isn't that the same as just testing the fitness on the whole period 2000-2014? I guess I'm asking is it even possible to implement a form of cross validation into your fitness function. GP isn't a parameterized algorithm like kNN, SVM, or Simple Linear Regression, so I'm confused as to how cross validation even works in this setting. I can always just evaluate the output of the naive GP on test data, but in my experiments I haven't found good performance, because of GP's tendency to overfit.


title: A Graphical Model for Fellegi-Sunter Record Linkage

text: I am trying to understand the Fellegi-Sunter Probability Model for Record Linkage problem. I am following the thesis at: http://www.inf.ed.ac.uk/publications/thesis/online/IM080663.pdf in order to learn about the model. According to my current understanding this model handles the linkage problem as a classification problem where we have record pairs $(a,b)$, their distance measures as the distance vector $\gamma(a,b)=(\gamma_1(a,b),...,\gamma_K(a,b))$ and $M(a,b)$ as the binary variable indicating whether the pair $(a,b)$ is a correct match or not. If I understand all the mechanism correctly, it should be something similar to a Naive Bayes Classification problem. I tried to model it using a Bayesian Network, where we have $N$ records, $q$ is a small prior probability for each $M_i$ to indicate a correct match and $\theta_1,...,\theta_K$ are the parameter vectors with uniform priors, for each distance $\gamma_1,...\gamma_K$. The model is like in the following: ![enter image description here](http://i.stack.imgur.com/R20PQ.png) Since the data does not have match-unmatch labels $M_i$ we have to run the EM algorithm in order to make a maximum likelihood estimation for the parameters $\theta_1,...,\theta_K$. My question is, is this model correct for this record-linkage problem? I could not find a proper graphical model for this problem so I turned to here for answers. Thanks in advance.


title: Interpretation of strucchange::breakpoints in R

text: OK. I've just read that a breakpoints object with "NA" means no significant structural breakpoint was found, which is fair enough. But then if you call summary() on that breakpoint object, it will report a breakpoint. Who to believe?! Also, my independent variable is not time. I'm searching for a change in the relationship with an independent variable, so the "breakdate" terminology is inappropriate. I can deal with that. I simply extract the nth value from the sorted vector of independent variable that the breakpoints were built on. But the Fstat at that value (as far as I can tell) is not at the maximum. I thought the best breakpoint would be located at the value with the maximum F stat. Is that not the case? Thanks, Eric EDIT: Here's how I got confused. I'm looking for breakpoints in a model Y~X, where X has a range of roughly 0 to 150. The summary() says there's a breakpoint at observation 511 (out of 910 data), corresponding to breakdate 0.56. But the value of X[511] is 120 or so. So, I assumed (oops) that I should use the observation number and not the breakdate. And that was when that breakpoint didn't line up with the maximum F stat from Fstats. (But nor would a breakpoint at 0.56). A similar investigation of a model where the X predictor has a range in 0 to 1 appears to yield more comprehensible outputs. The reported breakdate (interpreted as a value of the predictor) appears to line up with the FStat maximum. Is this a dangerous coincidence, or does breakpoints() require something of its predictor variable? And sorry, I don't know how I can attach a dataset to reproduce this. :-(


title: Hypothesis test for difference between domain means with weighted data

text: I have a large dataset that is weighted due to:   1. probability proportional to size sampling, &   2. disproportionate stratification.  I would like to compare the means for two different domains within the same sample, using the weighted means and variance estimates. I believe I could proceed using the Delta approximation of variance of a function of random variables: $$ {\rm V}(\mu_1 - \mu_2) = {\rm V}(\mu_1) + {\rm V}(\mu_2) - 2{\rm Cov}(\mu_1, \mu_2) $$ but I cannot find any precedent, nor a different solution, to this problem in the literature. Nor could I find any other questions on this forum that address this specific problem.


title: How can I represent R squared in matrix form?

text: This question is a follow-up to a prior question. * * * Basically, I wanted to study under what conditions when we regress the residuals to $x_1$, we will get $\small R^2$ of 20%. As a first step to attack this problem, my question is, how do I express $\small R^2$ in matrix form? Then I will try to express "$\small R^2$ of regressing residuals to $x_1$" using matrix form. Also, how can I add regression weights into the expression?


title: Data Augmentation Examples

text: I am looking for applied references to data augmentation (preferably with some written code). Either online references are books would be great. I found this book online: http://www.amazon.com/Bayesian-Missing-Data-Problems- Biostatistics/dp/142007749X/ref=sr_1_1?ie=UTF8&s=books&qid=1291905761&sr=1-1 But with no reviews I am hesitant on purchasing it. Thanks! Edit: I have two variables X and Y. Let's say X follows a mixture of normals and there is a logistic relationship between X and Y. There is measurement error when observing X. We observe 100 X Y pairs and need to estimate the function between the two. In a book on measurement error (John P. Buonaccorsi) the author recommends data augmentation (I believe the introduced variables are the true X means) for estimation. However no details are given. I am looking for simple examples (R code but doesn't really matter) to get started.


title: How do I see if the distributions of two variables are overlapping in SPSS?

text: This seems to be so straightforward that I can't believe I haven't found the answer in the SPSS help or online - I must be missing something. I have two scale variables, let's call them $X$ and $Y$. They are likert scale variables measuring agreement with attitude statements for two items that are thought to measure the same construct. It's two separate questions from the same instrument, given to the same sample (which was drawn from the same population). (*) These two variables have some distributions that have some overlap. My question is, how do we test the hypothesis that the two distributions are or are not different (to a p<0.05 level)? (Which, I think, could imply that the items are indeed measuring the same construct.) This seems to call for a t-test comparing whether the means $\mu_X$ and $\mu_Y$ of the two distributions are equal, testing the hypothesis: $\mu_X=\mu_Y.$ But I've looked through the documentation of the "Compare Means," "Independent samples T-Test," and "Paired samples T-test" in SPSS, both in the GUI and the command syntax, and I can't find how to do this type of test. I must be missing something? (*) I know there is some controversy about whether to treat Likert Scale variables as ordinal or scale variables, but for this preliminary analysis I'm comfortable treating them as scale.


title: Good GUI for R suitable for a beginner wanting to learn  programming in R?

text: Is there any GUI for R that makes it easier for a beginner to start learning and programming in that language?


title: Find out pseudo R square value for a Logistic Regression analysis

text: My name is Tuhin. I came up with a couple of questions when I was doing an analysis in R. I did a logistic regression analysis in R and tried to check how good the model fits the data. But, I got stuck as I could not get the pseudo R square value for the model which could give me some idea about the variation explained by the model. Could you please guide me on how to achieve this value (pseudo R square for Logistic regression analysis). It would also be helpful if you could show me a way to get the Hosmer Lemeshow statistic for the model as well. I found out a user defined function to do it, but if there is a quicker way possible, it would be really helpful. I would be very grateful if you can provide me the answers to my queries. Eagerly waiting for your response. Regards


title: Moving return of exponential moving average -- choice of alpha

text: I have a time series with an exponential moving average and I want to calculate a moving return of the EMA over the last m periods (something like a smoothed moving return). Let's say: Y(t) is the value of the time series at time period t S(t) is the value of an EMA of Y at time period t Now R(t) is the return of the EMA over the last m time periods: R(t) = S(t) / S(t-m) - 1 My question is: how many time periods should the EMA calculation use for a given m? Precisely, if the EMA is calculated using S(t) = alpha * Y(t) + (1-alpha) * S(t-1) and alpha is set by 2/(N+1), then how should N depend on m? I'm assuming that N should be sufficiently less than m to prevent 'overlap' of Y values that are used in the calculation of S(t) and S(t-m). Any theories or best practices about this?


title: Difference between naive Bayes & multinomial naive Bayes

text: I've dealt with **Naive Bayes** classifier before. I've been reading about **Multinomial Naive Bayes** lately. Also **Posterior Probability = (Prior * Likelihood)/(Evidence)**. The only prime difference (while programming these classifiers) I found between Naive Bayes & Multinomial Naive Bayes is that **Multinomial Naive Bayes** calculates likelihood to be **count of an word/token** (random variable) and **Naive Bayes** calculates likelihood to be following: ![enter image description here](http://i.stack.imgur.com/4Q5L1.png) **Correct me if I'm wrong!**


title: Variable frequency redistribution

text: I know that there is a way to "redistribute the frequencies of a variable" as stated here: Slide number 14 and 15 about redistribution and here: Dorian Pyle book chapter 7 section 2 paragraph 3 (7.2.3): > The easiest way to adjust distribution density is simply to displace the > high-density points into the low-density areas until all points are at the > mean density for the variable. Such a process ends up with a rectangular > distribution. Illustrating this passage are figures like these: ![Figures 7.7 and 7.8](http://i.stack.imgur.com/4cQpc.png) In both the references the problem is clear and solved but does not provide sufficient information for writing an algorithm and in particular does not say how to calculate the displacement (the way to move the frequencies on the left or right of the center of the distribution). Please can you help me to understand it ??


title: What is the distance between a finite Gaussian mixture and a Gaussian?

text: Suppose I have a mixture of finitely many Gaussians with known weights, means, and standard deviations. The means are not equal. The mean and standard deviation of the mixture can be calculated, of course, since the moments are weighted averages of the moments of the components. The mixture is not a normal distribution, but how far from normal is it? ![Mixture of Gaussians separated by 2 standard deviations vs. Gaussian with same mean and variance](http://i.stack.imgur.com/kJ841.png) The above image shows the probability densities for a Gaussian mixture with component means separated by $2$ standard deviations (of the components) and a single Gaussian with the same mean and variance. ![Mixture of Gaussians separated by 1 standard deviation vs. Gaussian with the same mean and variance](http://i.stack.imgur.com/VGtsR.png) Here the means are separated by $1$ standard deviation and it is harder to separate the mixture from the Gaussian by eye. * * * **Motivation:** I disagree with some lazy people about some actual distributions they haven't measured which they assume are close to normal because that would be nice. I'm lazy, too. I don't want to measure the distributions, either. I want to be able to say their assumptions are inconsistent, because they are saying that a finite mixture of Gaussians with different means is a Gaussian which isn't right. I don't just want to say that the asymptotic shape of the tail is wrong because these are just approximations which are only supposed to be reasonably accurate within a few standard deviations of the mean. I'd like to be able to say that if the components are well-approximated by normal distributions then the mixture isn't, and I'd like to be able to quantify this. * * * I don't know the right distance from normality to use: supremum of differences between CDFs, $L^1$ distance, earth-mover's distance, KL divergence, etc. I'd be happy to get bounds in terms of any of these, or other measures. I would be happy to know the distance to the Gaussian with the same mean and standard deviation as the mixture, or the minimum distance with any Gaussian. If it helps, you may restrict to the case that the mixture is of $2$ Gaussians so that the smaller weight is greater than $1/4$.


title: Where did the term "learn a model" come from

text: Often I have heard the data miners here use this term. As a statistician who has worked on classification problems I am familiar with the term "train a classifier" and I assume "learn a model" means the same thing. I don't mind the term "train a classifier". That seems to portray the idea of fitting a model as the training data is used to get good or "improved" estimates of the model parameters. But the would learn means to gain knowledge. In plain English "learn a model" would mean to know what it is. But in fact we never "know" the model. Models approximate reality but no model is correct. It is like Box said "No model is correct but some are useful." I would be interested to hear the data miners response. How did the term originate? If you use it, why do you like it?


title: Convergence of distribution

text: This is from _Probability and Measure_ by Billingsley, 3rd Edition. **27.21 (p. 370)** Let $X_1, X_2,...$ be independent and identically distributed, and suppose that the distribution common to the $X_n$ is supported by $[0,2\pi]$ and is not a lattice distribution. Let $S_n=X_1+\cdots+X_n$, where the sum is reduced modulo $2\pi$. Show that $S_n \Rightarrow U$, where $U$ is uniformly distributed over $[0,2\pi]$. Can someone provide some hints? Thanks! P.S. This problem refers back to other two problems. Namely: **26.1 (p. 353)** A random variable has a _lattice distribution_ if for some $a$ and $b$, $b>0$, the lattice $\{a+nb:n=0,\pm 1,\dots\}$ supports the distribution of $X$. Let $X$ have characteristic function $\varphi$. **(a)** Show that a necessary condition for $X$ to have a lattice distribution is that $|\varphi(t)|=1$ for some $t\neq 0$. **(b)** Show that the condition is sufficient as well. **(c)** Suppose that $|\varphi(t)|=|\varphi(t')|=1$ for incommensurable $t$ and $t'$ ($t\neq 0$, $t'\neq 0$, $t/t'$ irrational). Show that $P\{X=c\}=1$ for some constant $c$. **26.29 (p. 356)** **(a)** Suppose $X'$ and $X''$ are independent random varibles with values in $[0,2\pi]$, and let $X$ be $X'+X''$ reduced module $2\pi$. Show that the corresponding Fourrier coefficients satisty $c_m=c_m' c_m''$. **(b)** Show that if one or the other of $X'$ and $X''$ is uniformly distributed, so is $X$.


title: In a one sample t-test, what happens if in the variance estimator the sample mean is replaced by $\mu_0$?

text: Assume a one-sample t-test, where the null hypothesis is $\mu=\mu_0$. The statistic is then $t=\frac{\overline{x}-\mu_0}{s/\sqrt{n}}$ using the sample standard deviation $s$. In estimating $s$, one compares the observations to the sample mean $\overline{x}$: $s=\sqrt{\frac{1}{n-1}\sum_{i=1}^n (x_i-\overline{x})^2}$. However, if we assume a given $\mu_0$ to be true, one could also estimate the standard deviation $s^*$ using $\mu_0$ instead of the sample mean $\overline{x}$: $s^*=\sqrt{\frac{1}{n-1}\sum_{i=1}^n (x_i-\mu_0)^2}$. To me, this approach looks more natural since we consequently use the null hypothesis also for estimating the SD. Does anyone know whether the resulting statistic is used in a test or know, why not?


title: Selecting samples for calibation in regression model

text: I am developing a multiple regression model. Recently I came across a statement that for the calibration set it's "better" to select samples which are uniformly distributed in a response variable space -Y as opposed to "normally" distributed samples. I am trying to understand this. If for the calibration set we select samples which have roughly uniform distribution will it help to achieve better prediction for low and high values of Y? Thank you


title: Confidence interval of quantile / percentile of the normal distribution

text: What is the formula (if it exists) for the sample variance / confidence interval of a quantile / percentile of the normal distribution? For example, the 5th percentile for a standard normal population distribution is -1.64, but what is the 95% one-sided confidence interval if I have a sample of n=1,000? I.e. on average the 5th percentile of a standard normal sample will be -1.64 and 95% of the time the sample 5th percentile of a sample of n=1,000 will be below -1.54 (approximate from simulations). There seems to be some kind of formula as it's an option in SAS (`CIQUANTNORMAL`) and I think it may be related to the Probit but haven't found an explanation.


title: What does it mean to integrate over the posterior?

text: I have been reading a book that cites an example where a uniform distribution is the initial prior, and then a person scores 9/10 on a test. Then the resulting posterior becomes the prior distribution. The book provides the following explanation for what happens, but despite consulting the Wikipedia page on integration I am struggling to understand exactly what is going on, and why this process is necessary. ![enter image description here](http://i.stack.imgur.com/STt6A.png)


title: Data cleaning for large sample data set in multiple linear regression

text: I have 70,000 observations for my depentant variable. I have 12 independant variables. After removing zero value and error and missing value form my data set, my data reduced to 4000. Can I still do the multiple linear regression with this data set? I think 4000 data is more than enough for 12 indepndant variables, but I am not sure wheather removing almost 90% of observertions will harm my regression or not?


title: Why is the default cost function choice of a neuron quadratic loss?

text: I'm studying neural networks, and I'm trying to decide why the default choice of cost function for a single neuron seems to be quadratic loss: $$\sum_i(y_i- f_i)^2,$$ instead of: $$-\prod_ip_i^{y_i}(1-p_i)^{1-y_i},$$ as per logistic regression. Where both $f_i$ and $p_i$ are the sigmoid (activation) function. I understand that the neuron does not classify per say but instead modulates/dampens the output based on its activation function, so that when many neurons are connected together to form the network, classification only needs to be done on the output node using some cut-off value. Nevertheless, if classification is our goal, it seems like the right choice of cost function should be that which relates the strength of the firing directly to the probability that it would fire if it were only firing at total strength or not at all. And doing this in such a way that if we ran the neuron many times using this probability as the probability that it would fire, then the expected value of the neuron firing would equal the strength with which it fires every time upon minimizing our cost function. (sorry I may have worded that poorly) Is there a reason then that the logistic regression cost function is not the default choice?


title: Error Bars in a Monte Carlo Coin Experiment

text: An excel program runs a Monte Carlo coin experiment. `10,000` coins are tossed `1000` times, with the number of times `N` heads appearing recorded in a table. Call the number of times `N` heads appear during a run of the experiment `Q`. The average value of `Q` as well as its standard deviation is calculated at the experiment's end. My question is how do I calculate the error bars for such an experiment? I am not exactly sure how to do this, and I've been trying for the last few hours to find a solution, but all attempts have been fruitless. Any help would be greatly appreciated.


title: Graphical summaries of relationships between univariate distributions

text: I'd like to review published papers or book chapters (so I could formally refer to them) that graphically illustrate the parametric relationships between univariate distribution families. The papers by Lawrence M. Leemis are well-known, see e.g. here for 1986 and 2008 versions. This paper presents two more diagrams and mentions two more sources for previously published graphical summaries: Nakagawa, T. and Yoda, H. (1977) Relationships among distributions. IEEE Transactions on Reliability, 26(5), 352–353. Kotz, S. and Jogn Rene van Dorp (2004) Beyond Beta—Other Continuous families of Distributions with Bounded Support and Applications, World Scientiﬁc, NJ, p. 251. Here is another link to a paper with presumably original contribution. Are there any useful omissions in the above list? Is there a paper that reviews various graphical summaries of the known relationships between univariate distributions?


title: Why not report the mean of a bootstrap distribution?

text: When one bootstraps a parameter to get the standard error we get a distribution of the parameter. Why don't we use the mean of that distribution as a result or estimate for the parameter we are trying to get? Shouldn't the distribution approximate the real one? Therefore we would get a good estimate of the "real" value? Yet we report the original parameter we got from our sample. Why is that? Thanks


title: How is power for t test and Mann Whitney comparable when null hypotheses differ

text: In a previous question, I asked about comparing the power of a t test to a Mann Whitney test under different situations. One of the answers pointed out that the worst that the Mann-Whitney can ever perform relative to the t-test is that it would require 1/0.864x as much data to give the same power as the t test so long as the data sets being compared were from the same distribution. I guess I am confused about how power can be compared between these two tests given that they test different null hypotheses. For example, if I estimate the power for a t-test at 0.9 using effect size equal to a 20% difference between means, then that makes sense for a t test. But a Mann Whitney test does not test for differences between means. If the distributions were the same it would test for differences between medians. Am I right in thinking that the Mann Whitney test would require 1/0.864x as much data to detect a 20% difference in medians with a power of 0.9?


title: Time-series machine learning methods and R packages

text: I am trying to determine how to use machine learning models such as for eg., random Forest with (non-financial) time-series data. Using an example, suppose we wanted to find based on monthly scores on subjects for each student how well he/she will do in a month-end exam that occurs every month.               Month 1 Data          Name    State  EngScore  MathScore  HistScore  ...  MonthEndExamScore     John    NY     80        90         75              180     Jack    TX     78        65         90              170     John    CA     82        93         79              185      The same data is collected for the students in Month 2, 3 ... n. The task is to predict the MonthEndExam score of the current month using the student's historical data on performance during previous months. For the current month, the scores on all the individual subjects are known by mid-month, whereas the MonthEndExam is not known until the end of the month and it is what we would like to predict. I understand that one could use statistical methods such as ARIMA, etc, but I was wondering if there were similar methods in ML that can be applied here, ** in particular ** using packages in R (such as randomForest, caret, party, etc). Thanks in advance.


title: How to evaluate the quality of probability estimation?

text: I have a classifier that gives its decisions as probability estimations: for each datum it returns a set of probabilities $p_j$ for each known class $j$: $p_j(\vec{x})=P(c_{real}=j|\vec{x})$. I have a training dataset, on which I trained the classifier, and I have a test dataset to test it. In both datasets I have the data in form of $(\vec{x_i}, y_i)$, where $\vec{x_i}$ is the feature vector of each datum, and $y_i$ is the right human-assigned class of each datum (each datum can have only one right class assigned). I can easily evaluate the quality of classification with PR and ROC curves, but I also want to evaluate the quality of probability estimations, that makes my classifier. Now I'm using the next way: 1) I distribute all the $p_j$ in 10 bins: $(0 \dots 0.1), (0.1 \ldots 0.2), \ldots , (0.9 \ldots 1.0)$. 2) Then for each bin I look, how many decisions fall into each bin, and how many of them are right decisions: $P_{bin}=\frac{right-decisions-in- bin}{total-decisions-in-bin}$ 3) Having measured probability curve, I just look at its sum of squared deviations from $f(x)=x$ in each point and try to have a classifier when such sum is minimal. **But I do think that it not the best way of such an evaluation of probability estimation. Cannot you advice me something more reliable? Maybe, there are some standard techniques?** P.S. Some probability curves, measured in my way: ![Some probability curves, measured in my way](http://i.stack.imgur.com/4w1sA.png) The text marks on the plot show (#right decisions in bin / #decisions in bin) for one of the measured probability functions.


title: Residuals plots of residual vs predicted variable

text: I am interested to know why residual plots are plotted with residuals against predicted variable of y and not against y?


title: What makes a GLM estimate the means differently from the actual sample means?

text: _Disclaimer: I'll happily admit to not really knowing what I'm doing when it comes to General Mixed Models but hopefully this thread can bring me a little bit closer to some kind of understanding. If, however, this is an "OMG, you're such a n00b" question which should be closed and never spoken of again, I'll totaly understand :)_ I'm running a couple of General Linear Model (using SPSS) to try to get a better feeling for the method. My dataset consists of a number of stimuli presentations with a dependent scale variable (ranging from 0 to 1) and a number of independent nominal variables for each presentation (for example the gender of the participant, a drug condition, type of stimuli et cetera). Now, when I run a General Linear Model, I can generate a lot of pretty graphs showing both significant and non-significant main and interaction effects for all my data. When including just a couple of factors, the estimated means in these graphs are the same as if I do a descriptive analysis and just compute the means. However, when I start to add a lot of factors, the estimated means starts to deviate from the actual means in my sample. I'm a bit confused here. What is it that makes the GLM estimate the means to be the same as the sample means and when do these estimated means start to deviate from the sample means? Does it have something to do with the number of factors, or is it about the factors themselves (for example, how they divide the data with respect to the other variables).


title: Checking what attributes of a population made objects more likely to do certain action

text: I'm pretty lost at how to approach this, and particularly with regards to the terminology; so please feel free to point to resources, fix my question name, etc. Thanks. I used R back in college, I've got Orange installed. I've heard of MiniTab. I've got Excel installed. I've got access to SQL Server 2008 w/ Analysis Services. That's the breath of my toolkit. [I'm a developer, but that's all I know of for statistical analysis.] I'm marketing a product, and I've got a few datasets. Customers:   * Customer_id   * Campaign they clicked on, and the date.   * Campaign they converted on, and the date.   * The amount of biz they did with us in Q1 2010   * The amount of biz they did with us in Q2 2010   * The amount of biz they did with us in Q3 2010 ... etc. Campaigns:   * Day of the week   * Time of Day   * Specific pitch Results:   * Number of purchases per customer   * Value of purchases per customer I could easily go hogwild, guessing at patterns and testing them to see if there's any accuracy; but it's not a best use of time. What I'd like to know are a few things...   * A ranking of the various pitches in their effectiveness, and an indication of how much the pitch affects the result.   * Same for the day of the week and time.   * And, trickier -- How much of a correlation does their previous business have to do with who converts now? Is it the previous big spenders that convert? Are they the ones buying the most product?   * And finally, of all the attributes above, how do they rank in their importance to conversions?  I assume these would all come with something like an R-squared indicating their 'accuracy.' How would I go about starting? What reading would be essential to understanding this? Is there an approach/tool I can use that'll allow me to do a basic version of these sorts of analyses [and what would they be called] rapidly? I'm just looking for trends and hints of where it might be best to focus our energy; it can be rough. Finally, does anyone know of a resource for basic data mining? analysis? for advertising in particular? Much obliged,


title: Causation implication

text: I recently read an article about how you can increase longevity by sleeping less. This article, like many others I've read, references a statistical study and implies that causation was found between two events(sleeping less and living longer). Well I'm not a statistician, but I'm aware of the common fallacies about causation and I find it very hard to accept the information given in this article and other similar ones. Even further, I can almost never know from an article if the original study found causation or not, without taking the time and digging into it. Although I'd really like to, I'm not going to go into any depth bringing examples of how the causation might be wrong in the given article, because I would most likely be preaching to the choir. Luckily though, to my relief, the original article in NY Times only talks about correlation and even has a disclaimer in the end that basically says that the study doesn't imply causation. Also, I now read that, in accordance with my intuition, it is actually really hard to find and/or test for causation. Reading through the wiki page on the subject I saw two models which seemed to be the most common ones used to find actual links in causation: Granger causality test and convergent cross mapping. But I see, that both of these have some gaps in them. Now, I have a few questions regarding all this. Q1: Are there any models that can, with high accuracy, find actual causation in such widespread studies? If so what are the most common ones? Q2: How often do these widespread studies actually test for causation? E.g. if I read somewhere that a very broad study implies that A causes B, do I have any reason whatsoever to believe that there's something more than correlation behind those claims. Q3: Is there a good technique that can be used by a non-statistician to filter out good statistics from the bad ones in everyday life.


title: Compilation error in JAGS

text: With this data:               y <- c(1.105808,    1.000000,    5.304166,   33.665875,  139.865451, 109.033703,  176.639245,    1.000000,   28.521134,   44.281121 ,150.478570,   18.465554,   85.096431,   81.907537,  124.631226, 1.000000 ,  11.237294 ,  20.480519  , 68.642176 ,  30.047630,  54.051613 , 134.068889  , 72.215041 ,   1.000000 ,  31.254480,   6.226026  , 54.340496 , 161.667352 ,1345.948800 , 147.404744,  192.966923   , 1.000000  ,  1.150755 ,   1.000000  , 99.477430,  72.592598,  107.493014 , 130.201416 , 147.387423 ,   1.000000,   27.534944  ,  8.492657 ,  38.155558,   44.301978 ,  81.938633,   75.026848  ,144.926523 ,   1.000000 ,   1.075801  , 40.001560)          x1 <- c(7.878865,   9.117159 ,  9.998539 , 10.300563 , 12.683197 , 12.185060, 101.346385, 168.814861,   4.769977,   4.803769,   5.990192,   6.469412,   7.557664,   9.781595, 264.102447 ,702.321019 ,  5.663501 ,  6.319843,   6.643405 ,  7.147517 ,  8.154099 ,  8.811370 , 64.089236, 205.728163,   7.218225 ,  6.905615 ,  9.341990 ,  8.554343 , 15.873037 , 15.731554, 227.589294 ,398.435765 ,  6.217681 , 10.498929 , 11.088663 , 10.312797,  11.483123   ,9.276521 ,157.311069 ,391.279665 ,  3.985544 ,  4.389385,  4.663445 ,  4.934453  , 6.622184  , 7.770833 , 21.700911 , 33.470576,   5.405404 ,  6.531107)          x2 <- c(0.4225000, 0.6619411 ,0.5401000, 0.5138000, 0.5109000, 0.6325681, 0.6425919, 0.6466943, 0.4421000, 0.5870430, 0.5254000, 0.5525000, 0.5392000, 0.6330954, 0.6457942, 0.6582531, 0.4039000, 0.6154006, 0.4620000, 0.4875000, 0.5439000, 0.6494275, 0.6423681, 0.6450814, 0.5110000, 0.6060935, 0.5050000, 0.5200000, 0.5535000, 0.6383363, 0.6684222 ,0.6465682, 0.3963000, 0.5914320, 0.4819000, 0.5361000, 0.5886000, 0.6140674 ,0.6240171, 0.6150484, 0.4797000, 0.6211242, 0.5705000, 0.5709000, 0.6144000 ,0.6412593, 0.6611542, 0.6364444, 0.3599000, 0.6375195)      I'm trying to fit a lognormal random effects model in JAGS. Below my JAGS code:               # Lognormal Model     # N municipalities     # M years     # W Betas          model {              for(i in 1:N) {             for (j in 1:M) {                 k <- (i-1)*8 + j                  y[k] ~ dlnorm(mu[k], tau)                 mu[k] <- beta0 + delta[i] + inprod(x[k,], beta[])                 delta[i] ~ dnorm(0, tau2)             }         }              # Prior for betas         beta0 ~ dnorm(mbeta0, precbeta0)         for (l in 1:W) {             beta[l] ~ dnorm(m[l], prec[l])         }              # Prior for precision of Y         tau ~ dgamma(tau.a, tau.b)         sigma <- 1/sqrt(tau)              # Prior for precision of Delta         tau2 ~ dgamma(tau.a2, tau.b2)         sigma2 <- 1/sqrt(tau2)     }      Here my Rjags code:               br <- read.csv("Data/test.csv", header=T, sep=',')     br <- na.omit(br)     br$y <- br$y + 1          # Set up Data     y <- br[, 28]     x <- br[, c(29,30,12,22)]          # Data, Initial Values, and Parameters     N <- length(unique(br$ID))     W <- dim(x)[2]     data <- list(N=N, M=8, W=W, y=y, x=x, m=c(-2,2,4,3), prec=c(.20,.20,.20,.20), tau.a=1, tau.b=2, tau.a2=1, tau.b2=2, mbeta0=1, precbeta0=.01)     inits <- rep(list(list(beta0=0, beta=c(1,1,1,1), tau=1)),5)     param <- c("beta0", "beta", "sigma", "delta")          # Model     sim <- jags.model(file="Code/lognormal2.bug", data=data, inits=inits, n.chains=5, n.adapt=1000)      But I get the following ERROR:               Error in jags.model(file = "Code/lognormal2.bug", data = data, inits = inits,  :        RUNTIME ERROR:     Compilation error on line 11.     Missing values in subset expression of y      Any help?


title: Negative binomial regression does not converge

text: I am trying to use a Negative Binomial regression model for some count data in which the dependent count variable takes on the values 0, 1, 2, 3 or 4. I am using SAS and keep running into the warning: "The relative Hessian convergence criterion of X is greater than the limit of 0.0001. The convergence is questionable." where X ranges approximately from 0.001 to 250. I am just working with univariate models at this point, and this is an issue for all types of independent variables (e.g., whether the independent variable in the model is a continuous variable like age or a categorical variable like income level). I also noticed that the estimate for dispersion is always 0 and followed by either 0's all across for standard error and 95% Wald CIs or by a standard error of <0.01 and missing values for the 95% Wald CIs. The estimates and standard errors for the independent variables equal those if I run a Poisson regression rather than a negative binomial regression (I'm guessing because the dispersion estimate = 0). Is there a way to fix the convergence issue?


title: Non-metric multidimensional scaling with dichotomous independent variables?

text: I have a set of t-rflp fragment-length profiles from an investigation into microbial degradation of hydrocarbons as part of my undergrad dissertation. I want to study the effect of various variables on the microbial community of uncontaminated soils - specifically if the presence of surfactants and soil from oil-contaminated areas influences the microbial community under both aerobic and anaerobic conditions. Essentially the investigation has four t-rflp replicates for each of the following: Aerobic (No surfactant), (No Surfactant + Augmented), (Surfactant), (Surfactant + Augmented) Anaerobic (No surfactant), (No Surfactant + Augmented), (Surfactant), (Surfactant + Augmented) Problem is I've never done anything remotely similar to this before, and my supervisor is away for a few weeks. Is it appropriate to use dichotomous independent variables (+/- surfactant, +/- soil augmentation, Aerobic/Anaerobic conditions) for nMDS? I'm using the Paleontological statistics package (PAST) if it helps. On a related note, is it appropriate to run statistical analysis on DCA axis coordinates of the samples (GLM ANOVA?) Any help will be much appreciated!


title: ANOVA variable that is both BETWEEN & WITHIN

text: Does anyone know how / if can you treat an independent variable as BOTH a between AND within subjects variable (ideally in SPSS ANOVA-repeated measures). That is, we have a repeated measures study where the independent variable is DEVICE (6 levels). The catch is, each respondent will only be exposed to 3 devices and we want to do the analysis across all 6. We will taking 3 measures for each device.


title: What is the best way to Reshape/Restructure Data?

text: I am a research assistant for a lab (volunteer). I and a small group have been tasked with the data analysis for a set of data pulled from a large study. Unfortunately the data were gathered with an online app of some sort, and it was not programmed to output the data in the most usable form. The pictures below illustrates the basic problem. I was told that this is called a "Reshape" or "Restructure". Question: What is the best process for going from Picture 1 to Picture 2 with a large data set with over 10k entries? ![This is a sample, the actual file has over 10k entries. I've color coded to help with clarity](http://i.stack.imgur.com/g54tA.jpg) ![This is what it all needs to be converted into.](http://i.stack.imgur.com/QitXu.jpg)


title: Gamma bridge sampling

text: I am looking into gamma-bridge sampling, as described in for example [1]. What are the benefits of sampling in this method? Is it purely a computational aspect? Are there any problems with the way the function `rgamma()` in `R` samples? From what I can tell, it is an acceptance-rejection method. One of the things that can be better with gamma-bridge sampling is that in increment sampling, round-off errors from the increment will accumulate over the entire process. This will not be the case for gamma-bridge sampling. [1]: Avramidis AN, L'Ecuyer P, Tremblay PA. Efficient simulation of gamma and variance-gamma processes. Proceedings of the 2003 winter simulation conference. IEEE 2003, p. 319-326


title: Papers about quantitative association rules in data mining

text: I'm writing my diploma thesis about the interaction between warehouse management and data mining (specifically, association rules). I started the topic data mining by considering Boolean association rules, created by the a-priori algorithm (for example, market basket analysis). Then I found out about quantitative association rules. They consider not only the question 'what is together in the market basket', but also how much of one item is combined with how much of the second and so on. Maybe the most important author (and even sort of founder) of data mining is Rakesh Agrawal. He wrote a paper (pdf) about mining quantitative association rules. After I found that and considered it, I wanted to find newer papers or literature about this topic, because Agrawal released it 1996. I could find only 3-4 papers about quantitative association rules (respectively quantitative data mining) and don't understand why this topic is not discussed more. Do you know why it is so hard to find any papers considering it? Why aren't there advancement from Agrawal's idea or just didn't I found them? Also there exists a lot of software for data mining (for example WEKA) but they don't have the opportunity to mine quantitative association rules, only Boolean ones. So I want you to ask, if you know any software for mining quantitative association rules?


title: Why do we need clearly defined variables as response/dependent variables?

text: This is not a homework question, just sort of a mental block. Why do statistical analyses require clearly defined variables as response/dependent variable? Why can't my response variable be something [vague] like "Likely High Value Customer" or a yes/no question?


title: How to quantify redundancy of features?

text: I have three features that I use to solve a classification problem. Originally, these features produced boolean values, so I could evaluate their redundancy by looking at how much the sets of positive and negative classifications overlap. Now I have extended the features to produce real values (scores) instead, and I would like to analyze their redundancy again, but I am at a complete loss on how to do that. Can anyone provide me with a pointer or idea on how to go about that? I know this question is very vague, that is because I do not have a very strong grasp of statistics. So, if you do not have an answer for me, maybe you have some questions that can help me understand better myself. **Edit:** I am currently browsing Wikipedia on the subject, I have the feeling that what I want is a correlation coefficient, but I am still unsure if this is the right approach, and which of the many available coefficients is appropriate. **Edit 2:** In the boolean case, I first created for each feature the set of samples for which it was true. Then, the correlation between two features was the size of intersection of these sets over the size of the union of these sets. If this value is 1, they are completely redundant, because always the same. If it is 0, they are never the same.


title: How can I effectively summarize and visualize time series of employee activities?

text: I am managing many people entering data into a database. I have a log of user, date, time, table, and action that each person makes:               records <- data.frame(user = c('bob', 'bob', 'jane', 'jane', 'bob', 'bob', 'bob', 'jane', 'jane', 'bob'),                           date = c("2010-06-24", "2010-06-28", "2010-06-29", "2010-06-30", "2010-07-01", "2010-07-02", "2010-07-05", "2010-07-06", "2010-07-07", "2010-07-09"),                            time = c("01:40:08", "01:40:18", "01:40:28", "01:40:37", "01:40:44", "01:40:52", "01:40:59", "01:56:26", "02:16:37", "03:55:06"),                           table = c(rep('table1',5), rep('table2',5)),                           action = c('create', 'create', 'create', 'update', 'create', 'update', 'update', 'create', 'create', 'create'))       For a non-trivial example, the actual `records` dataframe with 10,000 entries can be downloaded as an .Rdata file here, and then:               load('records.Rdata')     library(ggplot2)     qplot(date, table, data = records, color = user, geom='jitter')      How can I visualize, overall and for each table:   1. the amount of time each person works per week   2. the type and number or frequency of actions that they made. ?


title: Question about combining hazard ratios - Maybe Simpson's paradox?

text: I’m reviewing an article, and can’t give details but here is the situation, and it’s got me puzzled Patients were divided into 4 categories (call them A B C and D), which were exhaustive and exclusive. Adjusted hazard ratios were computed for these four groups for all patients and for two subgroups of patients (call them Lung and Heart). The two subgroups were also exhaustive and exclusive. For the full sample, the hazard ratios for B C and D (compared to A) were 5.08, 4.39, and 1.81. For the lung group, they were 1.15, 1.16 and 1.07. For the heart group they were 2.55, 1.69 and 1.28. That is, they were lower for both subgroups than for the whole population. The lung group and heart group are not the same size; from the info provided, it’s not clear whether A B C and D are equally divided across heart and lung groups. So …. This seems like Simpson’s paradox. But I’ve not seen that term used for survival analysis. I don’t see why it could NOT be so applied. My feeling here is to suggest that heart and lung should ONLY be presented separately, but I’m not completely sure. (Excuse me if I am missing something obvious – dealing with some family health problems and not sleeping that well) Thanks for any insights or thoughts


title: Error creating parser in readHTMLTable

text: readHTMLTable seems pretty robust, but when I try to use it on this page, I get an error. Any ideas what it means or how I could get around it? The page's biggest table is nested inside another table... is that a problem?               > tables<-readHTMLTable(myURL, header=NA,a.data.frame=TRUE)     Error in htmlParse(doc) :        error in creating parser for http://www.ncbi.nlm.nih.gov/genomes/lproks.cgi


title: GLMNET or LARS for computing LASSO solutions?

text: I would like to get the coefficients for the LASSO problem $$||Y-X\beta||+\lambda ||\beta||_1.$$ The problem is that glmnet and lars functions give different answers. For the glmnet function I ask for the coefficients of $\lambda/||Y||$ instead of just $\lambda$, but I still get different answers. Is this expected? What is the relationship between the lars $\lambda$ and glmnet $\lambda$? I understand that glmnet is faster for LASSO problems but I would like to know which method is more powerful?


title: Meta-analysis of 3 regions' data for 5 years

text: Can we do a meta-analysis of data of 3 regions. A particular disease was treated with the same treatment but implemented thoroughly in 2 regions and not so thoroughly in 1 region. How to proceed with analysis? What is the best way to analyse these data. The data are of a disease (cases and deaths) as follows:                         Region 1      Region 2         Region 3                        Cases   Dead    Cases   Dead     Cases  Dead    Total cases     Total dead cases     2006    2320    528     1484    108       73    3       3877            639     2007    3024    645     1592    75        32    1       4648            721     2008    3012    537     1920    53         3    0       4935            590     2009    3073    556     1477    40       246    8       4796            604     2010    3540    494     1460    26       138    1       5138            521     Total  14969    2760    7933    302      492    13     23394            3075      _**Updated_** These are unpublished data.


title: How to use triple exponential smoothing to forecast in Excel

text: I have been burdened with the task of coming up with a forecast plan for my company. I have no experience and am VERY new to the whole forecasting scene. As of right now, my company has no plans of investing in any forecasting software so my only tool is Excel. I've tried to do some research online myself and it seems that this triple smoothing method would be a great asset, but I'm a little confused and I guess I don't really understand the equations. Below I have provided 3 years worth of sales for one item. We forecast in periods (4 weeks = 1 period). So there are 13 periods in one year. When we forecast, we have to forecast out 6 periods into the future, please help me use the triple smoothing technique to accomplish this.               Period 10 2009  69,088     Period 11 2009  83,400     Period 12 2009  75,735     Period 13 2009  79,526     Period 01 2010  81,005     Period 02 2010  94,013     Period 03 2010  90,567     Period 04 2010  94,568     Period 05 2010  101,687     Period 06 2010  93,540     Period 07 2010  84,249     Period 08 2010  91,280     Period 09 2010  78,531     Period 10 2010  89,465     Period 11 2010  83,341     Period 12 2010  87,106     Period 13 2010  65,636     Period 01 2011  79,632     Period 02 2011  89,722     Period 03 2011  87,483     Period 04 2011  99,228     Period 05 2011  113,215     Period 06 2011  96,057     Period 07 2011  95,475     Period 08 2011  92,466     Period 09 2011  103,529     Period 10 2011  94,515     Period 11 2011  76,146     Period 12 2011  81,736     Period 13 2011  80,174     Period 01 2012  81,437     Period 02 2012  102,695     Period 03 2012  120,775     Period 04 2012  97,058     Period 05 2012  119,921     Period 06 2012  102,311     Period 07 2012  109,498     Period 08 2012  110,318     Period 09 2012  98,103          Period 10 2012       Period 11 2012       Period 12 2012       Period 13 2012       Period 01 2013       Period 02 2013       Period 03 2013


title: Measuring representativeness of a sample

text: In the task of language identification, I would do the following:   1. Take a sample of my data   2. Prepare the ground truth    3. Train my classifier on this sample data   4. Test classifier accuracy on the other part of my data that I did not sample (and possibly cross-validate)  Let us say I am working on identifying the language of a tweet. Getting the ground truth of this is hard (either I need to manually annotate the tweets with their correct language or crowd-source this task). Therefore, what I did was to train my classifier on already available text (such as eBooks, news articles) of different languages. This is **plain wrong** - I am training my classifier for one domain (large well-written text) but testing it on an other domain (small potentially-garbled up text). Surprisingly, this worked with a 90% accuracy on 400 hand-tagged tweets. How will I know what is the sample size I need to pick (which is 400 in my case) that can show that the classifier will still behave at about 90% for a million tweets?


title: Change detection in circular data

text: Suppose I have a data set that represents circular data measured in degrees:               x <- c(rnorm(1000, 0, 10), rnorm(700, 110, 3), rnorm(1100, 230, 5)) %% 360      The R package `circular` provides a very nice way to represent that data, and a basic tool for detecting change points in it:               library(circular)     x <- circular(x, units='degrees')     cp <- change.point(x)      However, this particular algorithm is limited, because it's (in my experience) relatively inefficient, and it's limited to finding one change point at a time, so if multiple change points are present, a recursive approach is needed. This causes some difficulty in deciding when to terminate. If a linear change point algorithm is used, it will have difficulty with `x[1:1000]` because some values will be close to `0` and some close to `360`. For linear data, I like the 'PELT' algorithm of Killick, Fearnhead, and Eckley (2011) implemented by the R package `changepoint`'s `cpt.mean()` function. It's fast and seems to be pretty reliable. Has anyone looked at adapting this method to circular data? Or other recommendations?


title: Determining sample size

text: If I have a population and want to draw samples from there, what are the steps to determine the number of samples to draw from the population? In other words, if I have population N and want to draw k samples from them so that the samples are representative of the entire population, how can I determine k? I am having some insights from here but in the end of this document I find- > You can still use this formula if you don’t know your population standard > deviation which is my case (i.e., I do not know the population variance). How can I use that formula if I don't know the population standard deviation/ variance? n.b. I posted this question to that forum but still no luck.


title: Can you compare Cramer's V and Spearman's $\rho$ values?

text: I did correlation tests between a number of variables consisting of both interval and categorical scales ($n = 1274$, one sample). I used $\chi ^2$ test for categorical vs ordinal and Spearman's $\rho$ for interval vs ordinal data. The correlation between climate zone (zone $1, 2, 3, 4$) and daily exposure to air-conditioning (low, medium, high) is $0.35$ (Cramer's V). The correlation between exposure to AC and pro-environmental attitude score ($1-5$, interval) is $0.21$ (Spearman's $\rho$). Can I compare the strength of the relationship if I used different correlation tests, say the relationship between climate and exposure to AC is stronger than that between the pro-environmental attiudes and exposure to AC?


title: Sample size determination

text: Two tools completely agree on categorizing 85k sentences (my entire population). Now, I want to take a sample of sentences and give that to two humans. This is to check whether the tools agreed genuinely or both of them together miscategorized the 85k sentences. Now, I need to know how many sentences should I pick as a sample out of 85k sentences to check it statistically. Is there any formula or procedure to find out?


title: K-means Mahalanobis vs Euclidean distance

text: I currently am trying to cluster "types" of changes on bitemporal multispectral satellite images. I applied a thing called a mad transform to both images, 5000 x 5000 pixels x 5 bands. Each band is a "variable" as it is radiance information from a different spectrum of light. This transform is basically equivalent to PC applied to the substraction of both images. Naturally i can get up to 5 mad components. Now I would like to find ... excatly this... types of change on these components. If I use K-means on the components I would use an euclidean distance but I just wanted to know what could be the gain in using a mahalanobis distance if there is any... Just curious. Greetings and sorry for posting such a TERRIBLE quesiton before. I actually just wanted a type of overview of both options. Julián.