title: Understanding measure concentration inequalities

text: In the spirit of this question Understanding proof of a lemma used in Hoeffding inequality , I am trying to understand the steps that lead to Hoeffding's inequality. What holds the most mystery for me in the proof is the part where exponential moments are computed for the sum of i.i.d variables, after-which Markov's inequality is applied. My goal is to understand: Why does this technique give a tight inequality, and is it the tightest we can achieve? A typical explanation refers to the moment generating properties of the exponent. Yet, I find this too vague. A post in Tao's blog, http://terrytao.wordpress.com/2010/01/03/254a-notes-1-concentration-of- measure/#hoeff, might hold some answers. With this goal in mind, my question is about three points in Tao's post which I am stuck at and which I hope could give insight once explained.   1. Tao derives the following inequality using the k-th moment$$\displaystyle {\bf P}( |S_n| \geq \lambda \sqrt{n} ) \leq 2 (\frac{\sqrt{ek/2}}{\lambda})^k. \ \ \ \ \ (7)$$ If this is true for any k, he concludes an exponential bound. This is where I'm lost. $$\displaystyle {\bf P}( |S_n| \geq \lambda \sqrt{n} ) \leq C \exp( - c \lambda^2 ) \ \ \ \ \ (8)$$    2. Hoeffding's lemma is presented: Lemma 1 (Hoeffding’s lemma) Let ${X}$ be a scalar variable taking values in an interval ${[a,b]}$. Then for any ${t>0}$, $$\displaystyle {\bf E} e^{tX} \leq e^{t {\bf E} X} (1 + O( t^2 {\bf Var}(X) \exp( O( t (b-a) ) ) ). \ \ \ \ \ (9)$$ In particular $$\displaystyle {\bf E} e^{tX} \leq e^{t {\bf E} X} \exp( O( t^2 (b-a)^2 ) ). \ \ \ \ \ (10)$$ The proof of Lemma 1 begins by taking expectation over the taylor expansion $\displaystyle e^{tX} = 1 + tX + O( t^2 X^2 \exp( O(t) ) )$ .Why can the expansion be bounded by that quadratic term? and how does equation 10 follow?   3. Finally, an exercise is given:   Exercise 1 Show that the ${O(t^2(b-a)^2)}$ factor in (10) can be replaced with ${t^2 (b-a)^2/8}$, and that this is sharp. This would provide a much shorter proof than the one in Understanding proof of a lemma used in Hoeffding inequality , but I do not know how to solve this. Any further intuitions\explanations about the proof of the inequality or the reason we cannot derive a tighter bound are definitely welcome.


title: Problem solving a simple probability exercise

text: **I'm currently studying a Statistics course and I'm trying to solve a simple probability exercise:** A production facility employs:   * 20 workers on the day shift   * 15 workers on the swing shift   * 10 workers on the grave- yard shift A quality control consultant is to select 6 of these workers for in-depth interviews. Suppose the selection is made in such a way that any particular group of 6 workers has the same chance of being selected as does any other group (drawing 6 slips without replacement from among 45). The question is: _What is the probability that at least two different shifts will be represented among the selected workers?_ **My answer:** I'm going to answer this question using two methods. I have a problem with the second method. Though, the first method is more recommended by instructors; when practicing I try to use both methods to make sure that my answer is correct and that I understand the material. **Method 1:** Calculating **_the probability of 6 workers being selected from 45 workers_** , then subtracting it with **_the probability of 6 workers being selected exactly from one different shift_**. P(A) = The probability of 6 workers being selected from 45 workers is obviously one: (45P6 / 45P6). P(B) = The probability of 6 workers being selected from exactly one different shift is: 0.0054 (15P6 + 10P6 + 20P6 / 45P6) P(Solution) = P(A) - P(B) = 1 - 0.0054 = 0.9946. The solution 0.9946 is correct (verified from the answer key) and I have no problem with this method. **Method 2:** Calculating **_the probability of 6 workers being selected exactly from two different shifts_** and summing it with the **_the probability of 6 workers being selected exactly from three different shifts_**. P(C) = The probability of 6 workers being selected exactly from two different shifts: 0.29 (25P6 + 35P6 + 30P6 / 45P6). Note that 15+10=25 15+20=35, 20+10=30. P(D) = The probability of 6 workers being selected exactly from three different shifts is: 1 (45P6/45P6). Summing P(C) + P(D) will obviously give a wrong answer since the result will be above 1, which is clearly false. What did I do wrong in method 2?


title: Is it valid to preprocess data with PCA before running Locally Linear Embedding (LLE)?

text: I have some very high dimensional data, and performing Locally Linear Embedding (LLE) is very time consuming. I also have to perform several LLEs, with varying parameters, to compute the optimal number of neighbors via Spearman's Rho (a la Karbauskaite et al, _see_ citation). However, the data are highly linearly dependent--less than 10% of the features are necessary to preserve 98% of the variance. This does not preclude the data still being embedded in a lower dimensional manifold. Thus, my question is, is it valid to perform principal component analysis (PCA) to preserve the vast majority of the variance, followed by LLE, or will that radically change the results of the LLE? This would significantly reduce my compute time--so I'm hoping that can be done. As an aside, the Spearman's Rho is computed based on pairwise distances-- although this requires using the Geodesic distance across the manifold, NOT the Euclidean distances. Is there a way to compute Geodesic distances based on Scikit-Learn's implementation of PCA? Citation: Karbauskaitė, Rasa, and O. Kurasova G. Dzemyda. "Selection of the number of neighbours of each data point for the locally linear embedding algorithm." Information Technology and Control 36.4 (2007): 359-364.


title: Regression techniques for capped independent variables

text: Suppose you are trying to predict height ($y$) based on age ($x$). A straight linear regression won’t work well, since humans stop growing at a certain age. This suggests instead trying to model height as something like: $y = w\cdot \mathrm{squash}(x, a, b) + w_0$ where $\mathrm{squash}(x,a,b) = \max(a, \min(x,b))$ In general, when trying to predict some dependent variable $y$ from dependent variables $x_1, x_2, ..., x_n$, some of which you expect to display "squash"-characteristics like this, it seems a good idea to try to model $y$ by something like: $y = w_0 + \sum_{i=1}^n w_i\cdot\mathrm{squash(x_i, a_i, b_i)}$ I'm wondering what regularized techniques exist to estimate $\vec{w}, \vec{a}, \vec{b}$.


title: Permutation test comparing a single sample against a mean

text: When people implement permutation tests to compare a single sample against a mean (e.g., as you might do with a permutation t-test), how is the mean handled? I have seen implementations that take a mean and a sample for a permutation test, but it is unclear what they're actually doing under the hood. Is there even a meaningful way to do a permutation test (e.g., t-test) for one sample versus an assumed mean? Or, alternatively, are they just defaulting to a non-permutation test under the hood? (e.g., despite calling a permutation function or setting a permutation test flag, defaulting to a standard t-test or similar function) In a standard two-sample permutation test, one would have two groups and randomize the assignment of labels. However, how is this handled when one "group" is an assumed mean? Obviously, an assumed mean has no sample size in and of itself. So then, what is the typical way of working the mean into a permutation format? Is the "mean" sample assumed to be a single point? A sample of equal size to the sample group? An infinitely-sized sample? Given that an assumed mean is, well, assumed- I'd say it technically has either infinite support or whatever support you want to assume for it. However, neither of these are very useful for an actual calculation. An equal- sized sample with values all equal to the mean seems to be what is done sometimes with some tests (e.g., you just fill in the other half of the pairs with the assumed location). This makes a bit of sense, as it's the equal- length sample you'd see if your assumed mean was correct with no variance. So my question is this: In practice, do people actually emulate permutation- test style label randomization when the second set is a mean (or similar abstract assumed value)? If so, how do people handle label randomization when they do this?


title: How to handle new observations on HMM decoding?

text: I'm implementing the HMM algorithms described in Rabiner's tutorial. But there is several issues to considered when we apply HMM for real problems. One of this problems is how to consider new observations that hasn't emission probabilities. Anyone knows a way of considering a distribuition for new observations? I believe that the information about previous state can be useful for inferir the state of the the new observation.


title: How to get SGD to reach global optimal point in logistic regression?

text: I am trying to write a tool which involves implementing logistic regression. With the batch gradient descent method, the convergence is guaranteed as it is a convex problem. However, I find that with the stochastic gradient decent method typically converges to some random points (i.e., not very close to the minimum point resulted from the batch method). I have tried different ways of decreasing the learning rate, and different starting points of weights. However, the performance (e.g., accuracy, precision/recall, ...) are comparable (to the batch method). I understand that this is possible, since SGD (stochastic gradient descent) uses an approximation to the real cost each step. I guess this matters since otherwise the interpretation of the weights would not make much sense even the accuracy is comparable.   * Does it matter that SGD uses an approximation the real cost at each step?    * If it matters, how can it be made to converge or get close to the global optimal point?


title: Test the combined hypothesis "X is correlated with Y and is uncorrelated with Z"

text: If I have 3 variables, $X$, $Y$ and $Z$, what statistical approach might I use to ask the question "is in population X correlated with $Y$ but is uncorrelated with $Z$"? Obviously, I can assess if $X$ is correlated with $Y$ and separately if $X$ is correlated with $Z$ and have two p-values:               p1 <- cor.test(X,Y)$p.value     p2 <- cor.test(X,Z)$p.value      But I'm not sure how I'd go about sort of "combining" these two results **into one p-value** of my combined hypothesis. Or maybe there is a totally different approach for this question?


title: Adding up events in time series forecasting

text: Let us say that we have an event - variable (1/ 0) that denotes the occurence of an event on a daily basis e.g. a strike. Let us now say that we have a continuous variable (sales) that that we want to forecast. The interval on which the variable is measured is monthly. We want to incorporate the event variable in a model that will be used for forecasting sales but the event variable is measured in a more detailed interval. If there were 5 days in a month that a strike took place can we say that the event variable value for that month is 5? More generally can we create a derived event variable by aggregating the values of the original dalily event variabe? In this way we will have an independent variable with value ranging from 0 (no strikes occured the specific month) to 30 (if a strike occured every day in the specific month)? Does this process stand from a statistical point of view and is this variable going to have any predictive power? Thnaks in advance, Andreas * * * Thnaks for your answer. Let us make things less complicated assuming that the dependent variable is electricity consumption and the event is a day with extremely low average temperature (day with extreme cold). With this case we don;t care about weekends and the rest of the problems that are coused in the case of strikes. Let us say that we have electricity consumption per month period so 12 observations per year. Let us also say that we have recorder the days woth extereme cold daily by using a binary variable (so approximately 365 observations per year. What would you do if you wnated to incorporate the daily event in the model of forecasting the monthly series. Would you add up the 1, would you take the percentage of days over the month that exhibit the extreme cold or what else? What do you htink is the better solution?


title: Train NN or SVM to classify stock signals

text: I am applying Neural network and SVM to predict buy-hold - sell signals. I have trained nn and SVM in R. I used nnet function to train NN and svm to train SVM. I provided 20,000 data points to train and 2000 data points to test. training data set contain list of 10-15 technical indciators and buy- sell-hold signals. The problem I am having is it is not predicting the buy -sell - hold signals with good accuracy on the testing data. I have used sigmoid function in nnet and radial function in SVM. Any suggestion, how to improve the accuracy of the prediction?


title: Causal identification and penalized splines

text: I just got a rejection from an economics journal. Among the reasons cited for rejection were: > the benefits of using the semi-parametric method are not clearly brought out > compared to alternative simpler techniques with clean identification of > causal relationships It is certainly possible that I could have done a better job of motivating the methodology to a bunch of economists who generally stick to OLS. But have I violated "clean identification"? Please judge for yourself and let me know what you think: My main estimating equation is $$ y_{it} = \alpha_i + \beta_1 T_{it} + f\left(\begin{array}{l}Z_{it}\\\ Z_{it} \times T_{it} \\\ Z_{it}\times T_{it} \times X_t\end{array} \right) + \beta_2X_t + \epsilon_{it} $$ $Z$ is continuous, $X$ and $T$ are binary. I can justifiably assume that $$ E[\epsilon|\alpha,T] = 0 $$ Which is to say that the coefficient on $T$ is unbiased conditional on individual-level dummy variables ("fixed effects" in econometrics-speak). When I include continuous variable $Z$, I am simply looking at heterogeneity in estimated treatment effects over gradients of $Z$. So the average causal effect of the treatment $T$ is a an average of $\hat\beta_1 + \hat f_{Z\times T}$ for the various levels of $Z$ that I observe. The model is etimated by penalized quadratic splines (e.g.: Ruppert et al. 2003). Specifically: $$ y = \beta_0 +X'\beta + \displaystyle\sum_{1}^p (Z^{p})'\gamma + \displaystyle\sum_{j=1}^{\\#vars} \displaystyle\sum_{k=1}^{\\# knots_j}\delta_{jk}\left(\left(Z_j - \kappa_{jk} \right)^p \times \left(Z_j > \kappa_{jk} \right)\right) + \epsilon $$ This is solved by $$ \left[\begin{array}{c} \hat\beta\\\ \hat\gamma \\\ \hat \delta \\\ \end{array}\right] = (C'C + \lambda^{2p}D)^{-1}C'y $$ where $C$ includes the parametric terms and the knot terms, and where the ridge penalty $\lambda$ _only applies to the knot terms_ , and is chosen to minimize AIC. (I can't do full justice to the methodology -- see Ruppert et al, or Simon Wood's textbook on GAM's). Of course, I use these semiparametrics because I don't want to impose unfounded functional forms onto my data. Doing so would quite naturally bias my estimates as much as imposing a logarithmic fit onto a sinusoidal function would bias my estimates. But is there something inherent in penalized splines as I've described them that would inherently render the following statement untrue? $$ E[\hat\beta_1] = \beta_1 \text{ iff } E[\epsilon|\alpha,T] = 0 $$


title: First steps learning to predict financial timeseries using machine learning

text: I am trying to get a grasp on how to use machine learning to predict financial timeseries 1 or more steps into the future. I have a financial timeseries with some descriptive data and I would like to form a model and then use the model to predict n-steps ahead. What I have been doing so far is:               getSymbols("GOOG")          GOOG$sma <- SMA(Cl(GOOG))         GOOG$range <- GOOG$GOOG.High-GOOG$GOOG.Low          tail(GOOG)                          GOOG.Open GOOG.High GOOG.Low GOOG.Close GOOG.Volume GOOG.Adjusted     sma range     2013-05-07    863.01    863.87   850.67     857.23     1959000        857.23 828.214 13.20     2013-05-08    857.00    873.88   852.91     873.63     2468300        873.63 834.232 20.97     2013-05-09    870.84    879.66   868.23     871.48     2200600        871.48 840.470 11.43     2013-05-10    875.31    880.54   872.16     880.23     1897700        880.23 848.351  8.38     2013-05-13    878.89    882.47   873.38     877.53     1448500        877.53 854.198  9.09     2013-05-14    877.50    888.69   877.14     887.10     1579300        887.10 860.451 11.55      Then I have fitted a randomForest model to this data.               fit <- randomForest(GOOG$GOOG.Close ~ GOOG$sma + GOOG$range, GOOG)      Which seems to fit surprisingly well:               > fit          Call:      randomForest(formula = GOOG$GOOG.Close ~ GOOG$sma + GOOG$range,      data = GOOG)                     Type of random forest: regression                          Number of trees: 500     No. of variables tried at each split: 1                    Mean of squared residuals: 353.9844                         % Var explained: 97.28      And tried to use it to predict:               predict(fit, GOOG, n.ahead=2)      But this prediction ofc did not work. I try to predict the Close, should I lag the other variables by as many steps as I want the prediction, before fitting the model? Probably a lot of other stuff I should take into account as well but these are really my first steps trying out machine-learning. Thankful for any tips!


title: Probability Distribution

text: Good evening everyone, I'm faced with another question which I can't seem to be able to find a solution online. Thought I get some directions and advice here. I have the following distribution, which represents the case of the number of events that happen in a single day.   * When (x=0), P(X=x) = 0.2   * When (x=1), P(X=x) = 0.4   * When (x=2), P(X=x) = 0.3   * When (x=3), P(X=x) = 0.08   * When (x=4), P(X=x) = 0.02 Suppose I want to find the probability that the sample standard deviation is 0 if I randomly selected 4 days, how should I proceed? I have worked out the   * mean of x to be 1.32   * variance of x to be 0.8976   * standard deviation to be 0.947 Appreciate advice and direction please


title: The difference of two discrete uniform variables taken modulo. Is still uniform?

text: I have a uniform distribution that generates from the $\mathbb{Z}_q$ (so I have integers from the interval [$-q/2$, $q/2$). Then I subtract one from the other and take the result modulo $q$ (so again I have an integer from the interval above). How to prove (or disprove) that the resulting distribution remains uniform in $\mathbb{Z}_q$?


title: What is the relation of the significance level alpha to the type 1 error alpha?

text: In statistical hypothesis testing we decide on and _set_ the acceptable probability of error or **significance level** α (alpha) to a value that fits our theory. Traditionally alpha is .1, .05, or .01. When we _calculate_ the power function g of the parameter we test for, we recieve the distribution of the probability of two errors: the **Type 1 error** α (alpha) and the Type 2 error β (beta). In a graphical representation of this function, alpha is the value below the graph, beta is the value above the line: α = g(p) and β = 1 - g(p), with p being the parameter we are interested in. Here is an example: ![sample power function](http://i.stack.imgur.com/dNGUU.png) The red line is αmax for H0: p ≤ 0.4 and H1: p > 0.4; the blue line is β for a sample p̂ = 0.5 **How do the probability of error (alpha) and the Type 1 error (alpha) relate to each other?** Or am I just getting confused over two unrelated values having the same name (alpha)?


title: Population genetics: analysis of allele frequency distribution along environmental gradients

text: Is there a software or a R package that does a simulation of spatial distribution of alleles? I'm interested in how a distribution of alleles of functional loci is shaped by a pattern of environmental variables, such as temperatures and water quality. I have data on geographical patterns of environmental variables, population connectivity (Fst)(microsatellite data), geographic distances, population size and allele frequency of functional loci at 17 sites. I want to know allele frequency of the functional loci at other sites in a map for which we don't have allele frequency data. Maybe using existing allele frequency data (from functional loci) and other information mentioned above, allele frequency at other sites could be predicted with an assumption that loci are under selection of targeted environmental variables? Thanks!


title: Analysing possible effect of categorical response on ordinal response?

text: I have recently begun to analyse data from a survey where participants rated 64 images of models on a Likert scale and then selected whether they consider the model in-group or out-group. I am trying to discern whether the fact that participants consider the model in-group or out-group has an impact on the ratings the participants gave for each image. So my data set consists of 64 Likert-style questions rated `1-7` for each participant and 64 categorical (in-group or out-group) questions and I am completely unsure how to begin analysing them. Which sort of test would be appropriate for this form of data?


title: How would you explain generalized linear models to people with no statistical background?

text: I always have a hard time explaining statistical techniques to audience with no statistical background. If I wanted to explain what GLM is to such audience (without throwing out statistical jargon), what would be the best or most effective way? I usually explain GLM with three parts -- (1) the random component which is response variable, (2) the systematic component which is linear predictors, and (3) the link function which is the "key" to connecting (1) and (2). Then I would give an example of linear or logistic regression and explain how the link function is selected based on the response variable. Hence it acts as the key connecting two components.


title: Summing bootstrap-derived confidence intervals?

text: I have a couple things to add together- two different pools of material, and I need to know the total material. Because of various things related to the sampling, I would like to bootstrap the sampling data from each pool. So I generate a list of 100 masses for one pool, and 50 masses for the second pool. Plus, I get the confidence intervals built off the bootstrap. I can just add the mean, to get the mean total mass, but what about the CIs? I'm still learning bootstrapping, obviously, so sorry if it's a simple question. At first, I figured I could combine them similar to SD's, like sqrt(x^2 + y^2), but that doesn't work. Does simple addition work? Or do I need to build the CI's from scratch using the standard errors from each pool?


title: Are these possible LeveneTest results?

text: I have experimental data of response times that differ according to the factor `TrialType`: `Regular` or `Random`. I also have two groups, `A` and `B`, which I expect to behave the same (they were created for counterbalancing purposes). I want to test the equality of variance assumption in my data. I am pretty sure the variances are _not_ equal, but I want to show it with a test anyway. My data is contained in a long data.frame called `AB`:               summary(AB)     Block          Trial         TrialType      Subject          variable      value             Block1:1930   Min.   :  1.0   Random :5889   Length:9852        RT:9852   Min.   : 348.0        Block2:1974   1st Qu.:128.0   Regular:3963   Class :character             1st Qu.: 595.0        Block3:1982   Median :252.0                  Mode  :character             Median : 692.0        Block4:1983   Mean   :251.9                                               Mean   : 765.7        Block5:1983   3rd Qu.:376.0                                               3rd Qu.: 834.0                  Max.   :500.0                                               Max.   :4412.0        Group         A:4929        B:4923        I tried the Levene Test using both Group and TrialType as grouping factor (the second argument):               leveneTest(AB$value,AB$Group)          leveneTest(AB$value,AB$TrialType)      As you can see, I am still unclear whether the second argument should be `Group` (`A` and `B`) or `TrialType` (`random` and `regular`). I just tried both. For group, the p-value it was 8.479e-11, and for TrialType, p was 2.413e-12. Also, the degrees of freedom for both cases is 1, and the F values are very high: 42,23 for Group and 49.24 for TrialType. So the 3 values I am getting from these tests seem highly unusual. Therefore I suspect I might be doing something wrong. Any comments or/and help is much appreciated. Also, any tips on how to report these results would be most welcome.


title: How do you select variables in a regression model?

text: The traditional approach to variable selection is to find variables that contribute the most to predicting a new response. Recently I learned of an alternative to this. In modeling variables that determine the effect of a treatment--as for example in a clinical trial of a pharmaceutical--the variable is said to be _qualitatively interacting_ with treatment if, leaving other things fixed, a change in that variable can create a change in which treatment is most effective. These variables are not always strong predictors of the effect but may be important for a physician when deciding on treatment for individual patients. In her PhD thesis Lacey Gunter developed a method for selecting these qualitatively interacting variables that could be missed by algorithms that base selection on prediction. Recently I have worked with her on extending these methods to other models including logistic regression and Cox proportional hazard regression models. I have two questions:   1. What do you think about the value of these new methods?   2. In the case of the traditional methods what approach do you prefer? Criteria such as AIC, BIC, Mallows Cp, F tests for entering or dropping variables in stepwise, forward and backward... * * * The first paper on this came out in Gunter, L., Zhu, J and Murphy, S. A. (2009). _Variable selection for qualitative interactions_. **Statistical Methodology** doi:10, 1016/j.stamet.2009.05.003. The next paper appeared in Gunter,L., Zhu, J. and Murphy, S. A. (2011). _Variable selection of qualitative interactions in personalized medicine while controlling the familywise error rate_. **Journal of Biopharmaceutical Statistics** 21, 1063-1078. The next one appeared in a special issue on variable selection Gunter, L., Chernick, M. R. and Sun, J. (2011). A _simple method for variable selection in regression with respect to treatment selection_. **Pakistan Journal of Statistics and Operations Research** 7: 363-380. You can find the papers at the journal websites. You may have to purchase the article. I might have the pdf files for these articles. Lacey and I have just completed a monograph on this topic which will be published as a SpringerBrief later this year.


title: A non-statistician reference for confidence-interval Vs credibility-interval interpretation

text: As a beginner statistician, discussing the need to be accurate with the interpretation of statistical results with non-statistician is not an easy task. In particular, I am trying to convince some colleagues not to dive into the common mis-interpretation of confidence-interval (maybe I should not?). For some problem, there is some illustrating and highlighting (and funny) paper as this one (http://www.jsur.org/ar/jsur_ben102010.pdf) about the need of multiple testing correction for whose I extract a small citation: > [...] we completed an fMRI scanning session with a post-mortem Atlantic > Salmon as the subject. The salmon was shown the same social perspective- > taking task that was later administered to a group of human subjects." that is a terrific argument to encourage practioner to use multiple testing corrections. Are you aware of such a paper about the common mis-interpretation of confidence interval (that are typically interpreted as credibility-interval) ? or any other than can help me in my arguments.


title: How to incorporate costs (into logit model) of false positive, false negative, true positive, true negative if they are different costs?

text: How to incorporate costs (into logit model) of false positive, false negative, true positive, true negative responses, if they are different costs ? Is it possible to do that on the level of likelihood function ? Edition : I see know that likeligooe function could be quite easily modified to incorporate costs, but then likelihood function become discontinuous : y_i == 1 oraz f(x_i*B) > 0.5 cost = cost11   y_i == 1 oraz f(x_i*B) < 0.5 cost = cost10   y_i == 0 oraz f(x_i*B) > 0.5 cost = cost01   y_i == 0 oraz f(x_i*B) < 0.5 cost = cost00 non-modified likelihood function :   f(x_i*B)^y_i * ( 1 - f(x_i*B) ) ^ (1 - y_i) modified likelihood function : (cost11*positive(f(x_i*B)-0.5) + cost10*negative(f(x_i*B)-0.5) ) ^ y_i   *   (cost00*negative(f(x_i*B)-0.5) + cost01*positive(f(x_i*B)-0.5) ) ^ (1 - y_i) where : positive(x) = 1 if x > 0   positive(x) = 0 if x < 0 negative(x) = 1 if x > 0   negative(x) = 0 if x < 0


title: Shape detection for time series data

text: I have a large collection of time series - measurements taken every 15 minutes (96 measurements in a day) over the span of 1 year at various different locations. I've broken up each time series into 365 separate smaller time series, 1 for each day of the year. Looking at these time series, there are certainly many distinct shapes for a single day. Some look sinusoidal, some are constant, some look like a random stochastic process, some look parabolic, and some look like U's. What I would like to do is use an algorithm that can find these common shapes. I thought about clustering, and using the cluster centroids to define common shapes, but wanted to check with the community if this is right. So far, I've looked at Dynamic Time Warp as a metric, but it seems like that metric requires a lot of computation. I've also found http://mox.polimi.it/it/progetti/pubblicazioni/quaderni/13-2008.pdf from SE. I also saw Is it possible to do time-series clustering based on curve shape? but this question was from 2010 and might be outdated. Another idea I had was to take eigendecompositions of matrices that were formatted as: Matrix $M_i$ is a matrix of all time series observed on day $i$. Every row of matrix $M_i$ is a time series of length 96. Then, I would do 365 eigendecompositions, and use the eigenvectors as common shapes. Does this sound reasonable? Thanks!


title: Calculating the Kolmogorov-Smirnov coefficient

text: I need to calculate the coefficient of the Kolmogorov-Smirnov test for any given null hypothesis rejection level. For example, have a look at the table in this wikipedia entry. It only gives $c(\alpha)$ for $\alpha=0.10, 0.05, 0.025, 0.01, 0.005, 0.001$. How can I calculate $c(\alpha)$ for a wider range of $\alpha$s? If it is not trivial, can any one guide me on a place I can find a table for $c(\alpha)$ with larger $\alpha$ values? Thanks in advance.


title: Interpretation of lasso recovery results

text: When people say that lasso regression can under certain assumptions recover "the support", i.e. non-zero regression weights, what does this mean? This cannot mean _causal_ recovery, because Pearl has taught us you cannot infer causality from correlation (and lasso is just glorified correlation). So are these non-zeros in some way the minimum set of variables that cannot be explained by others? Because people tend to speak in terms of removing "spurious" and recovering "non-spurious" correlation, but really how can correlation be spurious, it's just the association between two random variables, there's no room for causal language here. I am confused.


title: How to evaluate a regression model's sensitivity to noise

text: How can I investigate the sensitivity of a regression model to noise? A basic idea is to add some (Gaussian) noise to the dependent and/or independent variables and (re)evaluate the RMSE. However, the problem with this approach is, that one is mostly interested in modeling the true data, and not the noisy data. Hence, is it reasonable to measure the RMSE on the original dependent variable, not the noisy one, even if the latter has been used for model building?


title: Categorical response variable prediction

text: I have the following kind of data (coded in R):               v.a = c('cat', 'dog', 'dog', 'goat', 'cat', 'goat', 'dog', 'dog')     v.b = c(1, 2, 1, 2, 1, 2, 1, 2)     v.c = c('blue', 'red', 'blue', 'red', 'red', 'blue', 'yellow', 'yellow')     set.seed(12)     v.d = rnorm(8)     aov(v.a ~ v.b + v.c + v.d) # Error      I would like to know if the value of `v.b` or the value of `v.c` has any ability to predict the value of `v.a`. I would run an ANOVA (as shown above) but I think it does not make any sense since my response variable is not ordinal (it is categorical). What should I do?


title: Proof that there is no autocorrelation

text: I was wondering whether my reasoning is correct and whether I've written it down correctly: Given the following equation $bread_i=\beta_{0} + \beta_{1}wage + u_{t}$, where $u_{t}$ is normally distributed with a mean of zero and a variancce of $\sigma^{2}$ and is independent of $wage$. Furthermore we know that there is no autocorrelation, $bread$ and $wage$ are not cointegrated, however there is a unit root present in $bread$ and $wage$. Given a first difference model $\Delta bread_{t} = \beta_1 \Delta wage_{t}+\epsilon_{t}$, where $\epsilon = \Delta u_{t}$. Is the error term $\epsilon_{t}$ in this equation uncorrelated over time? Please motivate your answer with a formal derivation. Well I thought the following: Since $\epsilon = u_{t} - u_{t-1}$ and we know that there is no autocorrelation between them, after first differencing there will be no autocorrelation afterwards. Question is of course how to show it? Using COV maybe? I am currently trying to derive it, but I'm a bit stuck. Thanks in advance!


title: Is my interpretation of Neural Network results correct?

text: I use a neural network with a topology of 17-30-1 (sigmoid, atan error function, mse as cost function, 5-fold cv) for text classification. (It's closely related to a previous question of mine.) The input data is quite noisy thus I could live with not a "perfect" classification score, but the results I get are probably too bad (or even just random) and thus I ask for your opinion.   1. The training error is around 0.06-0.09 (MSE), i.e. in average each classification differs approx. 0.25-0.3 from the predict label; in this binary case with a class threshold of 0.5 this might be acceptable. What do you think?   2. The test error (MSE) is unfortunately around 0.20 sometimes even 0.25; i.e. the effective error for a test sample is around 0.5, which to me means that the network a) suffers from high variance and b) is just as good as random guessing.  I don't need a perfect classification, but the network should however represent the patterns of the input data. But with this results I think the neural network is more or less useless or rather the input features are crap.


title: Interpretation of anova output with non significance

text: I have completed analysis on the effects of two drug treatments over a period of time on the CD4 cell count of a number of patients. I have taken the square root of the initial CD4 count as a covariate and I have taken a summary measure of the 'slopes' for each patient. My model is the following:               > slopes.aov <- aov(individual.slope.trans[-29] ~ sqrt(initialCD4)[-29] + treatment.fac[-29])          > summary(slopes.aov)                        Df Sum Sq  Mean Sq F value Pr(>F)     sqrt(initialCD4)[-29]   1 0.0060 0.006027   0.990  0.322     treatment.fac[-29]      1 0.0082 0.008184   1.344  0.249     Residuals             109 0.6638 0.006090               I am quite new to data analysis I am not quite sure how to interpret this model? Can I still use the regression coefficients to describe my summary measure even though we have no significance. I am really struggling with this. Also how can we describe to effect of the covariate.? I understand that I have no evidence to suggest that the treatments have an effect on the 'slopes' (my summary measure) and so I have no evidence to say that one treatment is performing better than another.


title: Distribution for number of trials before a fixed sum is reached

text: I'm trying to figure out the distribution over a number of trials until a stopping condition is met. In particular, imagine we are observing samples of a uniformly distributed random variable, $X \in \\{a,a+1,...,b\\}$. We continue to observe samples until the sum of these numbers passes some threshold $Z$ then we stop the process. That is, we have the set $\\{x_1, ..., x_T | \sum_{i=1}^T x_i \geq Z\\}$. I'm interested in the distribution over the random variable $T$ (the number of trials). This is a similar to a multinomial distribution, or a negative multinomial, but not the same. Assuming that $b << Z$, we can ignore the small rounding error right at the end assume that $\sum_{i=1}^T x_i = Z$. We can then make some basic observations, like $\frac{Z}{a} \leq T \leq \frac{Z}{b}$. Also, from experiments it seems that $\mathbb{E}[T] = \frac{Z}{\mathbb{E}[X]}$, not sure about the variance though, or the general shape of the distribution. edit: This question can be answered in part by observing the following, first denote $Y = \sum_{i=1}^N X_i$, we of course, have $Y/N \overset{d}{\to} \mathbb{E}[X]$. Naturally, for any large $N$ we expect $Y \approx N\mathbb{E}[X]$. In fact, this could also seen from Hoeffding's inequality. So, for a large enough $Z$, if we constrain $Y=Z$, we have $Z \approx T\mathbb{E}[X]$ or $T \approx \frac{Z}{\mathbb{E}[X]}$ as I observed.


title: Can survival analysis be used to predict earthquakes?

text: Given survival analysis relates to an analysis/prediction of time to an event, I was wondering if it was possible to be used to predict eathquakes. If so, how would one go about carrying out that analysis on say one particular region say Japan or California? So I guess a question related to that say...for example if we took a univariate time-series based dataset looking at the tremors over time in Japan, could this be converted in a manner that would applicable to be useable for survival analysis... this is assuming that potentially prior small tremors might be indicative of a much larger earthquake to come.


title: In a hidden Markov model, how do all observations and one state give you all states?

text: $Y^n$ are the observations of our HMM, where $Y_i=a_i$ is a single observation, where $a_i \epsilon \\{0,1\\}$. For example, $Y^n = k^n$ where $k^n=\\{0,1,1,0\\}$ $X^n$ are the actual states of our HMM, of which we know some (any) $X_i$. How can I both 1.intuitively and 2.carefully (using actual math) reason about the following fact: If we have $Y^n = k^n$ where $k^n=\\{0,1,1,0\\}$ and **any** $X_i$ then we know the entire $X^n$. I think it has something to do with the Markov implication that > Conditioned on the present, the past and the future are independent but I'm not very clear on how to proceed.


title: Different random effects in nlme and nlmer

text: I estimated a nonlinear exponential $f(t)=\alpha-(\alpha-\beta)*\exp(-\gamma*t)$ random effects model in `R`. I estimated the same model form using the identical data set but with `nlme()` from the `nlme` library and `nlmer()` from the `lme4` library. The fixed effect estimates are quite similar comparing `nlme` with `nlmer`: **nlme:**               Fixed effects: alpha + beta + gamma ~ 1                Value Std.Error   DF  t-value p-value     alpha 29.508227 0.4166492 3498 70.82272       0     beta   5.082116 0.1584696 3498 32.06998       0     gamma  0.313008 0.0070209 3498 44.58244       0      **nlmer:**               Fixed effects:           Estimate Std. Error t value     alpha  30.0432     0.5432   55.31     beta    4.9541     0.2110   23.47     gamma   0.2551     0.0135   18.90      However, the random effects estimates along with their correlation matrix are different. This is a consistent finding. It also holds for other nonlinear functions and other data. See e.g. the correlation between alpha and gamma. It is also striking how large the variances in `nlmer` are -- they actually appear to be far too large **nlme:**               > VarCorr(model)     Subject = pdLogChol(list(alpha ~ 1,beta ~ 1,gamma ~ 1))               Variance    StdDev   Corr              alpha    42.19955753 6.496119 alpha  beta       beta      4.84969745 2.202203  0.494            gamma     0.01055961 0.102760 -0.643  0.202     Residual  1.94836594 1.395839       **nlmer:**               Random effects:      Groups   Name  Variance  Std.Dev. Corr                Subject  alpha 72.551249 8.51770                               beta   9.840010 3.13688   0.458                       gamma  0.043274 0.20802  -0.195  0.142       Residual        1.716048 1.30998                     Number of obs: 3750, groups: Subject, 250      What is the problem here? Why are the variances in `nlmer` almost twice as large as those given in the `nlme` solution?


title: Treatment of Ordinal Predicting Variable

text: I am trying to perform an regression analysis where the response variable is ordinal and the 15 out of the 16 predicting variables are also ordinal. Besides treating all these ordinal predictors as factors, is there any other option? What is the best way to treat these ordinal predictors? Thank you.


title: Structure of data and function call for recurrent event data with time-dependent variables

text: I'm attempting to estimate the effect of 2 drugs (`drug1`, `drug2`) on the likelihood of a patient falling (`event`). The patients can fall more than once and can be put on or taken off of the the drugs at any point. My question is how the data should be structured with regard to the time period (days), specifically whether there needs to be overlap between the days. There are two reasons why I think my structure is wrong, the first being a seemingly incorrect `N`. I am also getting some errors where the time period is a single day (i.e. `time1=4`, `time2=4`) and am unsure how these should be coded. Should the start time of subsequent entries be the stop time of the previous entry? I've tried it both ways (with and without overlap), and while having overlap gets rid of the warning, the `N` is still incorrect.               Warning message:     In Surv(time = c(0, 2, 7, 15, 20, 0, 18, 27, 32, 35, 39, 46, 53,  :       Stop time must be > start time, NA created      Right now I have the data set up where the beginning of the next entry is the next day. Unique patients are identified by their `chart numbers`.               Time1    Time2    Drug1    Drug2   Event    ChartNo         0        2        1        0       0        123         3       10        1        1       1        123        11       14        1        1       1        123         0       11        0        1       0        345         0       19        1        0       1        678         0        4        0        1       0        900         5       18        1        1       0        900      Patient 123 was on drug1 at the start to day 2, after which point they had drug2 added. They went from day 3 to day 10 on both drugs before falling the first time, then fell a second time on day 14 while still on both drugs. Patient 345 went 11 days on drug2 without falling (then was censored), etc. The actual estimation looks like this:               S <- Srv(time=time1, time2=time2, event=event)     cox.rms <- cph(S ~ Drug1 + Drug2 + cluster(ChartNo), surv=T)      My main concern is that the `n` for my analysis is reported to be `2017` (the number of rows in the data), when in actuality I only have `314` unique patients. I am unsure if this is normal or the result of some error I've made along the way.               > cox.rms$n     Status     No Event    Event          1884      133       The same is true when using `coxph()` from the survival package.                n= 2017, number of events= 133      The number of events is correct however. This Post seems to have it set up with the 'overlap' I described, but I am unsure about the `N`, and they don't seem to be clustering by `ID`.


title: Good way to plot a variable between conditions with multiple measurements

text: I've a frame containing > Participant----Condition----Duration   >  Participant1,One,2000   >  Participant1,One,2780   >  Participant1,One,200   >  Participant2,Two,2000   >  Participant2,Two,2340   >  Participant2,Two,2120   >  Participant2,Two,2420   >  Participant2,Two,2990   >  Participant3,One,2120 and so on... basically for each participant there are multiple measures of the duration. In order to compare the duration vs the condition i'm using > lme(Duration~Condition,random=~1|Participant) and > plot.lme(model,Condition~fitted(.)) to have an idea of the model. I'm wondering if there is a way to plot the real data, in a summary style like plotmeans. I mean I can use plotmeans with the relation is 1:1 between participant and measure, otherwise the mean is "pulled" by the participnat with more entries. Is there something similar for repeated measures? Grand mean maybe? Thanks


title: Comparison test in unit root

text: I am not clear how to compare all the tests in Unit Root. Which one is the most powerful among:   * ADF   * PP   * KPSS   * Elliott, Rothenberg & Stock Unit Root Test


title: Logistic regression: maximizing true positives - false positives

text: I have a logistic regression model (fit via glmnet in R with elastic net regularization), and I would like to maximize the difference between true positives and false positives. In order to do this, the following procedure came to mind:   1. Fit standard logistic regression model   2. Using prediction threshold as 0.5, identify all positive predictions   3. Assign weight 1 for positively predicted observations, 0 for all others   4. Fit weighted logistic regression model What would be the flaws with this approach? What would be the correct way to proceed with this problem? The reason for wanting to maximize the difference between the number of true positives and false negatives is due to the design of my application. As part of a class project, I am building a autonomous participant in an online marketplace - if my model predicts it can buy something and sell it later at a higher price, it places a bid. I would like to stick to logistic regression and output binary outcomes (win, lose) based on fixed costs and unit price increments (I gain or lose the same amount on every transaction). A false positive hurts me because it means that I buy something and am unable to sell it for a higher price. However, a false negative doesn't hurt me (only in terms of opportunity cost) because it just means if I didn't buy, but if I had, I would have made money. Similarly, a true positive benefits me because I buy and then sell for a higher price, but a true negative doesn't benefit me because I didn't take any action. I agree that the 0.5 cut-off is completely arbitrary, and when I optimized the model from step 1 on the prediction threshold which yields the highest difference between true/false positives, it turns out to be closer to 0.4. I think this is due to the skewed nature of my data - the ratio between negatives and positives is about 1:3. Right now, I am following the following steps:   1. Split data intto training/test   2. Fit model on training, make predictions in test set and compute difference between true/false positives   3. Fit model on full, make predictions in test set and compute difference between true/false positives The difference between true/false positives is smaller in step #3 than in step #2, despite the training set being a subset of the full set. Since I don't care whether the model in #3 has more true negatives and less false negatives, is there anything I can do without altering the likelihood function itself?


title: Comparing negative-binomial and poisson standard errors

text: I have this question to answer for a class and am completely stumped. "The GLM log link and a dummy variable for gender (1=male, 0=female) has gender estimate 0.308. The SE is 0.038 assuming a Poisson distribution and 0.127 assuming a negative binomial model. Why are the SE values so different?" Any thoughts would be much appreciated.


title: How to prove that $X^T$e = 0

text: Hi I need to prove that $X^T$e = 0 where e is the residual in multiple linear regression model in matrix algebra? Need some guidance on how to do it... Is there any good pdf for the proofs for multiple linear regression model for matrix alegbra?


title: Is it always bad to retrain your model to include predicted data?

text: I understand intuitively why this is a horrible idea - you assume your model is correct and then increase your number of observations which will likely result in a poor fit on future data. I'm wondering if there is some mathematical/statistical property to describe this, or if there is any rare case where this may not be as fatal as I am thinking?


title: Comparing frequency distributions using Mann-Whitney test

text: I have two frequency distributions; which each are frequency of an event. I would like to compare these distributions, and conclude that these events and conclude if these events are statistically 1 event or they are two distinct events. however, it seems the size of these distributions are playing an important role in the final conclusion. meaning when I compare these events by looking at 1000 data point, the p-value is statistically more significant; meaning two distributions are two distinct distributions. However, when I narrow it down to 100 data point my test shows less significant and I can conclude these 2 distributions are sampled from 1 distribution. Here are my questions: 1) How can I correct my p-value to be robust given the sample size? 2) Is Mann-Whitney a right statistical test for this case? If not, can you point me to the right test?


title: Variance of marginal posterior distribution

text: Suppose $Y_1,\dots,Y_n\mid\mu,\sigma^2 \sim \text{ iid } N(\mu,\sigma^2)$ and suppose the priors $\mu \mid \sigma^2 \sim N(\mu_0, \sigma^2 / \kappa_0)$ and $1/\sigma^2 \sim \text{gamma}(\nu_0/2, \nu_0 \sigma_0^2 / 2)$ are placed on the unknown parameters. Then \begin{align*} p(\sigma^2 \mid y_1,\dots,y_n) &\propto p(\sigma^2)p(y_1,\dots,y_n\mid\sigma^2) \\\ &= p(\sigma^2) \int p(y_1,\dots,y_n \mid \mu,\sigma^2)p(\mu\mid\sigma^2)d\mu \end{align*} It turns out (see Hoff) with these priors that the marginal posterior distribution on $\sigma^2$ is $$1/\sigma^2\mid y_1,\dots,y_n \sim \text{gamma}(\nu_n/2,\nu_n\sigma_n^2/2)$$ where $$\nu_n = \nu_0 + n, \text{ and } \sigma_n^2 = \frac{1}{\nu_n}\left( \nu_0\sigma_0^2 + (n-1)s^2 + \frac{\kappa_0 n}{\kappa_0 + n} ( \bar{y}-\mu_0 )^2 \right)$$ where $\bar{y}$ and $s^2$ are the sample mean and unbiased sample variance. Thus we have that $$ \text{Var}\left[\frac{1}{\sigma^2} \mid y_1,\dots,y_n\right] = (\nu_n/2) / (\nu_n\sigma_n^2/2)^2 = \frac{2}{\nu_n\sigma_n^4}. $$ Now suppose we fix all hyperparameters except $\kappa_0$ and analyze the effect of increasing the variance of $\mu\mid\sigma^2$ (by decreasing $\kappa_0$). Since $f(\kappa_0)=\kappa_0 n / (\kappa_0 + n)$ is monotonically increasing as a function of $\kappa_0$, we see that \begin{align*} \kappa_{0,a} < \kappa_{0,b} &\Rightarrow \kappa_{0,a} n / (\kappa_{0,a} + n) < \kappa_{0,b} n / (\kappa_{0,b} + n) \\\ &\Rightarrow \sigma_{n,a}^2 < \sigma_{n,b}^2 \\\ &\Rightarrow \sigma_{n,a}^4 < \sigma_{n,b}^4 \\\ &\Rightarrow \frac{2}{\nu_n\sigma_{n,b}^4} < \frac{2}{\nu_n\sigma_{n,a}^4}, \end{align*} i.e., increasing the variance of $\mu\mid\sigma^2$ results in an increase in the marginal posterior variance $\text{Var}\left[1/\sigma^2\mid y_1,\dots,y_n\right]$. Intuitively, this result makes a lot of sense because if $p(\mu \mid \sigma^2)$ is more spread out, the averaging (integral) over $p(\mu \mid \sigma^2)$ above makes $p(\sigma^2 \mid y_1,\dots,y_n)$ more spread out. This makes me think that this result should be able to be proved in general, i.e., given a sampling model $f_{Y\mid\theta,\phi}$ depending on two unknown parameters $\theta$ and $\phi$, and priors $f_{\theta\mid\phi}$ and $f_{\phi}$, an increase in $\text{Var}\left[ \theta\mid\phi\right]$ should result in an increase in $\text{Var}\left[ \phi \mid y_1,\dots,y_n \right]$. Does anyone know of a general result like this?


title: Applied Analysis Question

text: I'm a newbie analyst and I'm facing a machine learning / regression problem which I cannot solve. The data I need to use in my analysis consists of information about press subscriptions of some companies. Each company has its unique ID but it further splits into separate company's "employees" with different values of response variable. I provide an example dataset: ![enter image description here](http://i.stack.imgur.com/VdDUA.png) Sometimes the situation goes into extremes - there are 20 observations of the same ID (company) which differ in some values but most of them remain the same (say size of company). My thoughts are when I use any statistical method repeatability of categorical values for a given company causes huge bias. I presume there must be some way of "weighing" of categorical observations according to how many times they appear for a given ID i.e. if there 20 employees under a given each the importance of each observation (for a given variable) is divided by 20 so the prevalence of some pattern in company which is very big does not bias my prediction. I hope I've managed to rely what I have trouble with. I'd be extremely grateful for references to any tools/methods/articles or your experiences which may be helpful in handling this matter.


title: How to split a data set to do 10-fold cross validation

text: Now I have a `R` data frame (training), can anyone tell me how to randomly split this data set to do 10-fold cross validation?


title: How are local and global plug-in bandwidths different in kernel smoothing regression?

text: I'm looking into an R package 'lokern.' It provides two bandwidths selectors, global and local plug-in bandwidth. I would like to understand the difference between two methods. My understanding is 'global' bandwidth uses one bandwidth for the entire data set, and 'local' bandwidth uses different bandwidths in one data set according to some criterion (such as density of the data etc). Is that correct? Is there any criterion when local bandwidth is more useful and when global is more useful?


title: Are the angles between iid training data and the classifier learned from them independent?

text: Suppose I have a training set $\mathcal{T}^N = \\{(x_i , y_i) \in \mathbb{R}^{d} \times \\{0,1\\}\\}$, where the data are i.i.d draws from some distribution. Now I learn a classifier $\hat{h}$ from these data using some deterministic algorithm. WLOG take the Euclidean norms of the $x_i$ and of $\hat{h}$ to be $1$ and denote by $\theta_i := \arccos(x_i^T\hat{h})$ the angle between $x_i$ and $\hat{h}$. Obviously it is true that for _any_ fixed $\hat{h}$ which is not a function of the data the angles $\theta_i$ are independent random variables, but here $\hat{h}$ is a function of the iid $x_i$: Is it still true that the angles $\theta_i$ are independent?


title: Correct structuring of random effects?

text: I have produced a mixed effects model as follows: `lmer(TotalPayoff~Type+Game+PgvnD*Asym+(1|Subject)+(1|Pairing),REML=FALSE,data=table)` each pairing contains 2 subjects and each Subject is in 2 different pairings with different subjects. Have I structured the random effects in the correct way to control for the dependence of `TotalPayoff` between data from the same subject as well as between data from the same pairing? (If one subject in a pairing has a high `TotalPayoff` the other will too).


title: Calulating a mean coefficient effect based on quantile regression estimates

text: I have used quantile regression to estimate a particular coefficient for the 10th, 20th, ... 90th percentiles. Now, I want to estimate the mean effect of the coefficient across the whole distribution. Using ordinary least squares to estimate the mean effect, however, is inappropriate given the heteroscedastic, non-Gaussian distribution of the model errors. Would it be appropriate to use the mean of the quantile regression estimates as an approximation of the coefficient's mean effect?


title: Finding uncertainty in coefficients from polyfit in Matlab (no Toolboxes)

text: I have a data set of x- and y-values, that I want make a linear fit of. Using `polyfit(x,y,1)` I get the coefficients `a` and `b` for a linear fit `ax = b` for this data, but I would also like to find the uncertainty or standard deviation for these coefficients. Does anyone know an easy way of doing this? My Google-fu only gave me this result, and seeing as the last answer in that thread is a correction to first answer, I don't know if I should trust any of the answers. The answer from that thread is:               [z,s]=polyfit(x,y,1);     ste = sqrt(diag(inv(s.R)*inv(s.R')).*s.normr.^2./s.df);      Here is the data set if anyone needs it:               x = [1.3979 1.6990 2.0000 2.3010 2.6021 2.9031];     y = [-1.5219 -1.6300 -1.8249 -2.2584 -2.5054 -2.6248];      (I realize it looks more like a third order polynomial, but they should be linear according to the theory behind my measurements.)


title: Is this a case control study or a retrospective cohort study?

text: From the echocardiography database at a tertiary care center, 1,625 subjects (mean age, 44 6 14 years; 47% men) with normal echocardiographic findings between 2000 and 2009 were identified. Gender dif- ferences and association with body surface area were assessed retrospectively for right atrial long- axis and short-axis dimensions, right ventricular short-axis dimension, end- diastolic and end-systolic right ventricular area, right ventricular fractional area change, and tricuspid annular plane systolic excursion. The impact of nor- mal values stratified for gender and body surface area was tested in 24 patients with moderate-sized to large atrial septal defects.


title: VIF values and interactions in multiple regression

text: I am running a multiple regression of the form y~a+b+c+ab+ac+bc I have checked the VIF values for the direct effects - should I check them for the interactions? I am assuming not as that would equate to looking at the multicolinearity between a variable and itself (albeit in an interaction) which will surely be very high?


title: Point-biserial correlation with unequal groups in dichotomous variable

text: The dichotomous variable I want to use for a point-biserial correlation has very unequal n's ($n_1=140$, $n_2=6$). Would this be breaking any rules? What would this mean for my results? Would it just limit the statistical power?


title: Linear model estimation in the presence of heteroscedasticity

text: Assuming a sample of random variables where the error terms for each random variable ($y_{i}$) are given by $\epsilon_{1}, \dots, \epsilon_{n} \sim N(0, \sigma^{2})$, a linear model is developed such that: $$ y_{i} = \beta_{0} + \epsilon_{i} $$ then, using OLS, finding, $$ \frac{d}{d\hat{\beta}_{0}} \sum_{i}^{n}(y_{i} - \hat{\beta}_{0})^{2} = 0 $$ and solving for $\hat{\beta}_{0}$ yields, $\hat{\beta}_{0} = \frac{1}{n}\sum_{i}^{n} y_{i} = \bar{y}$ However, now, let's assume that there are _two_ random variables where the error terms for the corresponding random variable ($y_{i}$) are given by $\epsilon_{1} \sim N(0, \sigma^{2})$ and $\epsilon_{2} \sim N(0, n\sigma^{2})$. Once again, the a linear model is developed such that: $$ y_{i} = \beta_{0} + \epsilon_{i} $$ Can the estimator for $\beta_{0}$ given above be used for this model? If not, in light of the heteroscedasticity involved in the error terms, how is an estimator for $\beta_{0}$ developed? I've seen linear models that, given a data point, use indicator functions to determine which terms in the model are non-zero. Perhaps this applies here?


title: Concordance vs. Sensitivity

text: I am confused between these two terms: _sensitivity_ and _concordance_. What I understand about these two terms: **Concordance:** the number of pairs where actual 1s have higher predicted probability of becoming 1s than 0s divided by total number of pairs. **Sensitivity:** percentage of 1s correctly predicted as 1. By their definitions, both seems to be the same. I have read somewhere that concordance = (true positive + true negative) / sample size. I have tried to calculate concordance using this formula. But it seems to be incorrect. First, I have run logistic regression in SAS using PROC LOGISTIC with the CTABLE PPROB=0.5 options. Then I have calculated true positive pairs from (sensitivity percentage $\times$ number of 1s) and true negative pairs from (specificity percentage $\times$ number of 0s). I added them up and divided by total responses. This is not matching up with concordance score.


title: Measuring the quality of an opinion mining algorithm

text: I am working on a piece of software to detect opinions in text. As a simplistic example, I would like my algorithm to tell me that Andrea has a positive opinion **`+1`** (rather than neutral **`0`** or negative **`-1`** ) of New York when I enter the following sentence: > **Andrea loves New York.** Now I have hundreds of millions of unstructured text files with sentences just like this one, or way more complex, so it is not an option to analyse these sentences by hand! Hence, the piece of software. But I have different methods, based on different papers, or on my own experience with opinion mining.                            method_a    method_b    method_c    labelled_by_hand(n)     positive     20%         30%         50%         30%     neutral      80%         30%          0%         20%     negative      0%         40%         50%         50%          sample size n = 100      I would like to know how to **assess the quality of a given algorithm**. I have no distribution, no statistical model that I can assume to be true so what I read here or there about AIC and BIC and other methods just seem way too elaborated (I may be wrong, but I don't have any _common variables_ to use to compare my models to each other). So my questions:   1. how many samples should I label by hand to achieve a given quality for my assessment test?   2. which measures are commonly used for a simple comparative task like this? Precision and recall?


title: How to perform T-Test in R including hypothesis

text: I would like to do a t-test in R. Unfortunately, statistics has never been my best subject, so I'm wondering what the output of t.test in `R` exactly mean.   For testing purpose I performed a t-test on two normal distributed vectors with equal parameters. Since both distributions have the same parameters, the difference of the variance is 0:               t.test(rnorm(1e7,mean=0,sd=1),rnorm(1e7,mean=0,sd=1),mu=0)          Welch Two Sample t-test          data:  rnorm(1e+07, mean = 0, sd = 1) and rnorm(1e+07, mean = 0, sd = 1)      t = 0.5014, df = 2e+07, p-value = 0.6161     alternative hypothesis: true difference in means is not equal to 0      95 percent confidence interval:      -0.000652202  0.001100571      sample estimates:     mean of x    mean of y      2.836860e-04 5.950136e-05       This obviously gives me the t of the t-test. In order to do something useful with it, I have to formulate a hypothesis which should be tested.   1. How / where (in terms of `R`) do I formulate a hypothesis?   2. How do I know if my hypothesis is significant? My idea was: An interesting hypothesis on two data vectors might be: $$ H_{0}: \mu_{X}-\mu_{Y} = w $$ vs. $$ H_{1}: \mu_{X}-\mu_{Y} \neq w $$ Where $\mu_{X}, \mu_{Y}$ are mean values of the underlying deviation of $X,Y$.   Since these might be different from the mean of the given data vectors, I have to use `t` in order to accept or refuse $H_{0}$ ($H_{1}$) (so far correct?). Looking at wikipedia, $H_{0}$ is refused if $|t| > t(1-\frac{1}{2}\alpha, m+n-2)$ where $\alpha = 0.05$ (standard paramter of `t.test` in `R`) and $df=m+n-2$. In this case t(..,..) is some giant table where I can look up the corresponding value.   1. Where do I get this value in `R` and at which point did is specify $w$ for $H_{0}, H_{1}$?


title: Violated assumptions for regression

text: I am not sure what to do with this regression residuals. ![Regression Diagnostics](http://i.stack.imgur.com/hXxz4.jpg) The Normal Q-Q plot is indicating that residuals are non-normal (not too extreme) and residual variance is not constant (not too bad). My input variables are factor variables so I am not able to normalize it. What do you guys think about this plot? Should I continue with this plot and say that assumptions are not violated and continue? In fact, I think assumptions are violated slightly but not extremely. What should I do if I think assumptions are violated?


title: Featuring Engineering from Trends in the Training Set

text: I have a predictive model with just OK performance, and I'm trying to improve it with feature engineering. My question: is it valid to create new features by looking at trends in the training set? For example, I noticed that the metric I am trying to predict tends to be higher during the months of December, January, and February. Would it be reasonable to include a feature that is '1' for those months, and '0' for the other months? Thanks!


title: Haar prior for von Mises distribution

text: Ok, Let me tell you that this is the very first time that I have no idea with the question below. I can not find a solution or anything that will lead me to it. I say this to prevent comments "what did you do?". I apologize in advance. This is one of the questions in Peter M. Lee's Bayesian Statistics: An Introduction book, it is Exercises section of Chapter 3. How can I find a Haar prior for a location parameter on the circle (such as $\mu$ in the case of von Mises's distribution).


title: Newey-West t-stats and critical values

text: When using Newey-West standard errors for my t-statistics of the slope coefficients in an OLS regression, can I still use the usual critical values for two-sided tests? (1.645 for 10% significance, 1.96 for 5% and 2.58 for 1%)


title: How to select next item in polytomous CAT / MIRT?

text: To make a Computerized Adaptive Test out of a sample of 20 dichotomous items a typical course of action would be to: 1) calibrate respondent data with a R package like mirt or ltm using Rasch, 2PL etc. 2) create an itembank based on obtained item parameters using catR then, for each item/test taken: 3) use catR's nextItem function and method (e.g. MFI) for selecting a next item based on previous items and answers 4) choose a stopping rule, e.g. stop if SEM < 0.2. 5) use last theta and SEM as 'test results' However, regarding polytomous items and multidimensional tests I'm a bit confused. For a polytomous, multidimensional CAT, e.g. 100 items (5 scales) of a personality test with a 5-point-likert scale, the test logic, item parameters (GRM/GPCM) and item selection are very different. Specific question: is there a polytomous version of catR's nextItem function available? Broader question: What methods/packages/steps would you recommend for making a working polytomous CAT with MIRT?


my question is: What is the difference between the variance of an individual forecast and the variance of a mean forecast? I know that the two variances are different, but what consequences hat the forecast on the data that is forecasted?


title: Which distribution should I use to describe this zero-enriched data?

text: I want to use a distribution to predict the p-value (whether it's just noise) of each element in this dataset. But I don't know what distribution should I use. Could anyone give me some tips about it? ![The histogram of the data](http://i.stack.imgur.com/96aNg.png)


title: What are $\ell_p$ norms and how are they relevant to regularization?

text: I have been seeing a lot of papers on sparse representations lately, and most of them use the $\ell_p$ norm and do some minimization. My question is, what is the $\ell_p$ norm, and the $\ell_{p, q}$ mixed norm? And how are they relevant to regularization? Thanks


title: Understanding basic identifiability

text: I'm having trouble understanding identifiability. Specifically, I'm not sure, in the following example, why $P\left(C\right)$ cannot be identified. Here's the example: You have 2 unfair coins with the following probabilities: $$P\left(H|C_1\right)=1/4$$ $$P\left(H|C_2\right)=3/4$$ Each coin has an equal probability of being chosen, so $P\left(C=C_1\right)=P\left(C=C_2\right)=1/2$. My question is, given that I only observe $H$, why can't I identify $P\left(C\right)$ using basic algebra. That is, \begin{eqnarray} P\left(H\right)&=&P(C_1)P\left(H|C_1\right) + P\left(C_2\right)P\left(H|C_2\right)\\\ &=& P\left(C\right)\left(P\left(H|C_1\right)+P\left(H|C_2\right)\right)\\\ P\left(C\right) &=& \frac{P\left(H\right)}{P\left(H|C_1\right)+P\left(H|C_2\right)} \end{eqnarray} In this case $P\left(C\right)$ would seem to be identifiable, so I'm not sure where I'm going wrong. Note that I'm not looking to identify $C$ itself (the actual coin), I want to know the distribution on $C$.


title: Time series models for projected sales

text: Suppose we are given various reports on projected sales of Product $X$ for $5$ years starting in $2013$. All of the reports just give the results and not the methodology as to how they arrived at the predictions. Suppose we are interested in the projected sales of Product $X$ for $10$ years. So we already have the projected prices from $2013-2018$ and are interested in the projected salesfrom $2019-2023$. Can we use a time series model for these years?


title: Non-normal distribution even with Kruskal-Wallis test

text: Thanks to anyone who can help! I`ve done a survey with Likert-scale questions with about 1,500 participants. Initial analysis in SPSS showed no outliers (using box-plot) but scores were not normally distributed for any of the factors, as assessed by Shapiro-Wilk's test (p < .001). Since the independent variable had 5 groups, I ran a Kruskal-Wallis H test, finding significant results for the tested variables (p < .001 for most), but an inspection of the boxplots that accompany these new results also showed that he shapes of the distributions are not similar (i.e. the medians were sometimes higher or lower than other variables being compared etc). My questions: Is this a serious flaw? Can I use these results? If so, should I mention that in at least one independent variable, the groups displayed non-normal distributions? thanks!


title: Inconsistent results from attempting to count matching values in lists within data frames in R

text: Pretty new to R and to stackexchange. My apologies if I make any faux pas. I'm working with the North Carolina births data from irss.unc.edu. I have a data frame `NCBIRTH1450` with several lists, all of which have numeric or NA values. I'm just trying to calculate some proportions from the data. For one of the lists, I'm able to get the proportion of male babies by using               pBoyBaby <- sum(SEX==1)/sum(!is.na(SEX))        The above returns a numeric value for `pBoyBaby`. My aim is to take the number of matching values and divide by the number of non-missing values, and I get an appropriate result. However, when I attempt the same function on another column of the same type,               pMomSmokes <- sum(SMOKE==1)/sum(!is.na(SMOKE))        The object `pMomSmokes` takes the value `NA`. How could I identify the difference between the natures of `SMOKE` and `SEX`, for instance, that explain this discrepancy? I really don't know why it's doing this. In case it's relevant, `SEX` has values of either `1` or `2` \-- with no `NA`'s -- while `SMOKE` has the values `0` or `1`, with a few `NA`'s. Thanks!


title: Can SVM do stream learning one example at a time?

text: I have a streaming data set, examples are available one at a time. I would need to do multi class classification on them. As soon as I fed a training example to the learning process, I have to discard the example. Concurrently, I'm also using the latest model to perform prediction on unlabelled data. As far as I know, a neural network is able to do stream learning by feeding examples one at a time and performing forward propagation and backward propogation on the example. Can a SVM perform stream learning one example at a time and discard the example straight away?


title: Generating a Gaussian dataset in MATLAB

text: I want to generate a bivariate Gaussian dataset. The dataset includes a total of 800 results drawn randomly from four two-dimensional Gaussian classes with means $(-3,0)'$, $(0,0)'$, $(3,0)'$, and $(6,0)'$, all with the same variance- covariance matrix $$\Sigma = \pmatrix{0.5 & 0.05 \\\ 0.05 &0.5}.$$ How can I do that in MATLAB? I'm not expert in MATLAB.


title: Computing probability of proportions of (possibly) overlapping sets

text: Consider $n$ _fixed_ subsets of a finite universe $U$: $A,B,C,D,\ldots$ Then let $Coll$ be the collection of sets we get from the closure of the set $\\{U,A,B,C,D, \ldots\\}$ under the operations $\cup,\cap$ and set-theoretic difference. (As an aside, to help picture Coll, note that $|Coll|≤ 2^{2^n}$ Imagine an n-ary Venn diagram partitioning U into 2ⁿ disjoint pieces--there are then 2^(2ⁿ) combinations of the pieces, which will be distinct if the 2ⁿ pieces are all nonempty.) I am looking for methods to estimate the size of a given element of $Coll$ given _some_ of the sizes of other sets in $Coll$. Specifically, I would like to know the probability distribution of the cardinalities of a given set in $Coll$, given _some_ precise values for the sizes of some of the other sets in $Coll$. I don't know what probability model is appropriate to use is no further information is known--what methods should I look to for this?


title: What's the recommended method for charting low values?

text: ![enter image description here](http://i.stack.imgur.com/AwHU7.png) I need to chart a number of indicators to evaluate emergency response times for a fire department. In our case, most values should be around 2 minutes, but the range can be between 0 and 10 minutes. I've been working with bar charts (and Few's bullet charts), however, since bar charts work on comparing lengths, and our target values are shorter (around 20% of the bar's potential total length), I'm curious if there's a better way to chart these values. I've considered charting the inverse value (bar is 8 minutes, instead of 2), but that seems counter-intuitive to read. Any help is appreciated. I'm trying to show the range in turnout time for each unit (blue bar, min to max) overlayed on excellent, mediocre and poor ranges (gray backgrounds). I'm working on adding a marker for the 90th percentile (our standard, but not on chart yet). My problem is as performance improves, the times reduce and consequently, the less length of the bar is used for comparison. I'm trying to determine if there's a better way to mark their performance indicators, since a lower number is better (as opposed to something like profit, where bigger is better). Sorry, it won't let me post an image to show you, but consider a basic bar chart with the x-axis as time and y-axis is categorical. Under each bar are indicators for 3 performance levels and a single marker is applied to each bar indicating the 90th percentile. Hope that helps. I was able to upload a picture of the chart. As values improve, they'll get smaller and move closer to the left edge of the bars, and the bars can shrink if the range decreases. I just want to make sure there's as much visual info as possible to make the data points legible.


title: How to correctly calculate linear model variable effects?

text: I'd like to estimate the unique effects for each variable in a linear model, however I am unsure if I calculating these correctly. I am using a model estimated using each of the variables in my data plus a binary interaction with each variable, e.g. `y ~ x + x:z` where x is a vector and z is two factors which is present in x. Doing this in R I believe a two levels factor is simply c(0,1). From what understand I can calculate the individual effects of the variables for each case of the binary interaction term is as follows:   * When the interaction term is 0, the variable effect of x on y is the estimated model coefficient for x   * When the interaction term is 1, the variable effect of x on y is the estimated model coefficient for x plus the coefficient for the interaction term of x:z Assuming this is correct I would also like to calculate the standard errors associated with the each effect. I looked at this question which states that the associated standard errors can be calculated as follows `sqrt(x + x:z + 2*cov(x,x:z))`. This however produces much smaller standard errors compared to the non-interaction term effects. This makes me assume I am doing something wrong. Could you tell me the correct way to estimate these effects?


title: What sort of data would be visualized with these colormaps?

text: I have a good feel for the types of data that go with most colormaps. But I don't know what kind of data one would want to use the following types of colormaps: Flag: ![Flag in matplotlib](http://i.stack.imgur.com/T3uMo.jpg) Prism: ![Prism in matplotlib](http://i.stack.imgur.com/yZ0zA.jpg) Source: http://dept.astro.lsa.umich.edu/~msshin/science/code/matplotlib_cm/ I can't see any use for them except as a fill pattern. What kinds of data benefit from these types of colormaps? And does this type of map have a specific name?


title: Effect of uncertainty in a parameter of a probability model on decisions

text: I am exposing several experimental organisms to a toxin in the food. If the dosage is high enough, they will die. I am examining the number that are dead after a day. For practical reasons, the dosage they get is not directly under my control. I have decided to model the number that are dead after one day as a binomial distribution, Bin(n,p). In order to figure out what p is for my model, I did an experiment with 300 organisms and counted the number that had died after a day. In the experiment of interest, I have exposed them to a potential enhancer of the the toxic effects. I would like to know if substantially more of the organisms have died. For a particular result m, I can easily compute the probability of m or greater. If this is less than 1% then I accept that the putative enhancer works otherwise I reject the hypothesis that the putative enhancer is an actual enhancer of the toxic effects. I am wondering if my analysis of the putative enhancer needs to take into account the uncertainty in the value of p. If so, how would I do it? (I am good with R and able to program if necessary.) Update: I just want to emphasize that I only care if my enhancer works i.e. more organisms die than I'm expecting. Any answer should leave out cases where fewer than expected organisms die.


title: Why does one get a 0.000 Mann-Whitney U value (not p-value) etc.?

text: 1) Is it okay to run a Mann-Whitney U test on unequal sample sizes i.e. two groups (N=11, N=12)or is there an alternative test? 2) What does it mean when one obtains a zero value for a Mann-Whitney U test (bearing in mind I obtained zero values before with equal sample sizes)               Ranks                    **Site    N   Mean Rank Sum of Ranks**     1         12    18.50     222.00     3         12     6.50      78.00     **Total   24**                          **Test Statistics**                 Mann-Whitney U     .000      Wilcoxon W       78.000      Z  -4.157      Asymp. Sig. (2-tailed) .000      Exact Sig. [2*(1-tailed Sig.)] .000b      Exact Sig. (2-tailed)  .000      Exact Sig. (1-tailed)  .000      Point Probability  .000           a Grouping Variable: Site      b Not corrected for ties.


title: Is it possible to make quantitative comparisons using income inequality metrics?

text: I've calculated income inequality metrics (Gini/Hoover/Theil coefficients) for several populations. I know I can make claims like "this population has a higher Gini coefficient, so its distribution in this variable is more unequal than this other population." But is it possible to make more quantifiable claims? For example, to say that "since the Gini coefficient for this population is twice as large, the variable is distributed eight times more unequally?" Or something like that? UPDATE The data isn't economic data - this has to do with relative numbers of contributions to a volutneer project. There are 5 datasets, with ~10,000 rows each. Each row corresponds to a single Internet user who did one or more tasks to contribute to that project. Some did just one, some did thousands. The question is: of the five projects, which ones are such that a small core group has done a lot of tasks each, and which ones are such that many people have done a few tasks each to add up to the total? FURTHER UPDATE Thinking about this more, maybe the Gini coefficient doesn't have enough information to make this kind of comparison. Maybe I need to compare the actual distributions of the variable among population members, using something like a K-S test?


title: Comparing sets of non-normal distributions

text: Apologies in advance if this question seems trivial -- I am a biologist and have little experience in statistics. I have six non-normal distributions (they're positively skewed): three are replicates from a drug treatment and three are replicates of an untreated control. I want to determine if the distributions from drug treated sample are significantly different from the control distributions. If I had only a single treated and control replicate, I would perform a Mann- Whitney U test to determine whether they were different. My confusion lies in how to appropriately handle the replicates (e.g. do I pool them?). Thanks for your time.


title: Forecasting unemployment rate

text: I have a data set of 100 geographic regions for which the unemployment rate has been observed during the last 9 years. Now, I want to simulate/forecast from this data the next year unemployment rate for all the 100 regions. What would be a suitable approach? I was thinking to use an auto-regressive mixed effect model, but I fear that I do not have enough data... Regards


title: Subsample size - when is it too small?

text: I ran into a problem (question) with biological data analysis and don't know how to solve it. Will try to explain problem as simple as possible: * * * **Dummy example:** There are 10000 textbooks in a library;   There are 1000 physics textbooks (10% of all the books);   I noticed that some of the physics textbooks have naughty drawings made by students;   Made a hypothesis that physics textbooks are enriched in naughty drawings, hence I did a permutation test (randomly resampled all the books in the library (10000) and got p value ~ 0). Indeed, physics textbooks are enriched in naughty drawings.   However, only 10 physics textbooks have those naughty drawings (0.1% of all the books in the library).   **Question:**   Is this subsample (0.1% of all the books) rational enough (statistically valid) to continue my analysis and draw conclusions about those naughty drawings physics textbooks (what field are those books from (e.g., molecular, optical physics), year of publication, etc.)? * * * **Actual problem:** There are 10000 genes in a genome;   My sample is 1000 genes (they are responsible for embryo development);   I noticed that some of those genes (10) have specific characteristic (sequence that is responsible for protein binding);   Did permutation test (randomly selected genes from all the genes in a genome) and got p value ~ 0. This lets my state that genes that are responsible for embryo development are enriched in specific characteristic compared to all the other genes in a genome.   Is it worth continuing to analyze those 10 genes? My sample is 1000 genes and this subsample is only 10 genes? Is this enough to draw any conclusions about biology when their are 10000 genes in a genome?


title: PCA on log transformed data

text: I am about to conduct two Principal Component Analyses (PCA) on species abundance data and species composition data. I have about 12 different locations where abundance data for over 50 invertebrate species has been collected. The software I am using has an option to log transform and scale the data prior to analyses. My question is, without spending too much time examining the data, is it a fairly common approach to both log transform and standardize the data prior to PCA analyses? Or is it best to just leave the data alone and run the PCA?


title: Data manipulation in R for functional data analysis

text: I am newbee and i am trying with functional data analysis. I have a 8x11 matrix data, how can i input into R as an object in this form:               $hgtm               boy01 boy02 boy03 boy04 boy05     1     81.3  76.2  76.8  74.1  74.2      1.25  84.2  80.4  79.8  78.4  76.3      1.5   86.4  83.2  82.6  82.6  78.3       1.75  88.9  85.4  84.7  85.4  80.3      2     91.4  87.6  86.7  88.1  82.2       3    101.1  97.0  94.2  98.6  89.4      4    109.5 104.6 100.4 104.4  96.9      5    115.8 112.3 107.1 111.0 104.1           $hgtf               girl01 girl02 girl03 girl04 girl05      1      76.2   74.6   78.2   77.7     76        1.25   80.4   78.0   81.8   80.5     80        1.5    83.3   82.0   85.4   83.3     83       1.75   85.7   86.9   87.9   87.0     86        2      87.7   90.0   89.6   90.3     89        3      96.0   94.9   97.1   98.6     96       4     103.8  102.1  109.2  106.3    103        5     110.7  109.2  116.2  113.9    110             $age           [1]  1.00  1.25  1.50  1.75  2.00  3.00  4.00  5.00      very much appreciated,


title: Difference between Excel's RAND(), RAND()*RAND(), etc

text: I plotted below the standarized results of:   * RAND()   * RAND() * RAND()    * ...   * RAND() * RAND() * RAND() * RAND() * RAND() * RAND() It seems that the results are getting to zero, is that because you're multiplying a bunch of numbers close to zero together, or is there another explanation? ![alt text](http://i.stack.imgur.com/5asaM.jpg) **EDIT** There is no real motivation behind this. Working on a spreadsheet earlier today, I was in a situation where i was multiplying a constant * Rand() * Rand() .. each rand() corresponded to a different variable (column), so I wanted to find out what goes on if you multiply a bunch of random numbers together.


title: How should I define missing values due to skip questions in SPSS?

text: I have a questionnaire that contains some skip questions. Like, say, the 3rd question is a yes/no type question. Only those who answered "yes" to the 3rd question are requested to answer the 4th, 5th and 6th question and those who answered "no" are requested to skip these three questions. I am sorry that I have no experience about how to make SPSS know that the missing values in the variables corresponding to the 4th, 5th and 6th question are due to the answer "no" in the previous (3rd) question. I should mention that some respondents didn't even answer the 3rd question and thus the answers to the 4th, 5th and 6th questions are automatically missing. So not all the missing values in these three variables (corresponding to the 4th, 5th and 6th question) are solely due to the answer "no" in the previous question. I need to know, if I need to use multiple imputation, then how can I avoid the imputation of the missing values due to the answer "no". I have seen something about user-defined missing values and I do not have much idea about it. Can anyone explain what should be done in my case? I will basically try a regression and later if I have time, a factor analysis with the data.


title: Why is the coefficient of variation not valid when using data with positive and negative values?

text: I can't seem to find a definitive answer to my question. My data consists of several plots with measured means varying from 0.27 to 0.57. In my case, all data values are positive, but the measurement itself is based on a ratio of reflectance values that can range from -1 to +1. The plots represent values of the NDVI, a remotely derived indicator of vegetation "productivity". My intention was to compare the variability of values at each plot, but since each plot has a different mean, I opted for using the CV to gauge the relative dispersion of NDVI values per plot. From what I understand, taking the CV of these plots is not kosher because each plot can have both positive and negative values. Why is it not appropriate to use the CV is such instances? What would be some alternatives?


title: How can I estimate coefficient standard errors when using ridge regression?

text: I am using ridge regression on highly multicollinear data. Using OLS I get large standard errors on the coefficients due to the multicollinearity. I know ridge regression is a way to deal with this problem, but in all the implementations of ridge regression that I've looked at, there are no standard errors reported for the coefficients. I would like some way of estimating how much the ridge regression is helping by seeing how much it is decreasing the standard errors of specific coefficients. Is there some way to estimate them in ridge regression?