title: Bootstrapped t test in SPSS: Why are some p value cells blank?

text: Using an independent samples t test, I am comparing the performance of 2 groups on an activity, based on the percentage of items correct (the percentage of total items correct, as well as percentage of items correct for four sub-domains of the activity). The n of one group is 23, and the n of the other is 34. Due to a ceiling effect, my data are non-normally distributed for both groups. At the recommendation of a statistician/professor, to correct for this problem I bootstrapped the t tests in SPSS using 95% bias corrected confidence intervals based on 1,000 replications. The smaller group has a larger variance, so I am interpreting the analyses using Welch's correction ("equal variances not assumed"). The SPSS output for the bootstrapped t tests produces a table with some p values left blank under the Sig. (2-tailed) column...Why is this? And what do I report instead of these missing bootstrapped p values? I just want to be consistent in how I explain which results were significant. All of the Lower and Upper BCa 95% confidence intervals ARE listed in the table, but I'm not sure how to succinctly communicate significance using only the confidence intervals (Though I know that if the interval does not include 0 then it can be concluded that the difference between the 2 groups is significant).


title: Building Hedonic model / regression

text: I'm not an economist and I am little confused. I want to build a hedonic model and regression for analyzing price evolution of some product (e.g. computers) on the market and hedonic price/quality indices. The analysis shoud be based on two points:   1. before entry of a new player on the market (period of two years) and    2. after entry of a new player on the market (period of three years). Some of the products in the "after" period have characteristics which products in the "before" period don't have. How should I do it: to gather computer data (prices, characteristics) for every year of the sample, then make regression equation based on different models (double-log, semi-log, linear), right? Then, this equation will serve to estimate computer prices that eventually were not part of the observed sample? How to find prices for products that have different characteristics (e.g. with CD ROM or without)? What exactly hedonic price/quality index represent? Please help me with some guidelines and advices. I've been studying hedonic model/regression for the last days and I'm all confused. Many thanks. J.


title: What could be a statistical test for comparing funnel data before and after important changes?

text: I have some funnel data from a website I manage. After some time I did some improvement on a multi-step purchase process. I wonder if it's possible to formulate a test that with statistical significance can assess if the improvements in the process have actually been effective. In other words I would like to test if the following two columns represent significantly different behaviours or not: \begin{array} {|r|r|r|} \hline \text{step} & \text{before changes} & \text{after changes} \\\ \hline \text{product page}&121 &1349 \\\ \text{clicked on purchase}&32 &412 \\\ \text{step 3}&29 &391 \\\ \text{step 4}&28 &386 \\\ \text{step 5}&26 &371 \\\ \text{transaction completed}&19 & 345 \\\ \hline \end{array} In that case you can see that the number of initial users of the first column is very different from the second column. But this is not what I would like to test: I would like to "demonstrate" that the subsequent behaviour contains relevant differences. I would also repeat the test with data that might come from different segments of the users. Any ideas?


title: Using glmer to estimate treatment interactions

text: In my data, I have two treatment conditions with repeated measures for each subject. I would like to run a mixed logistic regression separately for each of my two conditions where my binary outcome DV (dependent variable) is regressed on my IV (independent variable) and also have a random slope and intercept fitted for each subject. So, I run the following:               modelT0 <- glmer(DV ~ IV + (1|subject) + (0 + IV|subject), data = D0, family = binomial)     modelT1 <- glmer(DV ~ IV + (1|subject) + (0 + IV|subject), data = D1, family = binomial)      In the above, D0 and D1 are data sets restricted to treatment conditions 0 and 1, respectively. What I would like to do is compare the estimated fixed effects coefficient on IV across conditions to see if it significantly changes. To do this, I pool D0 and D1 into a single data set, D, and create a treatment indicator that takes value 0 in D0 and 1 in D1. I then run:               model <- glmer(DV ~ IV + treatment + treatment:IV + (1 + treatment|subject:treatment)                     + (0 + IV + treatment:IV|subject:treatment), data = D, family = binomial)      I should be able to look at the fixed effects coefficient on treatment:IV to get my answer, but the issue is that for whatever combination of random effects I seem to specify, the coefficients from the pooled regression are slightly different from the regressions specified by treatment. So for instance, the fixed effect coefficient on treatment:IV plus the one on IV in model is not equal to the coefficient on IV in model1. Any idea about what I might be doing wrong or how to answer the question I have? Thanks! EDIT: As per Henrik's suggestion, I'm copying the random effects output of the models below: summary(modelT0):                   Random effects:         Groups    Name        Variance  Std.Dev.          subject   (Intercept) 1.412e-07 0.0003758         subject.1 IV          1.650e+00 1.2844341      summary(modelT1):                   Random effects:         Groups    Name        Variance Std.Dev.         subject   (Intercept) 0.00378  0.06148          subject.1 IV          0.26398  0.51379       summary(model):                   Random effects:         Groups              Name         Variance  Std.Dev. Corr          subject.treatment   (Intercept)  0.0005554 0.02357                                    treatment    0.0066042 0.08127  -0.88         subject.treatment.1 IV           1.6500112 1.28453                                    IV:treatment 1.0278663 1.01384  -0.93


title: Triangular distribution and simulation monte carlo

text: How can I do to incorporate Monte Carlo simulation on a triangular function. I'm trying to do this using JavaScript and found a very interesting bibliioteca http://jstat.github.io/distributions.html # jStat.triangular.pdf however I do not know how should put, I think you should use value a, b, c and random that would be generated by Mersenne Twitter, but I'm a little lost


Can we use clustering output as predictor variable for classification?   1. I have a set of data and I do clustering analysis on it, it divides the data into different clusters.   2. Can I use this cluster information i.e. cluster1, cluster2, cluster3, as one of the input variable i.e. predictor variable for my decision tree algorithm is it statistically OK? I found great increase in the model accuracy when I do this


title: How to convert nominal dataset into numerical dataset?

text: For my work, im using the multilabel dataset from this webpage. Few dataset which are listed in the page (for, e.g _bibtex_ ) have nominal attributes, i.e attribute values are 0 and 1. My queries are given below   1. Is it valid to run kMeans clustering algorithm on these nominal dataset to get centers and target label which are meaningful?   2. Otherwise, to run kMeans algorithm (forget abt the target label), i need to convert this nominal dataset into numerical dataset. What is the standard procedure of doing it. I can normalize each instance, but it just gives me a real number with equal value for an instance.   3. I would also like to reduce the dimension of nominal dataset such as rcv1v2. How do i go about it. I can use any Feature selection technique but it requires an optimization criteria. But in my case, i need to compare the result of different algorithm on this dataset which have different optimization criteria, so i got into trouble of choosing which criteria. Is there any technique of selecting a top features?


title: How to compare joint distribution to product of marginal distributions?

text: I have two finite-sampled signals, $x_1$ and $x_2$, and I want to check for statistical independence. I know that for two statistically independent signals, their joint probability distribution is a product of the two marginal distributions. I have been advised to use histograms in order to approximate the distributions. Here's a small example.               x1 = rand(1, 50);     x2 = randn(1, 50);     n1 = hist(x1);     n2 = hist(x2);     n3 = hist3([x1' x2']);      Since I am using the default number of bins, `n1` and `n2` are 10-element vectors, and `n3` is a 10x10 matrix. My question is this: How do I check whether `n3` is in fact a product of `n1` and `n2`? Do I use an outer product? And if I do, should I use `x1'*x2` or `x1*x2'`? And why? Also, I have noticed that `hist` returns the number of elements (frequency) of elements in each bin? Should this be normalized in any way? (I haven't exactly understood how `hist3` works either..) Thank you very much for your help. I'm really new to statistics so some explanatory answers would really help.


title: Simulation of Random Variables (Proof)

text: Given U is an uniform random variable on [0,1] and F is a cdf with inverse $F^{-1}$. I want to prove that $F^{-1}(U)$ has distribution F. Is the following proof correct: $F_{X}^{-1}(u) = x$ <=> $F_{X}(x) = u$ Proof: $F_{X}(x) = F_{U}(u)$ => $\int_{0}^{u} P_{u}(u)du = \int_{0}^{1} 1 du = u.$ Is this correct?


title: Layperson's guide to valid polling via Internet?

text: Sites such as Survey Monkey allow a layperson to easily create and administer a questionnaire. However, the only thing I know about polling is that you can't just create a questionnaire, let anyone answer it, and get valid results from it. Is there a layperson-accessible guide to the problems and (hopefully) solutions to creating Internet-based polls from which valid conclusions can be drawn? If not, what would be the best route for a small-business owner unable to afford hiring a national polling company? (I'm not sure that this is on-topic for Cross Validated, but it seemed like the closest match...)


title: Inference from linear regression slope and Pearson

text: Sorry if this has been asked before but I've already done quite a bit of work here and I feel like I'm quite close to an answer. I am interested in testing whether the PHP function array_key_exists depends on the size of the array. If it did, we should expect that the growth of the computation time for the function would depend on the size of the array. If it does not, then the computation time should be constant. In other words, I want to test the strength of a possible linear relationship. I measured the computation time for varying sizes of arrays. A few things I see, in ascending order of importance:   1. There lots of points (almost half of them) where the time actually decreases as the size increases. This basically gives me grounds to accept the null hypothesis immediately, because _any_ linear relationship would entail a monotonically increasing function of the size.    2. Very low standard deviation in the time, indicating that the values are tightly clustered.   3. The slope of the simple linear regression is quite small.   4. Exceedingly low Pearson coefficient and its square has practically vanished.  If my null hypothesis is that the _real_ Pearson coefficient is zero, namely there is no relationship between the array size and computation time, how small does my _sample_ Pearson coefficient or its square need to be to reject it? What assumptions do I need to make (perhaps acceptable probabilities of Type I and Type II) to do so? Is it better to formulate the null hypothesis from the Pearson coefficient or from the linear regression slope? Are they about the same as far as explanatory power in this circumstance?


title: "Since $x$ is near-gaussian, its PDF can be written as..."

text: **Short question:** Why is this true?? **Long question:** Very simply, I am trying to figure out what justifies this first equation. The author of the book I am reading, (context here if you want it, but not necessary), claims the following: > Due to the assumption of near-gaussianity, we can write: > > $$ p_0(\xi) = A \; \phi(\xi) \; exp( a_{n+1}\xi + (a_{n+2} + > \frac{1}{2})\xi^2 + \sum_{i=1}^{n} a_i G_i(\xi)) $$ Where $p_0(\xi)$ is the PDF of your observed data that has maximum entropy, given that you had only observed a series of expectations, (simple numbers) $c_i, i = 1 ... n$, where $c_i = \mathbb{E}\\{G_i(\xi)\\}$, and $\phi(\xi)$ is the PDF of a standardized gaussian variable, that is, 0 mean, and unit variance. Where all this is going is that he uses the above equation as a starting point for making the PDF, $p_0(\xi)$ simpler, and I get how he does it, but I do not get how he justifies the above equation, ie, the starting point. I have tried to keep it brief to as not to obfuscate anyone, but if you want additional details please let me know in the comments. Thanks!


title: Help me fit this non-linear multiple regression that has defied all previous efforts

text: **EDIT: Since making this post, I have followed up with an additional posthere.** **Summary of the text below: I am working on a model and have tried linear regression, Box Cox transformations and GAM but have not made much progress** Using `R`, I am currently working on a model to predict the success of minor league baseball players at the major league (MLB) level. The dependent variable, offensive career wins above replacement (oWAR), is a proxy for success at the MLB level and is measured as the sum of offensive contributions for every play the player is involved in over the course of his career (details here - http://www.fangraphs.com/library/misc/war/). The independent variables are z-scored minor league offensive variables for statistics that are thought to be important predictors of success at the major league level including age (players with more success at a younger age tend to be better prospects), strike out rate [SOPct], walk rate [BBrate] and adjusted production (a global measure of offensive production). Additionally, since there are multiple levels of the minor leagues, I have included dummy variables for the minor league level of play (Double A, High A, Low A, Rookie and Short Season with Triple A [the highest level before the major leagues] as the reference variable]). Note: I have re-scaled WAR to to be a variable that goes from 0 to 1. The variable scatterplot is as follows: ![scatterplot](http://i.imgur.com/IJNowem.jpg) For reference, the dependent variable, oWAR, has the following plot: ![dependentvariableplot](http://i.stack.imgur.com/AMXDm.jpg) I started with a linear regression `oWAR = B1zAge + B2zSOPct + B3zBBPct + B4zAdjProd + B5DoubleA + B6HighA + B7LowA + B8Rookie + B9ShortSeason` and obtained the following diagnostics plots: ![linearRegressionDiagnostics](http://i.imgur.com/U6ABt5a.jpg) There are clear problems with a lack of unbiasedness of the residuals and a lack of random variation. Additionally, the residuals are not normal. The results of the regression are shown below: ![linearRegressionResults](http://i.imgur.com/3kW3ZAO.jpg) Following the advice in a previous thread, I tried a Box-Cox transformation with no success. Next, I tried a GAM with a log link and received these plots: ![splines](http://i.imgur.com/hWQAS2d.jpg) **Original** ![diagnosticChecksGAM](http://i.imgur.com/y6KVJtB.jpg) **New Diagnostic Plot** ![GAMDiag](http://i.imgur.com/UKD2ycB.jpg) It looks like the splines helped fit the data but the diagnostic plots still show a poor fit. **EDIT: I thought I was looking at the residuals vs fitted values originally but I was incorrect. The plot that was originally shown is marked as Original (above) and the plot I uploaded afterwards is marked as New Diagnostic Plot (also above)** ![GAMResults](http://i.imgur.com/78BObfU.jpg) The $R^2$ of the model has increased but the results produced by the command `gam.check(myregression, k.rep = 1000)` are not that promising. ![GAMResults2](http://i.imgur.com/XoubLbh.jpg) Can anyone suggest a next step for this model? I am happy to provide any other information that you think might be useful to understand the progress I've made thus far. Thanks for any help you can provide.


title: Number of folds for K-fold

text: Computing power considerations aside, are there any reasons to believe that **increasing the number of folds** leads to better model selection/validation? (the higher the number of folds the better) Taking the argument to the extreme, does **leave-one out** necessarily lead to better models than **K-fold**? Some background on this question: I am working on a problem with very few instances (e.g. 10 positives and 10 negatives), and am afraid that my models may not generalize well/would overfit with so little data.


title: Unsupervised classification of curves (in R)

text: I'd like to know if there is a method (in R) to classify curves based mainly on their appearance ie one cluster will represent a specific form of curves. Thanks


title: Using text mining/natural language processing tools for econometrics

text: I am not sure whether this question is fully appropriate here, if not, please delete. I am a grad student in economics. For a project which investigates issues in social insurances, I have access to a large number of administrative case reports (>200k) which deal with eligibility evaluations. These reports can possibly be linked to individual administrative information. I want to extract information from these reports that can be used in quantitative analysis, and ideally more than simple keyword/regex searches using `grep`/`awk` etc. How useful is Natural Language Processing for this? What are other useful text-mining approaches? From what I understand this is a large field, and most likely some of the reports would have to be transformed to be used as a corpus. Is it worth investing some time to become acquainted with the literature and methods? Can it be helpful and has something similar been done before? Is it worth it in terms of the rewards, i.e. can I extract potentially useful information using NLP for an empirical study in economics? There is possibly funding to hire somebody to read and prep some of the reports. This is a larger project and there is a possibility to apply for more funding. I can provide more details about the topic if strictly necessary. One potential complication is that the language is German, not English. Regarding qualifications, I am mostly trained in econometrics, and have some knowledge about computational statistics at the level of the Hastie et al. book. I know Python, R, Stata, and could probably get familiar with Matlab quickly. Given the libraries, I assume Python is the tool of choice for this. No training at all in qualitative methods if this is relevant, but I know some people I could reach out to. I am glad for any input on this, i.e. if this is potentially useful, if so, where to start reading and which tools to focus on in particular.


title: Creating a Normalization Factor

text: I have a relatively simple problem that I can't seem to find a satisfactory solution to. If I have three scales for three different sets of data. One varies from [-5,5] the other from [1,10], and the last from [-10,10] what would an appropriate way to make a single scale for data from all three sets of data?


title: what do you call a pie chart with equal arcs but different areas

text: I want to create a pie chart for four items, the arc of the slices would be equal, each 90 degrees, however the areas would be different since I will shrink or enlarge the slice depending on the data. Here is an that is close to what I want: ![\[example\]\(http://www.excelcharts.com/blog/creative-pie-charts-better- infographics/#!wp- prettyPhoto\[3934\]/0/\)](http://i.stack.imgur.com/HDOPk.png) Basically it would look like a circular bar graph that is filled in the center. What do you call this type of chart and how would I create it.


title: How do I interpret the results from the F-test in excel

text: I have a fairly simple question regarding the interpretation of the F-test in Microsoft Excel. Let't say these are the results of my F-test: ![enter image description here](http://i.stack.imgur.com/7Ds80.jpg) I am now wondering how to interpret it in order to choose the correct t-test (assuming equal or unequal variances) for my data-set. I have found guides telling me if F critical > F, then use unequal variances. However, some of the guides tell you to use only the p value, so I am unsure which parameters to look at when interpreting the results.


title: Difference between exploratory and confirmatory factor analysis in determining construct independence

text: Researchers often use two measures that have very similar items and argue that they measure different things (e.g., "I always worry when I am around cars"; "I am fearful of cars"). Lets call the hypothetical measures the Fear of Cars Measure and Anxiety from Automobiles Scale. I am interested in testing empirically if they indeed assess different latent constructs, or if they measure the same thing. The two best ways I can think to do this would be through exploratory factory analyses (EFA) or confirmatory factor analysis (CFA). I think EFA would be good because it allows all of the items to load freely without constraints. If items from the two scales load on the same factors, then I can conclude that the measures likely don't assess different things very well. I can also see the benefits in CFA, however, since I will be testing pre-defined models. For example, I could compare the fit of a model in which all items load onto a single factor (i.e., they don't assess different constructs) or the items are separated into the expected measures. An issue with CFA, I suppose, is that it would not really consider alternative models (e.g., a three factor model). For the purposes of discussion, lets also perhaps consider that there may be two other very similar measures out there (e.g., Car anxiety questionnaire and the Scales for the assessment of car fears) that I wish to throw into the mix! How can I best statistically determine if two measures assess different constructs?


title: Forward Stepwise selection

text: I am assuming the following model: $Y = \beta X + \epsilon$ Here both $X$ and $Y$ are matrices. I fit the least squares model without any regularization and get the matrix $\beta$. I would like to do the following now. I want to pick some subset $S$ of rows of $\beta$ using forward stepwise selection such that these rows minimize the squared loss. How should one go about it?


title: Supervised classification on different time series

text: I have 300 files, each file has a time series data with a class label(0 or1) for each data point.I want to build a classifier, which can predict the class of a new time series data.   1. How should I build my training and testing data set and what model I might use. Should I combine all the files and then do a 80/20 split or is there a better way?   2. Isn't this similar to classification for cross-sectional data, in other words How do I account for the time properties of the series.


title: Constant variance assumption in regression model

text: My question is how do we check the constant variance assumption in a **regression model**?


title: Can lognormal + constant be still considered lognormal?

text: In one of my computations (result of GLM with some model values replaced with observations) I resulted with a distribution that is a lognormal distribution **plus** some constant. Can I still speak of lognormal distribution? Again, this question is more about the statistical language and convention than about the statistics itself...


title: Visualizing predicted probabilities from logistic regression with a weighted sample in Stata

text: I was using the `predxcon` command to graph the predicted probabilities for logistic regression. Basically I want to see the pattern for predicted probabilities of the outcome by each age group. Using `predxcon` command, it provided the graph, as well as corresponding values. However, I realized that the sample is weighted and I need to use `svy` to adjust for the weights. I am wondering if there is any program like `predxcon` when using `svy`? I am planning to calculate the predicted probabilities based on coefficients provided by the regression result with adjustment by the weights. Is there anything I should be aware of, since this is `svy`? Do I need to consider weights or I can just use the coefficients for each predictor?


title: Find degradation rate on a scatterplot

text: I work for a commercial printing company, and am collecting data regarding the quality of our printing over time (density and dot gains). Essentially, for every job I have a date, a density, and a dot gain. I would like to be able to analyze a date range and check to see if our quality seems to be getting consistently worse. From a math perspective, I somewhat understand what I need to do. I need to get a best fit curve for our data, and then find the derivative of that curve. The slope of that derivative should give some insight as to if our quality is getting better or worse, and how quickly. My biggest hiccup here is how to find the best fit curve for a scatter plot. Certainly I can work out a way to find the derivative after that, but I can't think of any good way to find the original curve. So - how do you find the best fit curve of a scatter plot? And bonus points if you can give me a quick and (relatively) easy way to get the derivative/slope of that curve.


title: Markov chain long run probabilities

text: So I want to find the long run probability of getting two heads in a row when flipping a fair coin many times. I know this answer should be 1/6 (expected number of flips to get 2H in a row is 6) but I'm trying to figure out how to use a Markov chain to solve this. I think my transition matrix should be (3x3):                   0H 1H 2H     0H  .5 .5  0     1H  .5  0  .5     2H  0   0  1      Then I come up with a system of equations to find the limiting distribution.               .5x1 + .5x2 = x1     .5x1 = x2     .5x2 + .5x3 = x3     x1 + x2 + x3 = 1      But I can't seem to solve this equation correctly. Is my transition matrix/ approach correct?


title: Expected number of cards drawn

text: You draw cards from a standard 52-card deck without replacement until you get a **queen of spades** and stop. The cards have the values $2,3,4,\ldots,11(J),12(Q),13(K),14(A)$. The first card you draw you give it to your friend and from then on, whenever you draw a card with value smaller than the card you last gave to your friend, you give the new card to your friend . What is the expected minimum value among the cards you give to your friend and what is the expected number of cards you will give to your friend? The minimum value is the last card given to your friend. To find the expected number of cards will it be $4\frac{4}{5}+4\frac{4}{5}\frac{12}{13}+...$, since I can only give a 2 if it is the first card drawn, I can only give a 3 if 2 is not the first card drawn and 3 is before the Queen of Spades...? In this last case can I just multiply the 2 probabilities as if the events were independent?


title: Two-Way ANOVA with non-sig interaction effect: is it reasonable to argue for an interaction from simple effects?

text: This is in the context of marking an undergraduate sociology essay in which the student has not found the predicted significant interaction effect in a Two-Way ANOVA. The diagram below is schematic and doesn't depict the actual data. Imagine that the interaction effect is non-significant. The student has said that because the effect of task is significant in the treatment group but non- significant in the control group, there's some evidence that there really is an interaction effect. ![enter image description here](http://i.stack.imgur.com/c9Kxl.gif) I sense there is something fishy about this but I am struggling to come up with a precise explanation of why.


title: Proportion in regression?

text: How do I interpret a proportion in a regression model where the proportion is an independent variable (B1)? For instance, let's say B1 is 77 (or .77): how would one interpret that?


title: relation between median survival and HR

text: What would be the mathematical relation between 2 comparing survival medians (difference or relation of medians) and the HR? I know that the relation of medians of survival time (MSTc/MSTi), were "i" is intervention and "c", control, is not equal to HR(i/c), although there is a correlation between these two numbers. On other side, I found that Si(t)=Sc(t)^HR, where Si(t) is the proportion of survival in intervention group at a given time, "t". But could we calculate a theoretical MSTi, given an estimated MSTc and HR(i/c)? Or could we calculate MSTc/MSTi -the relation bewteen medians-, given the HR (i/c)? It would be interesting, as observed difference of medians is erratic. Moreover, Case et al (The Oncologist 2002;7:181-7), wrote: "Note that under the proportional hazards assumption, the ratio of the logs of the survival times would equal the HR". Does it mean that log(MSTi)/log(MSTc)=HR? I don't think so.


title: An inequality involving expectation

text: Let $f,g$ be two pdfs, and suppose $X$ is a random variable that has pdf $f$. Is it necessarily true that $E[f(X)] \ge E[g(X)]$? Although I doubt this will help, but I got this problem from studying the Kullback-Leibler divergence, which is defined as $D(f,g) = E[ln (f(X)/g(X))]$ (with $X$ having pdf $f$ like above). It can be shown using Jensen's inequality that $D(f,g) \ge 0$, which is equivalent to saying that $E[ln(f(X))] \ge E[ln(g(X))]$. But my question removes the $ln$ and I am wondering whether it will still be true. In fact, I hypothesize that not only is it true for the $ln$ function, but it will be true for any function that is increasing (which includes the identity function, which is my question here). So as a bonus, maybe we could ask if $E[h(f(X))] \ge E[h(g(X))]$ where $h$ is any strictly increasing function on positive reals. It seems like a continuous form of the rearrangement inequality to me.


title: Goodness-of-fit for very large sample sizes

text: I collect very large samples (>1,000,000) of categorical data each day and want to see the data looks "significantly" different between days to detect errors in data collection. I thought using a good-of-fit test (in particular, a G-test) would be a good fit (pun intended) for this. The expected distribution is given by the distribution of the previous day. But, because my sample sizes are so large, the test has very high power and gives off many false positives. That is to say, even a very minor daily fluctuation will give a near-zero p-value. I ended up multiplying my test statistic by some constant (0.001), which has the nice interpretation of sampling the data at that rate. This article seems to agree with this approach. They say that: > Chi square is most reliable with samples of between roughly 100 to 2500 > people I'm looking for some more authoritative comments about this. Or perhaps some alternative solutions to false positives when running statistical tests on large data sets.


title: What if P(A|B)=P(B|A)?

text: If we have P(A|B)=P(B|A), then what is this special case called and are there special properties? I'm interested in a simpler way of computing one of them and would like to take advantage of such properties.


title: Correct for multiple testing using Spearman

text: I am using a Spearman correlation for my dataset. I have three conditions and I want to know whether there is a correlation per condition. To do this I filter e.g. condition 1 and correlate variable X and variable Y. After I do the same for condition 2 and condition 3. Do I have to correct for the p-values? If yes, should I correct with Bonferroni? Thanks a lot!


title: How to use gsl_linalg_cholesky_ decomp to composite the lower triangular matrix?

text: As we know, Cholesky decomposition of $A = L*L^T$. I tried to write a simple function to decompose the lower triangular matrix $L$. I know there is a C++ function of `GSL/gsl_linalg_cholesky_ decomp` that can do it. I read its manual but I do not quite understand. Anyone can help? * * * I was asked to offer the manual of the function, which is http://www.gnu.org/software/gsl/manual/html_node/Cholesky-Decomposition.html# 1 index-gsl 005flinalg 005fcholesky 005fdecomp-1343


title: Best way to report "average speedup"?

text: I have a collection of benchmarks that we run against different implementations of a system. Is there a generally-accepted way to make comparisons between multiple sets of benchmark data? Today, we take the data for system $a$ (with data points $a_1, a_2, a_3,...$) and the data for system $b$ (with data points $b_1, b_2, b_3,...$) and take the geomean of the speedup (e.g. $\frac{a_1}{b_1}$) for each point. Is this a valid approach? In this case, it is safe to assume that the data points have equal importance.


title: Anisotropy in kriging for non gridded data

text: I have to perform a mapping of a DVB-T field (the TV signal), per every location I consider the median in time of the measurements, there are some issues, e.g. the variance seems to be proportional to the intensity of the field because of the measuring instrumentation. My first idea was to use kriging. Following R's gstat package vignette I have plotted the directional variograms and one of the variables that are measured shows a clear anisotropy.   1. How do I deal with it? Is there a favourite framework?    2. Looking at literature a change in variable is suggested but since my data are not in a grid I do not know how to deal with it. I have looked at this post Problems estimating anisotropy parameters for a spatial model   3. How do you know from the cloud variogram if there are enough points to consider aisotropy? In a second dataset which is less numerous a low number of samples for a range of distance, that is lower than the range corresponds to a drop in the variogram. I am not sure if this is a phisiological drop due to the field or if it is due to the low-sampling.   4. Is there any way to find out?   5. For doing this kind of analysis I am using gstat in R. Is geoR a better framework? Or is there a better package in absolute? There are two measurement campaigns, one with 110 locations, the other with 35, in a territory which is approximately 4 km * 4 km, but I have also simulated the field, so from the simulated field I can extract as many samples as I want. Any pointer to relavant literature is very welcome. At the moment I am reading [Webster, Oliver] _Geostatistics for Environmental Scientists_.


title: Avoiding p-values and reporting t-values instead. References?

text: I was wondering if I could get some opinions on an issue. I'm analyzing my data using mixed-effects modeling in R (lme4 package). My model has by-subject and by-item intercepts and slopes, and random correlation parameters between them. Since the current version of lmer() does not have MCMC sampling implemented, I cannot get a pvalue for the coefficients in the model. Therefore, I would like to report the **t-value** instead. I have often seen papers in my field (psycholinguistics) just say something like "In all models presented, |t| > 2 and |z| > 2 correspond to a significant effect at a significance level of .05". I was wondering whether there is some reference I can provide from this type of sentences? I understand that this tends to be the case, but I wonder whether this is something that has been shown (and I should give references) or whether it is ok to just state it and assume everyone will be ok with it. Suggestions welcome!


title: What's the value of information when we decrease the entropy of a probability distribution?

text: Suppose you have to choose between actions $A_1,\dots,A_n$. You have a probability distribution over each $U(A_i)$, i.e. over the utility of choosing each action. So you should choose the $A_i$ that maximizes $\mathbf E[U(A_i)]$. But now you are given the opportunity to reduce your uncertainty about one of these $U(A_i)$. I think the value of this new evidence should be $$\int_{e\in E}P(e)~\mathbf E[U'(A_i)]~\text{d}e$$ where $E$ contains every possible piece of evidence you could receive, $P(e)$ is how likely it is, and $U'(A_i)$ is the utility of action $A_i$ after you receive the evidence $e$. And I know I'm pushing it, but: can you think of some way to _roughly estimate_ the value of reducing entropy, given only the expectation and current entropy of each $U(A_i)$, and the amount by which you will reduce the entropy of one particular $U(A_i)$? For instance, if $U(A_1)$ has expectation 20 and entropy 10 bits, and you're given the same information about every other $A_i$, how much would you pay to reduce the entropy of $U(A_i)$ by 2 bits?


title: Could I compare two crosstabs and chi square tests?

text: For example, I have two date sets from city A and B. I constructed two crosstabs and perform two chi square tests for both cities respectively, about if they have significant associations between gender and their participation in recreation activities (eg. cycling, camping). The results shows there are significant associations between both variables. Is it logical if I: compare males between city A and city B in their participation of a certain activities, ie. cycling? compare if which city's male are more likely to participate more in a certain activities? Thank you.


title: Third order central moment of a positive linear combination of log-normal random variables

text: I asked this question on MathOverflow a couple of days ago. What is the sign ($\pm$) of the third order central moment of a positive linear combination of log-normal random variables? It seems to be a common notion that the skewness of random variables with longer tails to the right is positive. Is it correct? If so, how do you prove it?


title: Modeling time in multilevel logistic regression

text: I conducted an experiment in which participants listened to sentences while looking at pictures about the sentences on a computer screen. Whether at a given time point a participant looked at the left half of a picture or the right half of a picture was recorded via an eye-tracker. The study had a 2 by 2 within subjects and within items design. If time was not an issue here, I could use the model below (using `R` code). (Here, I assume varying random intercepts and random slopes for both subject and item.)               lmer(look ~ iv1 * iv2 + (1 + iv1*iv2 | subject) + (1 + iv1*iv2 | item),           family = "binomial", data=data)      However, given that participants listened to sentences each of which lasted a few seconds, where they looked at (left vs. right) would vary across time. So one way to model time might be to include `time` as a covariate (it may even be necessary to natural polynomials):               lmer(look ~ iv1 * iv2 + time + (1 + iv1*iv2 | subject) + (1 + iv1*iv2 | item),           family = "binomial", data=data)      But things get complicated since (1) it is possible that `time` interacts `v1` and `v2` and (2) I probably need to somehow model `time` in the random effect terms. An additional complication is that eye-tracking data is fairly large. The particular set of data that I am working with currently has 5 million rows, so running even the simplest multilevel logistic regression can be fairly time consuming. So my question is, given my design, what would be a good way to model time.


title: Proper proposal distribution in Metropolis-Hastings for a 2-state discrete process

text: This question is related to Understanding Metropolis-Hastings with asymmetric proposal distribution. I'm trying to draw sample paths from a 2-state markov chain using M-H. Unfortunately, I've only studied the Gaussian (random walk) proposal distribution. Is there a generally accepted (best) proposal distribution to use when drawing candidates from a discrete distribution? Given the prior estimates of the transition matrix, should I just compute the stationary distribution and draw unconditionally from them, thereby making it an independent chain? Or is it more typical to draw from a pre-specified symmetric transition matrix for the 2 states? Or is there another method which tends to work better? Note: I know that there are other ways to estimate the evolution of the 2 states. However, I use M-H in other steps of the model I'm working with> Therefore, I'd like to use it to estimate the time series of the 2-state process as well, if at all possible.


title: Non-parametric equality of means test

text: I'm looking to test for equality of means across different sample sizes of data, but know that the data is not normally distributed and heteroscedastic. Can anyone suggest anything?


title: How to get from $R$ to $R^2$ the hard way

text: If $r = \frac{S_{X, Y}}{S_X S_Y}$ and $r^2 = \frac{SS_{reg}}{SS_{tot}}$, how do I prove that squaring $r$ really gets me to $r^2$ using those definitions or derivations thereof? I guess the best way is to start by equating $\left[\frac{\sum{xy} - n\bar{x}\bar{y}}{\sqrt{\sum{x^2} - n\bar{x}^2}\sqrt{\sum{y^2} - n\bar{y}^2}}\right]^2$ and $\frac{\sum{(\hat{y}-\bar{y})^2}}{\sum{(y-\bar{y})^2}}$, but I can't get much past that. How can I further simplify this equation in order to prove that $\left(\frac{S_{X, Y}}{S_X S_Y}\right)^2 = \frac{SS_{reg}}{SS_{tot}}$?


title: can the same item be used in multiple constructs?

text: I am constructing a survey with around 10 multi-item scales. My piloting indicates that one of the most important multi-item scales (COP) has very poor reliability (.49), but there is an item in another scale (REX) which seems to hang together semantically with this COP scale. Including this item in the problem scale brings the reliability coefficient up to .74, but taking it away from the REX scale brings that alpha down to .53... Because I don't really want to go back to the drawing board and rewrite the items in the COP scale, and then re-pilot it, this seems to offer me a shortcut. I'm wondering if computing a variable (COP) with this item that will also be included in another variable (REX) is possible, or if it is just not done that way. Can single items be used multiple times - that is, in multiple scales of a single survey? Any opinions or help much appreciated.


title: How to test correlations between ordinal outcomes and various kinds of variables?

text: Running some analysis for 65 patients who have gone through a medical treatment post-diving. I have a dependent variable which is ordinal: 1-bad outcome, 2-okay outcome, 3-good outcome. I want to test some correlations with different independent variables:   1. Sex - nominal obviously. Should I use the Kruskal-Wallis test?   2. Dive depth - which is numeric in meters. Can I use univariate analysis?   3. Pre-treatment severity - ordinal variable as 0-no symptoms, 1-mild, 2-moderate, 3-severe. Kruskal-Wallis as well?


title: Fuzzy logic based object recognition

text: I need to classify objects using fuzzy logic. Each object is characterized by 4 features - {size, shape, color, texture}. Each feature is fuzzified by linguistic terms and some membership function. The problem is I am unable to understand how to defuzzify such that I may know which class an unknown object belongs to. Using the Mamdani Max-Min inference, can somebody help in solving this issue? Objects = {Object1, Object2, Object3, Object4} or denoted as {1,2,3,4} respectively. The fuzzy sets for each feature is : Feature : Size $\tilde{Size_{Large}}$ = {1//1,1/2,0/3,0.6/4} for crisp values in range 10cm - 20 cm $\tilde{Size_{Small}}$ = {0/1,0/2,1/3,0.4/4} (4cm - 10cm) Shape: $\tilde{Shape_{Square}}$ = {0.9/1, 0/2,0/3,0/4} for crisp values in range 50-100 $\tilde{Shape_{Cylindrical}}$ = {0.1/1, 1/2,1/3,1/4} (10-40) Feature : Color $\tilde{Color_{Reddish}}$ = {0/1, 0.8/2, 0.6/3,0.3/4} say red values in between 10-50 (not sure, assuming) $\tilde{Color_{Greenish}}$ = {1/1, 0.2/2, 0.4/3, 0.7/4} say color values in 100-200 Feature : Texture $\tilde{Tex_{Coarse}}$ = {0.2/1, 0.2/2,0/3,0.5/4} if texture crisp values 10-20 $\tilde{Tex_{Shiny}}$ = {0.8/1, 0.8/2, 1/3, 0.5/4} 30-40 The If then else rules for classification are R1: IF object is large in size AND cylindrical shape AND greenish in color AND coarse in texture THEN object is Dustbin (say Object1 = Dustbin from prior knowledge) or in tabular form just to save space               Object type  Size   Shape          Color        Texture     Dustbin :   Large  cylindrical       greenish   coarse     Can :       small  cylindrical       reddish    shiny     Bottle:     small  cylindrical        reddish    shiny     Cup :       small  cylindrical       greenish   shiny      Then, there is an unknown feature with crisp values X = {12cm, 52,120,11}. How do I classify it? Or is my understanding incorrect, that I need to reformulate the entire thing?


title: Homoscedastic and heteroscedastic data and regression models

text: How to understand the homoscedasticity and heteroscedasticity in context of regression models?   Is there a way to check these properties in R?


title: How can I find a Z score from a p-value?

text: How can I find a _Z_ score from a p-value? I know how to look up the p-value from a _Z_ score using a Normal distribution table, but I don't know how to calculate it. For example, a question says the alpha equals 5 percent. From this, I see in my handout that the _Z_ score is calculated to be 1.65. How do I determine this? Thank you.


title: Proof without Jensen

text: In the middle of an argument that I will present to my students I'll have to prove that $E[X^2]\geq E^2[X]$, but I don't want to use Jensen's inequality to do so. Is there any elementary way to go? Thanks!


title: Multiple testing problem - can these test results be caused by chance?

text: I conduct 30 independent significance tests on $p<.05$. I observe that in $X=3$ cases the test is rejected. Can this event be caused by chance due to a multiple testing problem? I understand that $X$ is binomially distributed, with expectation $E(X)=np=1.5$. What is the null hypothesis I should test, is it $H_0: p>.05?$. And how should this test be executed? EDIT: I think it is $H_0: p<.05$, so we evaluate $P(X>3)=1-P(X <= 3)=.061$. So this is close, but one would not yet reject. From @whuber's comment there seems to be a more sophisticated approach that I need to consider.


title: Factor analysis Vs PCA

text: Which is a better method of data reduction - factor analysis or principal component analysis?


title: Correlation coefficient for sets with non-linear corelation

text: What method can I use to test if there is a correlation between two sets of data? The correlation coefficient works if there is a linear association, but if I have two sets that are clearly (visually by graph) correlated in a non- linear way, how can I test that? Is there a coefficient or a special method?


title: SVM parameter dependence on number of samples

text: I need to do a grid search to optimize SVM parameters gamma, C and epsilon (svm from e1071 r package). The problem is that I have a fairly large data set, about 100000 rows and 40 variables. I have concluded that I can probably survive grid search and cross validation on 40000 samples of data. But can optimized parameters for subset of 40000 data be used on the final model of 100000 or are the parameters dependent on sample size, and how? Is any of the 3 parameters independent of other two, or at least in some stable range? For example, lets say I optimize just for gamma and C by holding epsilon constant. I know that C and epsilon influence the complexity of the model in different ways, and then I use best gamma from previous search, keep it constant and determine the other two arguments. Is this likely to give me a good result, or are all 3 parameters strongly depend on each other?


title: Correct Poisson Distribution Formulation

text: I have a Poisson formulation with the following parameters:   * I see 90 people in a city an hour   * 1.2% of people in the city are Indian How would I calculate:   * The maximum amount of time before I would see any Indian person with the highest probability   * The amount of Indian people I would see in an hour with the highest probability


title: Crazy divergence in mixed model - chains start well but go crazy later

text: I'm trying to compute a mixed model using jags (R2jags) and I got very strange divergence. The chains started so well, very well according to the results expected (also see `lmer` output of the same model below). But at certain point, the chains went crazy. Just look at the traceplot for **delta_tau** variable - the chains start so well, but then the green chain goes crazy, then blue and finally red... Any ideas why this happens? Can't be in initial values, because the chains started so well. Maybe the priors? Why is the system unstable? ![enter image description here](http://i.stack.imgur.com/rn08d.png) ![enter image description here](http://i.stack.imgur.com/v2jof.png) **EDIT:** Variables `gamma_tau` and `delta_tau` don't fall to exact zero, as you can see on these zoomed-in figures: ![enter image description here](http://i.stack.imgur.com/qkORJ.png) ![enter image description here](http://i.stack.imgur.com/gdEpz.png) This is the jags model:               model {          # likelihood     for (i in 1:N) {         logInd[i] ~ dnorm(mu[i], eps_tau)         mu[i] <- alpha[crit[i]] + (beta[crit[i]] + delta[species[i]])*year[i] + gamma[species[i]] # ekviv mix1b/c podle me     }          # priors     eps_tau ~ dgamma(1.0E-3, 1.0E-3)           for (j in 1:no_crit) {         alpha[j] ~ dnorm(0, 0.0001)         beta[j] ~ dnorm(0, 0.0001)      }          for (k in 1:no_species) {         gamma[k] ~ dnorm(0, gamma_tau)         delta[k] ~ dnorm(0, delta_tau)     }          gamma_tau ~ dgamma(1.0E-3, 1.0E-3)      delta_tau ~ dgamma(1.0E-3, 1.0E-3)     }      Code used to run jags (using R2jags):               no_crit = length(levels(crit))          win.data = list(logInd = mydata$logInd, crit = (as.integer(crit)),           year = mydata$Year, species = (as.integer(mydata$Taxon)),          N = nrow(mydata), no_crit = no_crit, no_species = length(levels(mydata$Taxon))     )          inits = function () { list(         alpha = rnorm(no_crit, 0, 10000),         beta = rnorm(no_crit, 0, 10000)     )}            params = c("alpha", "beta", "eps_tau", "gamma_tau", "delta_tau")          # ni: 1000 -> .. sec     ni <- 20000     nt <- 8     nb <- 8000     nc <- 3          out <- R2jags::jags(win.data, inits, params, "model.txt",         nc, ni, nb, nt,           working.directory = paste(getwd(), "/tmp_bugs/", sep = "")     )     R2jags::traceplot(out, mfrow = c(4, 2))      Here is output from the equivalent `lmer` model:               > summary(lmer(logInd ~ 0 + crit_i + Year:crit_i + (1 + Year|Taxon), data = datai2))     Linear mixed model fit by REML      Formula: logInd ~ 0 + crit_i + Year:crit_i + (1 + Year | Taxon)         Data: datai2        AIC  BIC logLik deviance REMLdev      8558 8630  -4267     8495    8534     Random effects:      Groups   Name        Variance   Std.Dev.   Corr        Taxon    (Intercept) 1.1682e-12 1.0808e-06                      Year        5.3860e-07 7.3389e-04 0.000       Residual             8.7038e-01 9.3294e-01            Number of obs: 2987, groups: Taxon, 103          Fixed effects:                    Estimate Std. Error t value     crit_iA      29.0539403  8.8116915   3.297     crit_iF       0.1848404  6.0286726   0.031     crit_iU      12.3405800 10.3326242   1.194     crit_iW       5.3248537  9.7416915   0.547     crit_iA:Year -0.0122717  0.0044174  -2.778     crit_iF:Year  0.0022365  0.0030222   0.740     crit_iU:Year -0.0038701  0.0051799  -0.747     crit_iW:Year -0.0003054  0.0048836  -0.063          Correlation of Fixed Effects:                 crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y     crit_iF      0.000                                               crit_iU      0.000  0.000                                        crit_iW      0.000  0.000  0.000                                 crit_iA:Yer -0.999  0.000  0.000  0.000                          crit_iF:Yer  0.000 -0.999  0.000  0.000  0.000                   crit_iU:Yer  0.000  0.000 -0.999  0.000  0.000  0.000            crit_iW:Yer  0.000  0.000  0.000 -0.999  0.000  0.000  0.000      Thanks in advance!


title: Using computer simulations to better understand statistical concepts at the graduate level

text: Hi I'm taking a graduate course in Statistics and we've been covering Test statistics, and other concepts. However, I am often able to apply the formulas and develop a sort-of intuition on how stuff works but I am often left with a feeling that perhaps if I backed up my study with simulated experiments I will develop better intuition into the problems at hand. So, I have been thinking of writing simple simulations to better understand some of the concepts we discuss in class. Now I could use say Java to:   1. Produce a random population with a normal mean and standard deviation.    2. Then take a small sample and try to try to empirically calculate Type-I and Type-II errors. Now the questions I have are:   1. Is this a legitimate approach to develop intuition?   2. Is there software to do this (`SAS`?, `R`?)    3. is this a discipline in Statistics that deals with such programming: experimental statistics?, computational statistics? simulation?


title: Modelling dependence between random variables

text: Is there any comprehensive reference on (or introduction to) how people have tried to model non-independent random variables? I already know about mixing processes, which express in various ways according to various coefficients how "future" events depend on "past" events, but that's about it...


title: How to test if the slopes in the linear model are equal to a fixed value?

text: Suppose we have a simple linear regression model $Z = aX + bY$ and would like to test the null hypothesis $H_0: a=b=\frac{1}{2}$ against the general alternative. I think one can use the estimate of $\hat{a}$ and $SE(\hat{a})$ and further apply a $Z$-test to get the confidence interval around $\frac{1}{2}$. Is this ok? The other question is strongly related to this one. Suppose that we have a sample $\\{(x_1,y_1,z_1),\ldots ,(x_n,y_n,z_n) \\}$ and we compute $\chi^2$ statistics \begin{equation} \sum_{i=1}^n \frac{(z_i-\frac{x_i+y_i}{2})^2}{\frac{x_i+y_i}{2}}. \end{equation} Can these statistics be used to test the same null hypothesis?


title: Is Cox regression analysis including demographic covariates required for observational studies instead of Kaplan-Meier?

text: I have been reading quite a few papers and recently have started analyzing some data on a review that I am doing. I am trying to examine the effect of two different treatment options on a group of patients. Using Kaplan-Meier analysis, the p-value is statistically significant at 0.03. However I believe that since this is an observational study that various confounding factors may be influencing my survival graphs in each arm. I have done a univariate analysis (using Kaplan-Meier and log rank testing) for three factors age, weight and height. My univariate analysis only identifies a statistically significant difference in survival among my patients in the age group (categorical age >45 or <45). Now I am not sure what to do.....I believe that I should now fit a Cox proportional hazard model and choose 'age' as a variable to control for and see if this affects my survival curves for treatment modality......is this the right thing to do? Or should I analyze each of the variables....weight, height, etc., with a different Cox regression model. And if I do, how do I interpret that data? For example, when placing 'age' in the variable column (I am using epi info) my p-value now is no longer significant in my analysis of survival between treatment modalities. Does that mean that age was one of the factors that caused my two groups to have differing survival times? Is this a multivariate analysis?


title: How many data points do we need for mixed effects longitudinal data?

text: I am collecting longitudinal data using for 4 time waves. Although the survey is administrated to the same population, different individuals may decide to complete it at each time point. As a result there are a number of individuals that only completed it once, others that completed it twice, some that completed it three times and others who participated in all four waves. For example at the moment there are about 2000 participants in time 1 and 1900 for time 2 but only 1200 that participated in both time 1 and time 2 (at the moment i am still collecting data for time 3 so i don't know yet what the final matched sample will be). The data are from different organizations so I would like to model this with using mixed effects with the lmer in R. e.g. ` lmer(outcome~"some repeated variables"+"organization level variables"+timewave+(timewave|subject)+(1|organizations)) ` My questions are * Do i need to remove individuals who completed it only once or twice to use a random slope for time?  * Is it meaningful to also try to fit a quadratic effect for time given that there are only 4 waves? (and would i need to remove subjects that have not participated in all four?)  Many thanks, George


title: Model for population density estimation

text: A database of (population, area, shape) can be used to map population density by assigning a constant value of population/area to each shape (which is a polygon such as a Census block, tract, county, state, whatever). Populations are usually not uniformly distributed within their polygons, however. Dasymetric mapping is the process of refining these density estimates by means of auxiliary data. It is an important problem in the social sciences as this recent review indicates. Suppose, then, that we have available an auxiliary map of land cover (or any other discrete factor). In the simplest case we can use obviously uninhabitable areas like waterbodies to delineate where the population isn't and, accordingly, assign all the population to the remaining areas. More generally, each Census unit $j$ is carved into $k$ portions having surface areas $x_{ji}$, $i = 1, 2, \ldots, k$. Our dataset is thereby augmented to a list of tuples $$(y_{j}, x_{j1}, x_{j2}, \ldots, x_{jk})$$ where $y_{j}$ is the population (assumed measured without error) in unit $j$ and--although this is not strictly the case--we may assume every $x_{ji}$ is also exactly measured. In these terms, the objective is to partition each $y_{j}$ into a sum $$ y_j = z_{j1} + z_{j2} + \cdots + z_{jk} $$ where each $z_{ji} \ge 0$ and $z_{ji}$ estimates the population within unit $j$ residing in land cover class $i$. The estimates need to be unbiased. This partition refines the population density map by assigning the density $z_{ji}/x_{ji}$ to the intersection of the $j^{\text{th}}$ Census polygon and the $i^{\text{th}}$ land cover class. This problem differs from standard regression settings in salient ways:   1. The partitioning of each $y_{j}$ must be exact.    2. The components of every partition must be non-negative.    3. There is (by assumption) no error in any of the data: all population counts $y_{j}$ and all areas $x_{ji}$ are correct.  There are many approaches to a solution, such as the "intelligent dasymetric mapping" method, but all those I have read about have _ad hoc_ elements and an obvious potential for bias. I am seeking answers that suggest creative, computationally tractable statistical methods. The immediate application concerns a collection of _c._ $10^{5}$ - $10^{6}$ Census units averaging 40 people apiece (although a sizable fraction have 0 people) and about a dozen land cover classes.


title: Sum of squares two ways, how are they connected?

text: The standard sum of squares as I know it is: $$ \sum(X-m)^2 $$ where $m$ is the mean. I ran into a different one which can be written two ways: $$ \sum(X^2) - \frac{(\sum X)^2}{n} = \sum(X^2) - m\sum X $$ I believe the latter is called the "correction term for the mean" (e.g. here). My algebra seems to be inadequate to show these are equivalent, so I was looking for a derivation.


title: Mean centering for PCA in a 2D array...across rows or cols?

text: I'm pretty new at this and I'm picking my way through the steps for running PCA on a 2D numpy array. Each subarray represents all pixels of an image (all rows & cols flattened). Example:               a = np.array([ [1,2,3], [4,5,6], [7,8,9] ])     # so, a[0], a[1], a[2] each represent a separate image      I'll use numpy.cov to calculate the covariance across these subarrays but first I need to mean-center the data, and that's where I'm getting confused. My novice question is: Should mean centering occur within each subarray? That is, should I calculate the mean of [1,2,3] and subtract it from each element, resulting in [-1,0,1], and then do the same to the next two subarrays (i.e., each subarray would get its own mean subtracted from each element)? Or, should mean centering occur across arrays? If so, across rows or cols? I've seen examples online of mean centering by calculating the mean along axis=0 (rows) (e.g., http://www.janeriksolem.net/2009/01/pca-for-images-using- python.html) and axis=1 (cols) (e.g., http://glowingpython.blogspot.it/2011/07/pca-and-image-compression-with- numpy.html). But I honestly don't know which is appropriate in this case.               a:     [[1 2 3]      [4 5 6]      [7 8 9]]          np.mean(a):     5.0          np.mean(a, axis=0):     [ 4.  5.  6.]          np.mean(a, axis=1):     [ 2.  5.  8.]          # which of the following mean-centered results makes sense?          a - np.mean(a):     [[-4. -3. -2.]      [-1.  0.  1.]      [ 2.  3.  4.]]          a - np.mean(a, axis=0):     [[-3. -3. -3.]      [ 0.  0.  0.]      [ 3.  3.  3.]]          a - np.mean(a, axis=1):     [[-1. -3. -5.]      [ 2.  0. -2.]      [ 5.  3.  1.]]


title: Regarding the convergence assessment on Markov chain monte carlo (MCMC)

text: During studying the convergence assessment on Markov Chain Monte Carlo, I once read the following statement: > A slowly converging sampler may be indistinguishable from one that will > never converge(e.g., due to nonidentifiability)! How to understand this statement? Thanks!


title: How to test a difference in proportions for significance when the underlying distributions are not normal?

text: I have a list of categories A, B, C, etc, each with a probability:               A   .015     B   .005     C   .02     D    ... etc.       I need to find if the difference between two of these categories, say               diff(C,A) = .02 - .015 = 0.005      is statistically significant. I have been using a Standard Error of Differences with the probabilities but since I don't have a normal distribution (and obviously not even continuous), what test would be recommended? Edit: The source of the data is essentially a list of the categories, these are not experiments I am performing but data that is given to me. The size of the list can vary but let's assume it's reasonable small (under 500). UPDATE: Accepted the answer below but you should read the comments attached to the answer for full understanding of the resolution.


title: Understanding the confidence band from a polynomial regression

text: I'm trying to understand the result I see in my graph below. Usually, I tend to use Excel and get a linear-regression line but in the case below I'm using R and I get a polynomial regression with the command:               ggplot(visual1, aes(ISSUE_DATE,COUNTED)) + geom_point() + geom_smooth()      So my questions boil down to this:   1. What is the gray area (arrow #1) around the blue regression line? Is this the standard deviation of the polynomial regression?    2. Can I say that the whatever is outside the gray area (arrow #2) is an 'outlier' and whatever falls inside the gray area (arrow #3) is within the standard deviation?  ![enter image description here](http://i.stack.imgur.com/XRurN.png)


title: Predicting percent of male and female party attendees

text: # Problem #1: easy, probably solved Suppose someone holds a party (A) and invites $n_A$ people. You're not attending a party, but you get full list of guests' identity card and are told that out of $n_A$ people $m_A$ were male (obviously, $n_A-m_A$ were female). After some time, there's another party (B) coming up. You also get a list of guests' identity cards ($n_B$ people in the list) and this time list includes $k$ people from the previous list, and $n_B - k$ new people, newer encountered before. **Problem** : Is it possible to predict the number of males and females attending the second party and give some reliability estimates? I believe that this should be a classic problem - but I haven't seen statistics for a decade and my skills became quite rusty. I've dug up some of my stat books for a while and couldn't find a reference. My thoughts:   * Edge case 1: if $k=0$, then it shouldn't be possible to predict anything   * Edge case 2: If $k = n_A = n_B$, then it's exactly $m_A$ males and $n_A - m_A$ females, with confidence = 1   * From what we know from the party A, for every guest who attended it, there's a $1/m_A$ probability that this guest's male and $1-1/m_A$ probability that this guest's female.   * It's easy to predict mean value of males for $k$ people, that would be $k/m_A$ and mean value of females - $k-k/m_A$.   * As for new guests, never seen before, it all depends on expected probability distribution. I see 2 viable alternatives:     * Assuming that there's an even distribution of males/females, $n_B-k$ new guests would be distributed 1:1 male:female, thus yielding mean values of $(n_B-k)/2$ males and $(n_B-k)/2$ females.     * Assuming the same distribution as it was at the party A, there would be mean value of $(n_B-k)(m_A/n_A)$ males and $(n_B-k)(n_A-m_A)/n_A$ females.   * Summing up, it gives us mean values:     * Version 1: $(n_B - k)/2 + k/m_A$ males and $(n_B + k)/2 - k/m_A$ females.     * Version 2: $(n_B-k)(m_A/n_A) + k/m_A$ males and $n_B - (n_B-k)(m_A/n_A) - k/m_A$ females Is there an error in my thoughts here? # Problem #2: harder, don't know how to solve If the solution to problem #1 seems okay, then a harder problem is when we have more than 1 party with known guest IDs and male/female distribution. We have full lists of guest and thus we can calculate any intersections between sets of guests at any parties. Probably, I should somehow adjust first probability estimation for any given guest (i.e. $1/m_a$ male probability) with additional data, but how to do that? Probably, it's also a classic problem, so if anyway can give me a pointer to the right direction, it would be most appreciated.


title: How to use Viterbi when there can be only one observation for a state

text: I want to calculate the best state sequence for a observation sequence. To do that, I want to use Viterbi algorithm. In my problem, there are two properties:   1. There can be only one observation possibility for a state. However the opposite is not true : an observation can be generated by more than one states.   2. States have a certain order. e.g. state Y can follow state X, but state X cannot follow state Y. More formally with a very simple case:   * States : $Q=q_{a1},q_{a2},q_{b1}, q_{b2}$   * State transitions : $A=[P(q_{b1}|q_{a1})=k, P(q_{b2}|q_{a1})=l, P(q_{b1}|q_{a2})=m, P(q_{b2}|q_{a2})=n]$, the rest ($P(q_{a1}|q_{a1})$, $P(q_{a1}|q_{b1})$, etc) is $0$ because of property (2)   * $k+l=1$, $m+n=1$   * Observations : $O=o_{a},o_{b}$   * where $length(Q) > length(O)$   * $P(o_{a}|q_{a1})=1, P(o_{a}|q_{a2})=1, P(o_{b}|q_{b1})=1, P(o_{b}|q_{b1})=1$ because of property (1)   * Observation likelihoods = $Q(q_{a1}|o_{a})=f, Q(q_{a2}|o_{a})=g, Q(q_{b1}|o_{b})=h, Q(q_{b2}|o_{b})=j$   * $f+g=1$, $h+j=1$ So,   * For example, for $q_{a1}$, observation is always $o_{a}$.   * But $o_{a}$ can be generated from $q_{a1}$ or $q_{a2}$. The question is : what is the most likely state sequence for given observation sequence $o_{a} o_{b}$ ? $q_{a1}q_{b1}$ or $q_{a1}q_{b2}$ or $q_{a2}q_{b1}$ or $q_{a2}q_{b2}$ To solve this problem I read that Viterbi is the way to go. However, I have a more specific model and applying Viterbi for this model is not really efficient. A lot of 0's and 1's in the calculation. Any directions for a better model / algorithm?


title: Dynamic conditional correlation in Stata

text: I am testing the co-movement between 2 exchanges by using the dynamic conditional correlation (DCC) developed by Robert Engle (2002). I want to apply this method in stata 12 and used this command:               Mgarch DCC (var1 var2=), arch(1) garch(1) distribution(t)      I read that it should give me a column with correlations per time unit. But i am not getting this colums. Does someone have any idea how to apply the DCC in Stata 12?


title: Model for predicting number of Youtube views of Gangnam Style

text: PSY's music video "Gangnam style" is popular, after a little more than 2 months it has about 540 million viewers. I learned this from my preteen children at dinner last week and soon the discussion went in the direction of if it was possible to do some kind of prediction of how many viewers there will be in 10-12 days and when(/if) the song will pass 800 million viewers or 1 billion viewers. Here is the picture from number of viewers since it was posted: ![PSY OGS](http://i.stack.imgur.com/tB74k.jpg) Here are the picture from number of viewers of the No1 "Justin Biever-Baby"and No2 "Eminem - Love the way you lie" music videos that both have been around for a much longer time ![Justin](http://i.stack.imgur.com/6mww9.jpg) ![Eminem](http://i.stack.imgur.com/xSUNA.jpg) My first attempt to reason about the model was that is should be a S-curve but this doesn't seem to fit the the No1 and No2 songs and also doesn't fit that there are no limit on how many views that the music video can have, only a slower growth. So my question is: what kind of model should I use to predict number of viewers of the music video?


title: Monte Carlo calculation of value at risk

text: I want to calculate the VaR with Monte Carlo Simulation, I am referring to this page: http://financetrain.com/calculating-var-using-monte-carlo-simulation/ which was created to the JP Morgan desrciption: http://www.jpmorgan.com/tss/General/Risk_Management/1159380637650 Now I have several questions, where I could not find any answer: The "random variable $\varphi$ generated from a normal distribution": Why is it from a normal distribution? I thought it should be from a uniform variable which is then transformed by plugging it in into the standard stock price model with the corresponding sample mean and sample variance? Standard stock price model: $R_i=\frac{(S_{i+1}-S_i)}{S_i}=\mu \Delta t + \sigma \varphi \sqrt{\Delta t}$ The Monte Carlo Simulation has mathematical connection to Kolmogorovs Strong Law of Large Numbers. Could you explain me this in this context? Is the following correct: I simulate 1000 final prices, the estimator of the mean (sample average of all final prices) will converge to the real unknown mean, this is K SLLN right? The final stock prices will be normal distributed, but due to the fact that the random number is from a normal distribution and the stock price model function is just a linear transformation or because of the central limit theorem? If this is not correct, what does Monte Carlo has to do with the Central Limit Theorem? Also I have at the end approx. 1000 stock pathes with the final stock prices. I am going to plot those and is it ok, if I "underlay" a normal distribution graphically at the end of the plot for the final stock prices? So these should be normal distributed right? Again: Is this a consequence of K SLLN or due to the random number generation from a normal distribution? A lot of questions also a lot of not understanding, so thanks for your answers if you answer.


title: Is confirmatory vs exploratory statistics "induction vs deduction"?

text: This webpage says: > Inferential Statistics - Deductive Approach > > Descriptive Statistics - Inductive Approach But I doubt it. If I understand correctly,   * Inferential Statistics is "given some data, find the probability model that generates the data", so it is a specific-to-general logic process, and therefore it is induction.   * Descriptive Statistics is "given some data, represent the data in another way", so it is specific-to-specific process of changing data's representation. So it is neither induction or deduction. Am I wrong?


title: Multivariate linear regression vs neural network?

text: It seems that it is possible to get similar results to a neural network with a multivariate linear regression in some cases, and multivariate linear regression is super fast and easy. Under what circumstances can neural networks give better results than multivariate linear regression?


title: Since the beta distribution is similar in form to the binomial, why do we need the beta distribution?

text: It appears that the binomial distribution is very similar in form to the beta distribution and that I can re-parametrize constants on either pdf to make them look the same. So, why do we need the beta distribution? Is it for a specific purpose? Thanks!


title: within subjects design with condition rotated across presentation lists (staggered cross-over design?)

text: In my field, it is common to rotate levels of an experimental condition across different presentation lists in a within-subjects design. This is done to prevent a particular participant from seeing every level of every experimental item, while at the same time, allowing them to still see equal instances of each type of condition level. Presumably, this reduces the chance of there being learning effects etc. To illustrate - say I have an independent variable, "Condition", with three levels (A, B, C). I also have six experimental Items (1-6), which appears in three versions, in accordance with the four condition levels (note: I actually use many more items in the real experiment):               1A     1B     1C     2A     2B     2C     3A     3B     3C     4A     4B     4C     5A     5B     5C     6A     6B     6C      I then distribute the levels of each item across three different presentation lists (three versions of the experiment) with a simple latin square:               List1     List2     List3     1A     1B     1C     2B     2C     2A     3C     3A     3B     4A     4B     4C     5B     5C     5A     6C     6A     6B      Each participant only sees one list - only one version of the experiment (List1, 2, or 3). Thus the design is "within-subjects", but each participant/subject will only see one level for each item (here, a total of two instances of each level of Condition). I used to analyze such data using repeated measures ANOVA. For this, I would generate Participant means for each Condition. Each mean would be based off of the same number of cells. I also used to do the same for Item (collapsing Participant information). However, now I am trying to use linear mixed models in R. However, there is no aggregation for means with this approach, and there will thus be a good number of empty cells in the output. I have very little understanding as to how this statistical approach is treating my data, and what some potential drawbacks (or flat out violations of statistical assumptions) might be. Some specific questions: 1) Can lmer accomodate this somewhat regulated "missing data" that is brought on by my list distribution, or will the analysis somehow be critically compromised? 2) How critical is it to have equal numbers of participants across all participant lists (with ANOVA, this was important, but I am not sure now with mixed models)? I would appreciate any comments on this. I have had a hard time finding mention of this type of within-subjects design in the stats literature. Perhaps this is because I don't know what to really call this type of design. This post is somewhat similar: Is there literature on staggered within subjects designs? What are the consequences of such a design? but it has a different end question, as I am interested in how lmer can handle this. Also, I am not entirely convinced that the design I am discussing is a "crossover" design. Since I have the added dimension of there being a number of experimental "Items" (in my case, these are sentences). Thanks


title: Definition of statistical model in case of hierarchical model

text: In wikipedia the definition of a parametric model is the following: > A parametric model is a collection of distributions, each of which is > indexed by a unique finite-dimensional parameter: > $\mathcal{P}=\\{\mathbb{P}_{\theta} : \theta \in \Theta\\}$, where $\theta$ > is a parameter and $\Theta \subseteq \mathbb{R}^d$ is the feasible region of > parameters, which is a subset of d-dimensional. I wonder what is $\Theta$ in the case of a hierarchical model. Is it composed of all the latent variables of the model or only the one at the top level? Does this include the hyper-parameters? My concern is that this definition has an influence on the definition of model identifiability.


title: Correlation between the proportion derived from the data of different size

text: I want to obtain a correlation between the accuracy of two linguistic features by language learners. The learners wrote a number of essays, and the number of correct and incorrect use of the article and the plural marker was tallied in each essay. The data for a learner look like the following, though the number of essays varies across learners.               Essay   ArticleCorrect   ArticleIncorrect   PluralCorrect   PluralIncorrect       1          7                0                  2               1       2          4                1                  0               2       3          6                2                  3               1       4          7                0                  2               1       5          8                2                  1               3      ...        ...              ...                ...             ...      10          4                1                  2               0      What I would like to know is whether the performance of the learner on the article use is correlated with his/her performance on the plural marker (e.g., whether the accuracy of the plural marker tends to be high when the accuracy of the article is high). Straightforwardly applying the Pearson correlation on the accuracy (i.e., $\frac{correct}{correct + incorrect}$ ) is not adequate as the number of attempts varies across essays as well as across the linguistic features. Weighted correlation, I believe, only allows weighting by essays but not by linguistic features. Is there a way (preferably in R) to do the above? Any help is appreciated.


title: Is there a test for comparing principal components decompositions between samples?

text: Is there any test for comparing the directions, magnitudes, etc of PCA results for different samples? I'm leaving the nature of the test deliberately vague because I'd like to hear all the various possibilities... e.g. there might be (and I'm speculating here) a test comparing the sizes of the first principal components, or a test comparing the directions of the principal components, or there's some kind of distance measure between PCA outcomes and a test statistic for their equality. As far as a use case goes, I don't have on in mind. Just out of curiosity, maybe as an exploratory technique.


title: How to test for presence of trend in time series?

text: Apart from detecting trend from a time series plot, how do you test for its presence before removing the trend using moving average? I fitted a mathematical trend to the data and the slope was approximately zero, does that mean there is no trend?


title: How to transform leptokurtic distribution to normality?

text: Suppose I have a leptokurtic variable that I would like to transform to normality. What transformations can accomplish this task? I am well aware that transforming data may not always be desirable, but as an academic pursuit, suppose I want to "hammer" the data into normality. Additionally, as you can tell from the plot, all values are strictly positive. I have tried a variety of transformations (pretty much anything I have seen used before, including $\frac 1 X,\sqrt X,\text{asinh}(X)$, etc.), but none of them work particularly well. Are there well-known transformations for making leptokurtic distributions more nearly normal? See the example normal Q-Q plot below: ![enter image description here](http://i.stack.imgur.com/IgSjh.png)


title: Heat kernel and laplacian

text: Is there any relation between heat kernel and Laplacian? I know what each of them is but I am not sure about the relation between them.


title: Interpreting plm output in R - number of observations used with very unbalanced panel

text: I'm running a fixed effects model with the `plm` function, and I'm looking for help interpreting an aspect of the output. If the output reads:               Call:     plm(formula = dependentvariable ~ independentvars, data = data, model = "within",          type="time")          Unbalanced Panel: n=176, T=1-2, N=211          Residuals :        Min. 1st Qu.  Median 3rd Qu.    Max.      -0.0654  0.0000  0.0000  0.0000  0.0654           Coefficients :         Estimate Std. Error t-value Pr(>|t|)        x1  -0.4219101  0.1662230 -2.5382 0.020054 *      x2  -0.0072536  0.0069678 -1.0410 0.310933        x3  -0.2221514  0.0574869 -3.8644 0.001044 **     x4   0.1118861  0.1247960  0.8966 0.381177             ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Total Sum of Squares:    0.18087     Residual Sum of Squares: 0.045841     R-Squared      :  0.74656            Adj. R-Squared :  0.067226      F-statistic: 5.08807 on 11 and 19 DF, p-value: 0.00099507      How do I know how many observations were used here, since my panel is extremely unbalanced? With both a small n and large N provided in the output, I'm confused about how to find how many observations were included.


title: How does maximum likelihood estimation have an approximate normal distribution?

text: I have been reading about MLE as a method of generating a fitted distribution. I came across a statement saying that maximum likelihood estimates "have approximate normal distributions." Does this mean that if I apply MLE repeated times over my data and the family of distributions that I'm attempting to fit to, the models I get will be normally distributed? How exactly does a sequence of distributions have a distribution?


title: Comparing 2 independent counts of the same data

text: I am looking for an appropriate statistical method to test my hypothesis that there is no significant difference between counts conducted by 2 individuals of the same data. The data is as follows: I have videos of fish migrating up a river, and 2 different people counting the videos. I have 960 videos, so each video is counted once by each person. I am unsure what type of statistical test will evaluate the difference between counters, as this is ultimately what I am testing: How well do counters agree with one another, in other words, is one counter significantly different than the other counter when presented with the same data. If anyone could suggest the appropriate statistical test that would be awesome. Thanks!


title: Quantitative research comparing perception of preparation to evaluation

text: I am trying to complete my PHD and I am using a survey to determine students' perception of their preparedness to teach and want to compare their answers to what their final summative evaluation says about their ability to teach. I created a questionnaire for the students to answer that reflects each of the areas of the summative evaluation. What statistical test can I use to compare the results? I also want to do another analysis and break the data down into the three different category of schools the students teach at to see if that makes a different in their perception of preparation. I did a field test to establish a level of construct validity of the questionnaire would that suffice? I only have 114 for the study and don't think breaking the group into 2 groups would give me strong enough numbers for the dissertation study? Thank you for your assistance! Jane


title: Probability in hidden Markov model

text: I'm considering a Hidden Markov Model as follows: $X_{n+1} = F_n(X_n,\Theta)$ $Y_n = G_n(X_n,\Theta)$ where the $Y_n$ are the observations, the $X_n$ the hidden states and $\Theta$ the parameters. At some point, I have to deal with this probability $p(X_k = x_k | X_{k-1} = x_{k-1}, Y_{0:N} = y_{0:N})$ where $0 < k < N$ and $Y_{0:N}$ denotes the $Y_0, \dots, Y_N$. If it were $p(X_k = x_k | X_{k-1} = x_{k-1}, Y_{0:k} = y_{0:k})$ I would know how to deal with it (I think) but I wonder what it changes to add the observations for all future times. Any idea how I could handle this probability, maybe express it with respect to usual and known probabilities? Is there something I'm completely missing here? Thanks a lot,


title: Fixed effect/mixed model

text: Unit of analysis: States (39)   Time: 3 years   Total observations: 117   Dependent variable: % children graduating within 4 years from high school   Independent variable: 8 composite scores ranging from 1-30. These scores did NOT change throughout the period of examination (i.e. 8 measures are time- invariant).   I am trying to determine whether to proceed as a fixed effect or mixed effect model. Does the fact that I have only 39 states (and not the population of all 50) mean I have to go mixed effect? A followup: if it is mixed, I assume 1) my fixed variable would be the IVs, and 2) my random would be state? Thanks!


title: An alternative form of $L$-estimators

text: $L$-estimators are based on the ordered observations $X_{(1)} \leq X_{(2)} \leq \ldots \leq X_{(n)}$ of the random sample $X_1, X_2, \ldots , X_n$. The general $L$-estimtor can be written in the form: $$ T_n = \sum_{i=1}^{n}{c_{ni} h(X_{(i)})} + \sum_{j=1}^{k}{a_j h^{*}(X_{( [n p_{j} ] +1)})} $$ where $0< p_1 \ldots < p_k < 1$ and $[ \cdot]$ represents greatest integer function. $h( \cdot )$ and $h^{*}(\cdot)$ are given functions, coefficients $c_{ni}$ ($i=1,2, \ldots, n$) are generated by a bounded weight function $J: [0,1] \longrightarrow \mathbb{R}$ in the following way, either: $$ c_{ni} = \int_{\frac{i-1}{n}}^{\frac{i}{n}}{J(s)}ds, i=1,2,\ldots,n $$ or Approximately, $$ c_{ni}= \frac{1}{n} J(\frac{i}{n+1}), i=1,2,\ldots,n $$   1. Show that $T_n$ can be expressed in an alternative way $$T_n= \int_{0}^{1}{J(s)h(Q_n(s))}ds,$$ where $Q_n$ can be defined as either $Q_n(s)=\\{X_{(i)} , \ \ if \ \frac{i-1}{n} < s \leq \frac{i}{n}$, $i=1,2,\ldots,n$} and $\int_{0}^{1}{J(u)}du=1$ or $Q_n(s) = \inf_x \\{x : F_n(x) \geq s \\}$ where $F_n$ is the empirical distribution function of $F$.   2. Show that influence function of $T(F)$ at $x$ i.e. $$IF(x,T;F)= \int_{- \infty}^{\infty}{F(y)h'(y)J(F(y))}dy - \int_{x}^{\infty} {h'(y)J(F(y))}dy$$ where $T(F)=\int_{0}^{1}{J(s)h(Q_s))}ds$


title: PCA, randomness of component?

text: I'm trying to decide if a component of a PCA shall be retained, or not. There are a gazillion of criteria based on the magnitude of the eigenvalue, described and compared e.g. here or here. However, in my application I know that the small(est) eigenvalue will be small compared to the large(st) eigenvalue and the criteria based on magnitude would all reject the small(est) one. This is not what I want. What I am interested in: is there any method known that takes the actual corresponding component of the small eigenvalue into account, in the sense: is it really "just" noise as implied in all the textbooks, or is there "something" of potential interest left? If it is really noise, remove it, otherwise keep it, regardless of the magnitude of the eigenvalue. Is there some kind of established randomness or distribution test for components in PCA that I am unable to find? Or does anyone know of a reason that this would be a silly idea? Thanks! Addition: histograms (green) and normal approximations (blue) of components in two use cases: once probably really noise, once probably not "just" noise (yes, the values are small, but probably not random). The largest singular value is ~160 in both cases, the smallest, i.e. this singular value, is 0.0xx - way too small for any of the cut-off methods. What I'm looking for is a way to formalize this ... ![probably really "just" noise](http://i.stack.imgur.com/2Zcq1.png) ![probably not noise but may contain interesting bits](http://i.stack.imgur.com/HnmGA.png)


title: Confusion over lmer and p-values

text: I was under the impression that the function `lmer()` in the `lme4` package didn't produce p-values (see `lmer`, p-values and all that). I've been using MCMC generated p values instead as per this question: Significant effect in `lme4` mixed model and this question: Can't find p-values in the output from `lmer()` in the `lm4` package in `R`. Recently I tried a package called memisc and its `getSummary.mer()` to get the fixed effects of my model into a csv file. As if by magic, a column called p appears which matches my MCMC p values extremely closely (and doesn't suffer the processing time that comes with using `pvals.fnc()`). I've tentatively had a look at the code in `getSummary.mer` and have spotted the line that generates the p value: `p <- (1 - pnorm(abs(smry@coefs[, 3]))) * 2` Does this mean p values can be generated directly from `lmer`'s output rather than running `pvals.fnc`? I realise this will no doubt start the 'p value fetishism' debate but I'm interested to know. I've not heard `memisc` mentioned before when it comes to `lmer`. To be more succinct: What is the benefit (if any) of using MCMC p-values over those generated by `getSummary.mer()`?


title: Logit link notation

text: Here is the whole question: > (a)A Simulation Study (Logistic Regression). Assume $y|x \sim > \mathrm{Binary}(p)$, where $p= \mathrm{E}(y|x)$, and > $\mathrm{logit}(p_i)=-1+5.1x_{1i}-0.3*x_{2i}$. Generate data with > $x_{1i}\sim \mathrm{Unif}(0,1)$, $x_{2i}=1$ for $i$ odd and $x_{2i}=0$ for > $i$ even, and sample size $n$=500. Try generalized linear model (GLM) with > logistic and probit links. What’s your finding? We are required to use R to do it. I am not familiar with the notations like $y|x\sim\mathrm{Binary}(p)$,where $p= \mathrm{E}(y|x)$, and $\mathrm{logit}(p_i)=-1+5.1x_{1i}-0.3*x_{2i}$, with $x_{1i}~\mathrm{Unif}(0,1)$ ,$x_{2i}=1$ for $i$ odd and $x_{2i}=0$ for $i$ even,and sample size $n$=500. Also how to use to generate $x_1$ and $x_2$? Many thanks.


title: How do I handle a normal distribution with peaks at each end?

text: This question may appear to be quite odd, but the following explanation should make it a little more comprehensible. I work on the analysis of population trends of birds. In our team, we work with a simple trend index, that is calculated as $$\frac{\rho_t - \rho_{t+1}}{\rho_t + \rho_{t+1}} = \frac{(\mbox{density at time }t) - (\mbox{density at time }t+1)}{(\mbox{density at time t })+ (\mbox{density at time }t+1)}$$ When the population size does not change, the index value is $0$. When the population decreases, the index is negative; when the species goes extinct, the index value is $-1$. An increasing population yields positive index values, with the extreme of $1$ when the species colonises a formerly uninhabited area. We calculate this index for a large number of monitoring plots. For abundant species the index values of all plots together yield a nice, more or less bell-shaped curve that is convenient for further analysis. For less abundant species, however, the number of plots where it disappeared, and the number of newly colonised plots can be large leading to a three-modal curve with modes at $-1$, $0$ and $1$. The rarer the species, the higher get the peaks at $-1$ (extinction) and $+1$ (colonisation) and the flatter gets the bump in the middle. Analyses of such data (e.g. in various types of regression analysis) are difficult since I am not aware of a statistical distribution that can describe them. I know that this is not a very precise question that can easily be answered, but I'd be thankful for any advice how to treat data like this or how to calculate a more "user-friendly" trend index for two points in time.


title: Can someone help me read this hierarchical regression table?

text: What do all the symbols mean? Thanks for your help!               Table 2     Summary of Hierarchical Regression Analysis for Variables Predicting Depression (N=108)                           Model 1                 Model 2     Variable          B     SE B    β          B    SE B    β     RSES            -.613   .119  -.449     -.600   .120   -.440     UES             -.027   .044  -.055     R2                  .201                    .194     Adjusted R2         .194                    .189      F                   26.705***               13.473***      F change            26.705***                0.395      *** p<.001


title: The curse of dimensionality? (linear SVMs)

text: How do you know whether you suffer from it? Let's suppose I have a 2 class problem - 2000 training examples and 30 features. While it works good for the most part, sometimes I get edge cases that are so horribly misclassified. I feel like it's because there is a lot of dimensions and too little examples. I feel like the dimensionality is affecting my results - is there a way to check if it's true? Do I need to have several examples for each combination of the feature values? With 30 features from which each can have hundreds of values, this would require ENORMOUS amount of data... How much improvement can I anticipate if I were to increase the training set size to 4000 examples, 8000 examples, 16000 examples etc. - is there a rule of thumb when trying to figure this out?


title: How do you study for statistics? What methods do you use to learn?

text: I am a grad student. I am struggling with core statistics courses: probability, statistical inference. In probability, for example, prove unbiased estimator for uniform. In statistical inference, find the UMVUE for an estimator. American stat dept has either highly qualified foreign students, or the couple of Americans who are struggling. I would like to learn this and have a good foundation. Looking to pass the Ph.D. exam in the future. All I do is study, but still get 60/100 on my midterms. My problem is, the prof never tests on what he taught, but there are always students that passes, and continue on. I don't know what to do anymore. What are your suggestions?


title: Use of latent semantic analysis in topic modeling

text: I am learning Latent Semantic Analysis. After we have the SVD decomposition of term-document matrix, how to infer topics from that? For example, how to get the following result: Topic 1: apple, orange, fruit Topic 2: computer, machine


How do I show that Fisher's and Normal Linear decision rule differ by a constant? Thank you for your help.


title: Visualizing the difference between uncorrelatedness and orthogonality

text: A previous question received a linked answer as below: > The first important distinction here is that linear independence and > orthogonality are properties of the raw variables, while zero correlation is > a property of the centered variables and > The key to appreciating this distinction is recognizing that centering each > variable can and often will change the angle between the two vectors. I'm having a hard time visualizing what "centering a variable" means and how it may change the angle. Could anyone produce an image?


title: Dataset similarity measurement

text: Suppose: For each input data, I could obtain true result as well as experiment result from my algorithm, one example is The true result is `A = {a, b, c},` Experiment result is `B = {b (0.3), c (0.8)}`, the numbers denote the probability of the element's presence. The are many pairs of `(A, B)` Could you suggest any distance/similarity measurement between A and B so that I could show that my algorithm works well? My friend suggested: _L2 Norm_ , _T-test_ , _KL-divergence_ (the last one I think inappropriate) Thanks


title: How to handle multiple comparisons in a three-arm clinical trial?

text: I have 3 arms in a trial. I want to compare results of a survey before during and on completion of a treatment. Data is not normally distributed. What test should I use?


title: Cross-validation for log-transformed x's and y's

text: How to cross-validate a model with log-transformed x's and y's and obtain the RMSE in the original units? For example, I got an error message when I used the cvTools package:                library(cvTools)      library("robustbase")      data("coleman")      set.seed(1234)      folds <- cvFolds(nrow(coleman), K = 5, R = 10)      fitLm <- lm(log(Y) ~ ., data = coleman)      repCV(fitLm, cost = rtmspe, folds = folds, trim = 0.1)      Error in eval(expr, envir, enclos) : object 'Y' not found


title: Basic misunderstanding with the Augmented Dickey Fuller test

text: Suppose we have the following model: ![enter image description here](http://i.stack.imgur.com/ogehC.png) ADF would test if $\gamma$ equals 0 (null hypothesis), and rejects it if the statistic is sufficiently negative in favor of $\gamma < 0.$ But what if, say, gamma equals -5? Wouldn't this mean that yes, there's no unit root, but since $\gamma = -5$ the process still won't be stationary, it would explode? If so, is it really correct to use ADF to conclude that the process is stationary when we get a low enough negative value (which is what I was taught) as opposed to just concluding that the process has no unit root? I must be missing something. Thanks!