title: ANOVA with very uneven sample sizes-- Can fit be improved through reference to a more complete dataset?

text: I am analyzing a set of costs for contractor services for a particular city. These data have been extremely difficult to collect due the difficulty of contacting contractors and reluctance of contractors to divulge their pricing information. Factors are "service type" and "contractor". We generally have quotes from more than one contractor on each service and a number of contractors perform most of the services. For a few services we have fairly good replication (i.e. ~15 quotes), for some of the services we have some degree of replication (i.e. 5-10 quotes per service). But for about a third of the services replication is very poor, i.e. <5 quotes per service. The data have proven extremely challenging and time consuming to collect and I doubt if there will be more. The upside to this story is that there are several construction cost estimation tools that can give us estimates of the costs of these services. And while these costs wouldn't be specific to the area we're interested in, we think we can assume that the prices of the services relative to each other would be similar to what we expect to see in the area we're looking at. Sometimes elegance is the first casualty of pragmatism, but this is what I'm thinking of doing-- I am planning to collect a parallel set of cost estimates from these cost estimation tools and compare the relationships between the results of our surveys to the relationships of costs within this set of parallel cost data. Essentially using the better replicated data points in our dataset as anchors and generating expectations for our poorly replicated data points based on the relationships from the parallel dataset. I'm not even sure what this approach would be called and it feels like a very specific type of issue, so it's difficult to research online. But I'm wondering if there's any theoretical basis for this type of approach that's been worked out. Shall I buy a book on Bayesian stats? Most of my experience in grad school was analyzing data from designed experiments, so this sort of situation is interesting, but definitely foreign to me. Thanks in advance for your help.


title: Comparing means with unequal variance and very different sample sizes

text: I am trying to compare the means of the same variable between men and women. This is the statistics:                    N        Mean        Variance    Coef. Var.     Gender            2000      26.12         10.89         0.13         Male                  50      56.10         25.01         0.09        Female      Neither variable is normally distributed but taking the log makes it pretty darn close. What is the appropriate way to test the means between males and females? Should I use the log or not? Any additional advice using Stata would be helpful. My initial reaction is that females fare better than men, but I want to be statistically rigorous.


title: Mann-Whitney U test with very large sample size?

text: I'm doing a Mann-Whitney U test to compare two very large samples (sample size 1 = 13250; samlple size 2 = 38871) originating from a raster image. I know t-tests are not recommended to compare rasters, because since rasters have so many values, they will almost surely detect a significant difference, no matter how small that difference might be (see this post in gis.stackexchange.com and point 3 of the first answer to this post). **My question is whether the same problem applies to Mann-Whitney or not**. I ran the test on SPSS and got the following results:   * Test Statistics    * selection frequency   * Mann-Whitney U 4094520,000   * Wilcoxon W 759591276,000   * Z -210,227   * Asymp. Sig. (2-tailed) ,000 While I did expect to find a difference between the groups, I don't know what to make of such large values. There is a paper which did the same thing I did. See the last sentence of the Material and Methods section. However, they obtained much smaller values http://icesjms.oxfordjournals.org/content/69/1/75/T5.expansion.html (they don't make clear if these are Mann-Whitney Us or z-scores). Granted they used a smaller sample (combined sample size = 5629), but the difference in the magnitude of the values still seem strange to me. **So, are my results simply the result of a very large sample, but still valid? Or should I use another test?**


title: How do you compute the annual standard error of a regression model when the model itself is based on monthly observations?

text: Let's say you regress monthly stock returns, and the regression model has a standard error around such monthly stock returns of 2%. Next, you forecast the next year annual stock return based on the regressed monthly returns. What is the standard error around this annual stock return? Some may think this is not a very good example. Don't get hung up on what is the proper way to model stock returns. This has nothing to do with the question. The question is simply figuring out the calculation to convert a monthly standard error into an annual one.


title: What is the difference between Rf_dpois in Rmath.h and the dpois that I use directly in R?

text: If I see `Rmath.h` in `/usr/share/R/include`, the signature of the function `dpois` or `Rf_dpois` is               double dpois(double, double, int);       However, if I do `?dpois` in `R`, I see:               dpois(x, lambda, log = FALSE)      Are these both the same thing? If yes, can someone please clarify? Also I tried to find the code for `dpois` at R-svn. Can someone please tell how to find the code for functions like these? Thank you for any help or pointers.


title: Splitting one variable according to bins from another variable

text: I have continuous data "A", binary categorical data "O", gender/sex and age for several participants in a study. A linear model in R shows no correlation between A and age. I would now like to group A into groups by age and see if there is a difference between the groups. I know about 'hist' and 'split' in R, but these do not do what I need. (1) How can I divide/split A into groups based on age (18 to 27, 28 to 37, etc.). (2) Once I've done that, I would use a $\chi^2$ test? (3) Could I also test O in the same groups, using counts?


title: PCA Based Filtering but only filter out small values

text: PCA based filtering is used to identify and eliminate noise in data. This would basically involve computing the PCs and using the top k PCs to denoise the data. What if I know for sure that only the extremely small values in my matrix are noise? Now, a value may be small w.r.t the entire matrix but not small w.r.t a particular row/column. Can I achieve this with some transformation of the input matrix followed by a PCA based filter?


title: Sufficient statistics and parametric bootstrapping

text: Does the resampling step of the bootstrap method require to have the sample entirely, or a sufficient statistic suffices? In general, the bootstrapping can be nonparametric or parametric. In the former case, we can simply keep the frequency of each distinct observation, and we have the empirical distribution to draw bootstrap samples from. In the latter case, the idea of only keeping a _sufficient statistic_ seems interesting. Are there any pitfalls which should be avoided in this case?


title: Closed form posteriors for a simple bivariate Bayesian regression

text: I'm analyzing a simple linear regression $Y_{i}$~$a+b*X_{i}+e_{i}$, with $e$ being normally distributed with known variance and where I have normal priors on $a$ and $b$. I'm trying to piece together an approximate closed form solution for the variance of my posterior of $a$ as a function of sample size for the purposes of an optimization I'm trying to do. I know that in general, there's a closed form conjugate prior for Bayesian regression in the multivariate case in matrix form. So I can do some sort of moment matching and then calculate out the matrix algebra for the bivariate case and probably can end up with a nice looking closed form. I've tried that, but I keep getting bogged in algebra, and I feel like this has to be a well- known problem that's been fleshed out already. Does anybody have a citation or link to point me to?


title: Random permutation test for feature selection

text: I am confused about permutation analysis for feature selection in a logistic regression context.   Could you provide a clear explanation of the random permutation test and how does it applies to feature selection? Possibly with exact algorithm and examples. Finally, How does it compare to other shrinkage methods like Lasso or LAR?


title: Smart date parsing in R?

text: Is there a function that does smart date parsing in R? I know the `strftime`/`as.POSIXct`/`as.POSIXlt` commands, but they require a date formatting string, or throw the error "character string is not in a standard unambiguous format". This happens even when I pass a string with more than enough information to be parseable, like "Fri Apr 29 16:43:20 GMT 2011". It would be really nice to not bother with reverse-engineering the format every time I import dates from a new source. Is there code that already does this? Thanks!


title: What variation can I expect between lab results if the coefficient of variation is 9%?

text: I know that correlation of lab to lab test results for clotting times of patients on anticoagulation treatment is not great. Clotting times are measured in INR (international normalized ratio) and a research paper I am reading asserts > For the larger hospital laboratories analysing external controls with INR > > 2, the mean coefficient of variation (CV) was reduced from 9.1% to 5.6% I have a patient with an INR of 2.5, what range of INRs could I expect to occur by chance alone if I sent his samples to a range of different laboratorys if the COV was 9% ? paper on interlab COV The question arises because at my medical practice we use a near patient testing for INR testing (different methodology to the whole blood test the hospital lab uses). We test our equipment quarterly against blind samples sent to us by a national reference lab. For 5 years we have always been within the "in range" test limits, so far, our equipment has been reliable. A patient had had an INR test by a hospital lab about 10 days later than our test, our test result was 2.5, the hospital test was 1.7 and the patient was worried this meant that someone was "wrong". paper looking at coaguckek S vs hospital lab INR My instinct was that given the 10 day interval between tests, my experience that day to day variation of INR on the same patient can be 10%, the different assay techniques of hospital to near patient device, and what I had looked up about inter hospital COV , that this was probably chance alone. Out of curiosity, and to reassure my patient and to double check for safety I will do a near patient test and send a hospital sample from the patient on the same day. My gut feeling was that if the INRs were within 0.5 of each other that would be fine. However, even better than gut feeling, I thought I would ask the experts what range would likely represent chance variation , very happy to take 95% confidence :-)


title: Is it necessary to use cross-validatation to avoid overfitting when applying random forest algorithm?

text: > Is it necessary to use cross-validatation to avoid overfitting when applying > random forest algorithm? This is a question asked in my recent data scientist interview. Can anyone give any ideas? Thanks


title: Quantile regression analyzing the conditional quantiles of one of the regressors?

text: Given response $Y_t$ and predictor $X_t$, we can use OLS to analyze the conditional mean; $E[Y_t | X]$. Quantile regression can be used to analyze the conditional quantile function; $Q(Y_t(\tau)|X)$. Specifically; $$(1) \,\,\,\,\,\,\,\,\,\,\,Q(Y_t(\tau)|X) = a(\tau) + b(\tau)X_t + e_t$$ In this one variable case, we can easily look at the conditional quantile function of $X_t$ by simply estimating: $$(2) \,\,\,\,\,\,\,\,\,\,\,Q(X_t(\tau)|Y) = \alpha(\tau) + \beta(\tau)Y_t + e_t$$ However, consider the following specification: $$(3) \,\,\,\,\,\,\,\,\,\,\,Q(Y_t(\tau)|X,Z) = a(\tau) + b(\tau)X_t + c(\tau)Z_t + e_t$$ In this case, I can't see a way to do something similar to what I did in $(2)$. _My question is_ : **Can I do something similar to what I'm proposing in (4), where I analyze the response, $Y_t$, with respect to the predictor $Z_t$ as well as the _conditional quantiles_ of the predictor $X_t$?** $$(4) \,\,\,\,\,\,\,\,\,\,\,Y_t = \alpha(\tau) + \beta(\tau)Q(X_t(\tau)|Y,Z) + \gamma(\tau)Z_t + e_t$$ (excuse poor notation)


title: Spectral clustering of graph

text: I am trying to cluster the graph using spectral clustering. However I am unaware of the number of classes that exist in the data. Will it be a good idea to do PCA on the adjacency matrix to find number of actual clusters in the dataset? Are there any other options?


title: Test for difference of distributions on a torus

text: I have two circular dependent variables and would like to test for a difference in the distributions (presumably circular means) between multiple treatment groups. There are a number of multivariate tests such as MANOVA that will work on $DV \in \mathbb{R}^2$ but I'm looking for a test of differences of means (or something) on $DV \in \mathbb{S}\times\mathbb{S}$. If it can generalize to higher orders spaces ( _e.g._ , $\mathbb{S}\times\mathbb{S}\times...\times\mathbb{S}$) that would be a plus. The test employed should consider the multivariate interactions, so we can't use something like a Kuiper test from circular statistics on each DV by itself. We've considered transforming the data into some linear metric to use in a test, such as a Wilcoxon rank-sum on within vs. across condition cosines. However, it would be greatly preferred if we can cite an article describing the test rather than describe it ourselves. As an example dataset consider hatching fruit flies. We measure the time of hatching and the direction towards which they initially fly away. We then want to look for an effect of a drug applied to the eggs.


title: Study design involving 4 variables measured every month for 5 years

text: I want to design a study on hemodialysis patients. I know from basic science knowledge that there are four variables (suppose a, b, c, d) which effect the rates of heart attacks/ strokes and amputation in hemodialysis patients. My objective is to prove that all these have a correlation with heart attacks/ strokes and amputation in dialysis patients. After proving I want to find out which variable has the strongest link with ANY of these events. Now I have 88 dialysis patient in my hospital and have the variable data for a, b, c, d in excel form monthly for 5yrs. How can I use these variables a, b, c, d which will be 12 times in the year and 60 in 5yrs. How is time averaging done? Let me explain more. I will only take one adverse event heart attack for now. Increase in variable A and B increases risk of heart attack and increase in variable C and D decrease risk of heart attack. How should I design a restrospective study over 5yrs to maximize the chance of finding a significance between the above mentioned hypothesis? Variables a to d are not binary but rather continuous with limits ofcourse. Example variable C ranges from 1 to 4mg/dl in the blood. The optimum level is 4mg/dl, the lower the level the higher risk of heart attack. Now how should I analyse my data to maximize the chance of finding an association between heart attacks and low level of variable C. I will use the similar strategy for all variables. What test will I use and where can I get the software to do those tests? I am sorry I dont know the meaning of time average. There is a study that I heard about with a similar idea. They used a variable measured every month for a few years. I am not sure how they time averaged? Perhaps it is the arithmetic mean? This is the study that uses time averaging- Dialysis Modality and Correction of Uremic Metabolic Acidosis: Relationship with All-Cause and Cause-Specific Mortality Tania Vashistha, Kamyar Kalantar-Zadeh. It appeared in the journel CJASN. The link is here- http://cjasn.asnjournals.org/content/early/2012/11/21/CJN.05780612.abstract Is there a way to attach pdf then I can with this question. Thanks to all in advance


title: Graphical notation for chaining nodes over time

text: Is there any plate-notation-like notation used for chaining nodes in a graphical model? For example how to represent a Hidden Markov Models graphically and compactly without using repetitions: $Z_t\rightarrow Z_{t+1}, Z_t\rightarrow Y_t$ I do not like the common graphical representation of Markov Models with the $1, 2, ..., n$ and each step is a graphical model in itself.


title: What is needed for a line to be an appropriate model for the data?

text: My answer key responds with a large $R^2$ and a residual plot with no pattern. The residual plot with no pattern makes sense. However, is having a large $R^2$ necessary? A linear model can be appropriate, but weak, with a low $r$ value, thus it follows that a model can be appropriate, although weak with a low $R^2$ value. Is the answer key padding with extra information? (Note, I'm currently in AP Statistics, and I'm sorry if this is a really basic question! Thank you for taking your time, though!)


title: Different histograms for same data in R and Oracle Crystal Ball

text: I am analysing small sample size data:               >dput(dat.demand2050.unique)       c(79, 56, 69, 61, 53, 73, 72, 86, 75, 68, 74.2, 80, 65.6, 60, 54)        for which the histogram and density function look like this in R: http://db.tt/1uVSufat But when I used the Oracle Crystal Ball Simulation tool (it is an Excel plugin, but I can't post the link due to low reputation) to get the probability distribution, it gives a different result: http://db.tt/ugrpn140 Does anyone have experience with this tool? I am pretty sure the "probability" shown on the Y-axis is the same as the "Density" for the plot obtained in R, but I could be wrong. Another question is if the plotted quantities are indeed the same, what is the reason for the difference? Is it due to a difference in the algorithm?


title: Calculating mean and variance with logarithmic sample weights

text: I have run into a problem that must be pretty simple, but I keep getting snagged somewhere. I have an algorithm that returns a sample and the logarithm of the sample weight (which get themselves pretty large, residing in the range 1e2 - 1e4 for typical applications). Now, to take the sample weighted mean, one can use the formula: $\langle Q \rangle = \frac{ \sum_i w_i Q_i }{ \sum_i w_i }$ Taking the log of both sides: $\ln \langle Q \rangle = \ln \sum_i w_i Q_i - \ln \sum_i w_i$ and applying the identity to both RHS terms (see Wikipedia): $\ln (a+b) = \ln a + \ln \big[1+\exp(\ln b - \ln a)\big]$ in the general form for a sequence of numbers: $\ln \sum_i a_i = \ln a_0 + \ln \bigg[1 + \sum_i \exp( \ln a_i - \ln a_0 )\bigg]$ would allow one to calculate the mean using only the logarithm of the sample weights, and not the sample weight directly (and similarly for the variance). The problem here, however, is obviously numerical stability. After some research, I have come across two parallel suggestions: (1) For the case of two points, $a$ and $b$, it's easy to order the terms such that $b > a$ and apply methods for calculating $\ln(1+\epsilon)$, with $\epsilon$ between 0 and 1. I'm failing to see, however, how this simplification generalizes to the case of more than two points, where the value of the sum (in the general formulation) will probably be above 1. (2) Using an algorithm for calculating the weighted incremental variance (as here). I'm also failing to see here how to apply this algorithm to the case of logarithmic sample weights, as the mathematical formula on which to apply the logarithm (and identities) is not clear. My apologies for the long question, but after some time of going through the research on tricks for making this calculation more stable, I am still getting quite confused, as the tricks either don't seem to apply generally, or lead to contradictions. Any thoughts on this problem?


title: If a measure is valid (but not necesarily reliable), can it be consistently replicated?

text: I had this question on an exam, and I was positive that the answer was A. > The Beck Depression Inventory is a scale intended to measure depression > levels, with higher scores indicative of higher levels of depression. If > this were a valid measure of depression, we would expect that: > > A) the results of the inventory cannot be consistently replicated. > > B) a person's score on the inventory is not related to his or her level of > depression. > > C) people who get higher scores on the Beck Depression Inventory are more > depressed than people who get low scores. > > D) people who get lower scores on the Beck Depression Inventory are more > depressed than people who get high scores. I don't think my answer was wrong. If the test was valid, it is not necessarily reliable. Just because you are taking the same test, you are not going to get the same score every time. You are going to get a score that reflects your depression level at the time of taking the test. Which means that the results cannot be consistently replicated.


title: A problem on correlation and independence

text: Here is a little challenge/problem for you guys. Let $(X,Y)$ be jointly discrete RVs such that each have at most two mass points (not necessarily $1$ and $0$, i.e. not necessarily indicator variables). Suppose $X$ and $Y$ are uncorrelated. Are they independent?


title: Expectation of the Sum of K Numbers

text: Given $n$ numbers, where the value of each number is different, denoted as $v_1, v_2, ..., v_n$, and the probability of selecting each number is $p_1, p_2, ..., p_n$, respectively. Now if I select $K$ numbers based on the given probabilities, where $K \leq n$, what is the expectation of the sum of those $K$ numbers? Note that the selection is without replacement, so that the $K$ numbers cannot involve duplicate numbers. I understand that if the selection is with replacement, the expectation of the sum of the $K$ numbers equals $K \times E(V)$, where $$E(V) = v_1 \times p_1 + v_2 \times p_2 + ... + v_n \times p_n.$$ Furthermore, what about the expectation of the variance of those $K$ numbers? I am a CS PhD student who is working on a big data problem, and I don't have any statistics background. I expect that someone can give me a formula as the answer. However, if the answer is too complicated to be described by a formula or intensive computation has to be involved, an approximate answer is totally acceptable. You can assume $n$ here is quite large, and the probability can vary a lot. In practice, the values of those probabilities come from a query log, which records a series of aggregation queries. The point is that the frequency of each number involved in the queries can be quite skew, i.e., some are rarely queried, while some are queried very frequently. You can assume the probability distribution is normal distribution, zipf distribution or any other reasonable alternatives. The value distribution is only a contiguous subset of any possible distribution. In other words, if you have a histogram that represents a certain distribution, the all the numbers involved in this problem are the numbers all within a single bucket. In terms of the value of K, you may assume it is always less than the number of frequently queried elements.


title: Is there an error in the one-sided binomial test in R?

text: I am not sure, but I think R is doing not what it is supposed to do. Using `binom.test()` my understanding of the parameter `alternative="greater"` is the following hypothesis. $$H_0: p \le p_0 \quad \text{vs.} \quad H_1: p > p_0$$ (I mean it says greater and not greater or equal). P-Value then should be calculated as follows: $$p(c) = P_{H_0}(Z>c) = 1 - P_{H_0}(Z\le c)$$ Here is the equivalent R-Code               binom.test(8, 10, alternative="greater", p=0.5)$p.value     #[1] 0.0546875          1-pbinom(8,10,0.5) # 1-P(c<=Z)     #[1] 0.01074219          1-pbinom(8-1,10,0.5) # 1-P(c<Z)     #[1] 0.0546875      So where is my mistake? Or is R just a little imprecise? And if I am right and R is realy testing the "wrong" (="not exactly what one would expect") hypothesis: What is with the other tests? Does "greater" always mean >= ?


title: exponential approximation at adaboost

text: I am reading about AdaBoost algorithm, ![ppt](http://i.stack.imgur.com/uEMxy.png) I cannot understand (7). How can it be like this? If you want full document, you can get it at http://cseweb.ucsd.edu/classes/fa01/cse291/AdaBoost.pdf Thank you for advance.


title: Which prediction algorithm fitted for prediction in R

text: I was used R and mongodb for finding predictions of next date outcome for that I write R code as below               library('RMongo')     mg1 <- mongoDbConnect('demo','localhost',27017)     query <- dbGetQuery(mg1,'hosts',"{'hostId' : '101.10.202.10'}")      date <- query$runtimeMillis         memory <- query$memoryUtilization     plot(date,memory)     lm1 <- lm(memory~date)      when I plot this graph it looks like below ![enter image description here](http://i.stack.imgur.com/57OkC.jpg) Now I want to predict given date memory utilization so I add following code for finding prediction                 new <- data.frame(date=1377843220)       prediction <- predict(lm1, newdata=new, interval="predict")      but given prediction not match actual values I also tried interval="confidence" it also shows me wrong results. How should I find out prediction of given data is linear regression algorithm fitted here or any other algorithm helps? Data excel file here Download


title: Usage of Linear optimization model

text: Young Energy operates a power plant. The power plant is a coal-fired boiler that produces steam which in turn drives a generator. The company can buy different types of coal, and then mix them to meet the demands placed on it which is fired in the boiler. The table shows the characteristics of the different types of coal are: ![enter image description here](http://i.stack.imgur.com/Jg23l.png) The requirement to be burned in the pan is:   * BTU/lb: 11900,    * content of the ashes max 12,2% and    * max moisture 9,4%. > How should I implement a linear optimization model in this context?


title: How to measure user trust in crowd-sourced data?

text: I am developing an application that will feed itself from users, all data will come from users crowdsourcing it with GPS positions. Said data must be pretty accurate or else my application will not work, not at all. I've setup a trust index subsystem where each point of data has as a trust index, and said index is bumped up each time data with similar characteristics is uploaded, in other words, if point A is located at 4.33, -75.33 (Made up point) and another user uploads a 4.25, -75.30 (B) point, i can be sure to bump the trust index for A instead of creating another point altogether (Assuming the physical distance from A to B is 5 meters as a maximum). Any suggestions on how to enhance this algorithm? Any other algorithm i could use to cross-validate all the new data with previous data? Any paper i can read about this topic? As for more information, this trust index will be used to create different routes (as in sets of points) to navigate a map so the system can take in account only those points with certain minimum trust index or outright discard points that never got the minimum trust in a period of time from the database to save space.


title: Metric for nearest neighbor method

text: Is there a requirement that the measure used in Nearest Neighbor methods be a proper metric distance? What will happen if I use an arbitrary function (e.g., one that does not satisfy the triangle inequality)? How can you convert an arbitrary function to a valid metric distance?


title: Is it valid to keep insignificant variable?

text: I'm building a model that predicts rather rare event. The data set consists of about 240.000 examples from which only 1410 are positive. One of the variables from the data set, let's call it $X_{1}$ is a binary variable with value 1 in only about 1% examples. From the knowledge about the event I know that it is **impossible** to occur an event when $X_{1} = 1$. When I include $X_{1}$ in logistic regression model then it has highly negative coefficient, but it is also highly non significant. The coefficient is -12 and the standard error is 122. **Question** :   Is it valid to keep the variable $X_{1}$ despite the low significance?


title: I want to show $E(B(t)-B(s))^4=3(t-s)^2$

text: Let $B(t)$ and $B(s)$ are brownian-motion I want to show $$E(B(t)-B(s))^4=3(t-s)^2$$ thanks for help.


title: Monte Carlo estimation of probabilities

text: I would appreciate some advice on how to use Monte Carlo for estimating probabilities. Generally speaking the problem I have involves running an experiment and counting the frequency of output (which depends on random inputs), and then using the relative frequency to estimate the probability of an output. However, the starting conditions is also randomly generated. So I am repeatedly doing this:   1. randomly generate starting conditions according to set rules,   2. randomly generate input,   3. record output. My questions are:   1. Is this a reasonable method?    2. How can I estimate errors for the probability estimate? Any help appreciated; I am new to this type of thing!


title: Covariance of conditionally independent random variables

text: I want to find the covariance between two random variables $X$ and $Y$, which are independent given another random variable $M$. I thought the calculation would be $$cov[X,Y]=E_M[cov[X,Y|M]] =E_M[E_{XY}[(X|M-E_X[X|M])(Y|M-E_Y[Y|M])]].$$ But since $X$ and $Y$ are conditionally independent, $cov[X,Y|M]=0$, and thus, $E_M[cov[X,Y|M]]=0$. Obviously this is false. Would replacing $E_X[X|M]$ with $E_M[E_X[X|M]]$ (and doing the same for $Y$) give me the correct calculation? A simple yes or no answer will suffice, but I would be especially appreciative of a simple derivation of $cov[X,Y]$ in this situation.


title: What strategy can be used to construct a multilevel model with a modest amount of groups?

text: Here is my analytic situation:               Dependent: retirement age     Independent: job mobility (categorical, dummified)     Control variables: gender, age, education     Levels: about 15 different countries.       Ultimately I am not interested in comparing countries among each other, instead I would like to know whether career mobility has an effect on retirement age. I need to remove the country effects in order to know that. This means removing all variation due to local pension regulations, culture, economical parameters etc. One option would be to have 15 different models (for each country), another would be to include a country control variable in a simple linear regression model. My question is whether there is a statistically sound way to use a multilevel model in this case, using 15 countries as groups (2nd level entities). Thank you.


title: Removing macro-level time variance

text: The title might be a bit misleading. Unfortunately statistics is not my area of forte, so gentle guidance along the right path is much appreciated. That said, here's my problem: I'm working on analyzing precise time measurements from an experimental process. I have 16 measurement groups. Of these, one should (on average) differ slightly in magnitude from the others. The remaining 15 should have similar or identical properties. I'm trying to determine which of these measurement groups is the odd one out. My first pass at the problem involved a naive data collection routine rotating between sample groups until enough samples from each group were collected. My data is clearly not normal. Furthermore, sample values seem to fluctuate as a function of time (non-scientific visual inspection shows that larger values and smaller values tend to cluster across sample groups). The variable under consideration should not be fluctuating like this, so I believe this effect is attributable to some other element in an admittedly noisy environment. My first pass at analysis was to run a KS test comparing each sample group to the set consisting of all other collected data. This worked as expected: the sample which was known to vary showed a much lower p-value than the other candidates. The KS test isn't ideal here. Despite a theoretical belief that the data is continuous, my measurement resolution limits dictate that there must be a number of duplicate data points (3 effective places of resolution and n between 10k and 1m per group). Unfortunately, this sample set was taken from a somewhat idealized experimental condition. I need to improve my analysis to include some larger elements of noise. I have done the experiment in my real-world condition, and have been having trouble finding reliable ways to analyze my data. The problem seems to be related to the fact that while my effect remains constant over much larger (and thus longer-running) sampling regimes, the macro-level time changes seem to be interfering. As I move to an order of magnitude more samples in my real-world condition, rather than converging on one single value as different, the KS test starts to show that ALL the values differ (and the ordering is inconsistent). I've re-written my data collection to log both the timestamp for each sample, and several timing values taken over the course of my experimental procedure (originally I was just logging a single aggregate timing result by calculating a difference using these values). With this expanded data in hand, I've been looking for appropriate analysis techniques. Here's where I'm stumped. I'm not sure what I should do with this timestamp data. I've turned a multiple sample comparison problem into a time sequence. If I can use the time data to even out my values, I'm still (at least potentially) left with a complicated multi-variable analysis. PCA assumes normality so that's out. ICA isn't appropriate because I know that most of my measured variables are not, in fact, independent of each other. Once I've sorted that out, I still need to compare to the overall group and figure out which one is different. Because each sample-per-group has a discrete timestamp, I don't have a simple way to compare across my sampling groups. It seems like the Wilcoxon–Mann–Whitney (Wilcoxon rank-sum) test is likely to be of use here, but only after I've done work on the underlying data. Suggestions? Edit: for what it's worth, this is a practical problem. Approximate solutions and things that narrow the candidate field are perfectly reasonable. No points for excessively pretty P values here.


title: Does sampling error include measurement error?

text: Gross sampling error (MSE) appears to be a composite of two errors sampling and measurement error. How do we assess measurement error ? can we find out net sampling error ?


title: If mean is so sensitive, why use it in the first place?

text: It is a known fact that median is resistant to outliers. If that is the case, when and why would we use the mean in the first place? One thing I can think of perhaps is to understand the presence of outliers i.e. if the median is far from the mean, then the distribution is skewed and perhaps the data needs to be examined to decide what is to be done with the outliers. Are there any other uses?


title: Calculating confidence intervals via bootstrap on dependent observations

text: The bootstrap, in its standard form, can be used to calculate confidence intervals of estimated statistics provided that observations are iid. I. Visser _et al._ in "Confidence Intervals for Hidden Markov Model Parameters," used a parametric bootstrap to calculate CIs for HMM parameters. However, when we fit an HMM on an observation sequence, we have already assumed the observations are dependent (in contrast to mixture models). I have two questions:   1. What does the iid assumption do with the bootstrap?   2. Can we ignore the iid requirement in a parametric bootstrap? * * * Visser _et al._ method is briefly as follows:   1. Assume we have an observation sequence $Y=o_1,o_2,...,o_n$ resulted from sampling an HMM with real but unknown set of parameters $\theta=\theta_1,\theta_2,...,\theta_l$.   2. The parameters can be estimated using the EM algorithm: $\hat{\theta}=\hat{\theta}_1,\hat{\theta}_2,...,\hat{\theta}_l$   3. Use the estimated HMM to generate a bootstrap sample of size $n$: $Y^*=o^*_1,o^*_2,...,o^*_n$   4. Estimate parameters of the HMM according to the bootstrap sample: $\hat{\theta}^*=\hat{\theta}^*_1,\hat{\theta}^*_2,...,\hat{\theta}^*_l$   5. Repeat the steps 3 and 4 for $B$ times (e.g. $B$=1000) resulting in $B$ bootstrap estimations: $\hat{\theta}^*(1),\hat{\theta}^*(2),...,\hat{\theta}^*(B)$   6. Calculate the CI of each estimated parameter $\hat{\theta}_i$ using the distribution of $\hat{\theta}^*_i$ in bootstrap estimations. Notes (my findings):   1. The percentiles method should be used to calculate CIs in order to have correct coverage (normality is a bad assumption).   2. Bias of the bootstrap distribution should be corrected. Meaning that the distribution mean of $\hat{\theta}^*_i$ should be shifted to $\hat{\theta}_i$


title: Ways of finding associations in time series

text: I am trying to develop a logic to identify association between different time series for association mining. I have a lot of series and need to find whether or not the association exists. I figured two ways to do this:   1. Constructive: try to find relation in given two series.   2. Destructive: try to prove that relation does not exist. Are there any existing mathematical parameters to identify such relations? if not, any suggestions?


title: Confusion related to variance and mse

text: I was reading this wikipedia article and it states that MSE of a predictor is equivalent to variance of the error. To test it I did something like this               a = randn(10,1);     b=randn(10,1);          sqrt(mse(a,b))     std(a-b)      But the outputs are different, I get 1.3430 and 1.4292. I am not sure what's going on here. Any suggestions guys?


title: Times series and starting point incidence

text: I'd like to **modelize** the height of a plant which depends on   * Its aging (in months).   * The month it has been planted. _For instance a plant of 3 months is higher if it has been planted in march than if it has been planted in November, but this difference tends to be null over the time._ _I.e a plant of 2 years which has been planted in march has the same height than a plant of 2 years which has been planted in November_ Which model fits whith these two requirements ? ps : I'd like to modelize it with R, but I think the major time series models are available in R.


title: Simulating a bimodal biased IV estimator

text: How can I simulate a bimodal biased IV estimator? The common unimodal heavy-tailed biased estimator would be something like this:               IV_data <- function(beta=10, lambda=-20, theta=0.5, n=50, mu=100, sd=2){       z = rnorm(n,mu, sd)       u = rnorm(n, 0, sd)       x = theta*z + u        y = beta*x + lambda*u + rnorm(n, mu, sd)       data.frame(y, x, z, u)     }          library(sem)     set.seed(20)     beta_IV <- replicate(1000, coef(tsls(y~x, ~z, data=IV_data()))[2])     d<-density(beta_IV)     plot(d)      I guess if I tweak the parameters a bit I could get a bimodal version. Or do I have to have a complete different model structure to get it?


title: Does the actual coverage of a 95% CI on $R^2$ get closer to nominal coverage with larger sample size?

text: If the answer is "it depends", what does it depend on? Does convergence depend on the ratio of predictor variables to sample size, or the size of $R^2$, or something else? I am mainly interested in CIs on Unadjusted $R^2$, but would also be interested to hear the answer to this question in relation to CIs on Adjusted $R^2$. Primarily I'm interested in finding out the conditions under which convergence happens or doesn't happen, but I would also appreciated an explanation of why. EDIT: I am at present using the calculator made here: http://www.danielsoper.com/statcalc3/calc.aspx?id=28 ![CI formula](http://i.imgur.com/A326ZYV.png) ![Standard error](http://i.stack.imgur.com/c3rFH.png) I'm interested to know whether this formula gets close to nominal coverage with larger sample size (or stays the same, or gets further away from nominal coverage). I'm also interested to know if this holds true for all $R^2$ CIs or just ones calculated using this formula.


title: Interpreting coefficients of Regression Model (Mincer Model)

text: ![Mincer Model](http://i.stack.imgur.com/UBJV6.png) Hi all, I am an undergraduate student who is currently doing an assignment. I am now facing a few problems which are:- 1) Age is usually a positive return to wage, but in my regression output, that's not the case. How should I interpret the signs and its coefficient then? And regarding the age^2, which is to test the concavity will also be odd. Any thought? 2) Higher is believe to be the proxy for innate ability and instrument variable (mother educated is used to solve endogeneity), is it correct the way I do it? 3) The t-value for certain variables are relatively low and insignificant, is interpreting those variables are still relevant and appropriate? Thank you in advance for your kind assistance! Stata/Statistics Guru, please help me out!


title: Best method of statistical analysis for large dataset of different points

text: I have generated a set of ground displacements at certain points using two different methods. I am now trying to find a good method of statistical analysis to compare how similar the results of the two methods are. I only have two measurements for each point but I have over 1 million different points. Is there a method of analysis that I can use to compare the similarity of the displacement generated by both of the methods at each of the different points which will then allow me to calculate the similarity of the overall data set?


title: Growing a classification tree manually

text: I want to create some CARTs in R. I used Clementine a while ago for this purpose and remember I could go into "manual mode" and grow the trees by hand. I could manually choose which variable to use for what split and so on. Can this be done in R?


title: Non-parametric for two-way ANOVA (3x3)

text: My dependent variable is continuous, non-normal (skewed left according to Shapiro-Wilk test). I have two independent variables (treatment group by colour, food type). There are 3 levels within each independent variable. The number of observations for each independent variables are not equal. I have looked up non-parametric tests such as Friedman's Test and Scheirer- Ray-Hare Test, neither of which seem suitable (due to unequal number of observations). Are there alternative tests that anyone could suggest? I am using SAS.


title: Multilevel modeling of response time data

text: I'm trying to figure out how to set up and analyze the following experiment. It's a basic reaction time-type experiment with 4 independent variables (2 levels each) and 1 dependent variable (RT). The 4 independent variables are:   1. Hand used for response (Right vs Left)    2. Side of visual field where target appeared (left vs right)   3. Exposure Time (400ms vs 800ms)   4. Priming of target (Primed vs Not-primed) Basically I want to see how these things interact in order to test a specific hypothesis about their interaction. The problem with this experiment is the fact that I will be generating about 30 observations per participant for each unique combination of those variables. From what I can gather repeated-measures ANOVA would not be the way to go as this would require reducing those multiple observations to a mean, losing valuable data. Some posts here (and elsewhere) suggest multilevel modeling for analyzing this type of experiment. What would be the most efficient way of setting up this experiment and analyzing the data? Should I look at grouping the variables a certain way or can I just chuck the data into long format, use the MIXED function in SPSS and call it a day? From the pilot testing it appears as if the distributions will be positively skewed as with most RT data of this nature. Should i be performing some log-like transformation before doing the analyses? p.s Is this thing even called a 2x2x2x2 factorial design?


title: Questionnaire analysis and grouping responses

text: I want to know if the creative process engagement (CPE) is related to creativity. I don't know how to analyze the questionnaire. I will be collecting data for CPE with Zhang and Bartal's (2010) 11-item scale and creativity with Shin and Zhoe's (2007) four-item scale. Responses will be recorded on 5-point Likert scales. I know how to input data into SPSS. What I don't know is how to analyze it. Should I enter data for each variable in separate files and analyze it (e.g. creativity q1-q4 in one file and creative process engagement q1-q11 in another file)? How will these questions tell me if respondent thinks there is creativity or not and if respondent engages in creative process or not? Will I have a binary value after analysis (0,1) for each respondent? (i.e. if n=10 will i have 10 values for variable say for creativity) Second, data on creativity is provided by supervisors and data on creative process engagement is provided by employees. Now this means I'll end up with 3-8 (depending on firm) responses for cpe for each response of creativity (i.e. if y=1 then x=3. supervisor having 3 workers). What should I do if I want to use regression? Can the answers of workers (x=3) grouped (x=1) so I have 1 response?


title: Term selection in a model

text: In selecting terms to include in a model, say a linear one, should we always test the significance of the main effects first, keep only the significant ones, then consider the possible interactions between them? For example, if I have five covariates $x_1$ - $x_5$, and I choose the significant main effects in some way. Assume the first three are significant, $\hat{Y} = x_1 + x_2 + x_3$. Now I consider the interaction between $x_1$ - $x_3$. Finally I may end up with something like $\hat{Y} = x_1 + x_2 + x_3 + x_1 \cdot x_3$. Is it justifiable to do this? A problem may be the significant interaction between $x_4$ and $x_2$ ?


title: What to do if transformed data still has unstable variance

text: After using Bartlett test statistics for equal variance to test my data, i concluded that there was significant difference in the variance at alpha = 0.05. I transformed the data using Box-Cox family of power transform and tested for equal variance and it still not equal. what next to do? I used all the different transformation, but still not stable. My estimating parameter is 0.86


title: Determining the quality of low-pass filter estimate

text: I'm taking a low-pass filter, specifically an exponential moving average, of samples coming in every second. Those samples are of electric field strength, and suffer sizable noise. How can I determine the quality of the estimate provided by the filter?


title: How to test whether a series data follow Ornstein-Uhlenbeck process (OU process)?

text: I have some measures which seems to have some mean-reverting properties and I'm wondering whether they can be modeled as Ornstein-Uhlenbeck process (OU process). And actually I quite expect it because if so, the following work will be much easier because it seems that there has been many useful results of OU process. But I don't know how to validate it. I used to think that I could calculate the KL distance but the probability density function seems to not have a closed form, so it can't be derived. So I'm wondering how to test whether a series data follow an OU process?


title: The impact of non-normality (contaminated normal) on type 1 error and power when sigma is known

text: When sigma is unknown, contaminated normal distributions tend to have (1) the issue of the actual Type 1 error rate being lower than the nominal level and (2) the issue of reduced power. This pattern has been confirmed by data simulations. The question then is: does knowing the population variance help with (2) keeping the type 1 error at the nominal level and (2) having good power? Thank you in advance for your help! * * * To just add some info: By **power** , I simply mean correctly rejecting the null hypothesis of mu=mu0 when this h0 is actually false. I thought that contaminated normal distributions tend to cause wider CIs, which can make correct rejections more difficult. I did some data simulations to check contaminated normal distributions. Say I repeatedly sample from a contaminated normal distribution with mu=0 and sigma=3 for 1000 times. Then I create 90% confidence intervals based on these samples and check how many of these 1000 CIs contain the population mean. As it turned out (I tried this several times), the percentage of CIs containing the population mean is approximately 90% (that is when sigma is known). So in this case, it seems that the the actual type 1 error rate is comparable to the nominal level of 0.1. I also tried this assuming that sigma is unknown (so using t-statistics, rather than z), I consistently got more than 90% CIs that contained the population mean, so in this case, the actual type 1 error rate seems to be **lower** than the nominal level of .1. But I thought the type 1 error rate and power would suffer even when sigma is known (e.g., the actual type 1 error can be lower the nominal level). But my simulations (assuming that I did it correctly) at least did not show this. So I want to check with you guys to see what's going on. Thanks :)


title: multiple change point analysis simultaneously for mean AND variance WITHOUT distribution assumption in R

text: Is there a multiple change point analysis for both the mean and the variance simultaneously, implemented in R without a distribution assumption? I know the changepoint packet http://cran.r-project.org/web/packages/changepoint/index.html but as far as I see, it only allows change point analysis (mean & variance simultaneously) with distribution assumptions (normal, exponential, gamma).


title: How to calculate the probability that an algorithm classifies seven wines out of ten correctly when the true error is 0.23?

text: I am considering the following problem. > Calculate the exact probability that an algorithm classifies seven wines out > of ten correctly when the true error is 0.23. Should I solve this with Bayes' rule?


title: Method to find 'principal components' of EMG recordings of reflexes

text: I have a series of electromyographic (EMG) reflex recordings, of which I would like to find out the principal components. My thoughts are that there are multiple processes that added together result in the EMG recording. Here is a series of reflexes of a single recording ![series of reflexes](http://i.imgur.com/V70bK.png) and here a single reflex ![single reflex](http://i.imgur.com/COKlh.png) As you can see all reflex responses are stimulus time-locked and have clear a beginning and end (Note the vertical line is the detected latency in does not belong to the recording). I already tried to apply a principal components analysis (PCA) and got good results but there are negative factor loadings, which I can not easily explain physiologically. A non-negative matrix factorization (NNMF) does not allow negative input values so I offset the the signal to be strictly positive. I think both methods deliver, on first sight, good results, but how interpretable are they in the sense of multiple processes adding up to the compound EMG signal. My questions now are: Which method should I use? Are there better methods than PCA and NNMF for this task and can I just offset the values for NNMF. Thanks in advance and please let me know if I should explain something in more detail. -- Edit: Here is a sample result of the NNMF with two factors ![sample result](http://i.imgur.com/WYzT5.png) The upper figure shows the two components and the lower image the factor loadings for each response in the series. What can be seen is that at the start of the stimulation the green factor is more present and is 'replaced' between the 10th and 15th response with the blue one. Edit2: NNMF seems to only return correct results for rectified responses. Is there a method to restrict PCA for non-negative loadings?


title: Comparing groups of different sizes in social network analysis

text: I've got social network data for three groups of chimpanzees and I'm looking to do some basic analysis across and between groups. Does anyone by any chance know how to compare social network data between groups (for example, scores of strength from enclosure 1 vs enclosure 2), and/or 2) combine data from different groups (for example, adults from all three groups vs youngsters from all three groups)? The three chimpanzee groups differ in size (11, 26 and 40), adult/youngster ratio etc. I'm just exploring the data, but I just wondered if/how data can be compared independent of group, because of course social network data output from group 1 cannot be directly/simply compared to group 2, for example. Is there a way to standardize such data. At the moment i'm using SOCPROG for my analysis, but am happy to investigate other approaches if necessary.


title: Multilevel covariance structure and more

text: I have some doubts about the covariance structure in a multilevel model fitted in R (using the `nlme` package). I'm not an expert (just starting to learn statistics...), so I apologize if some of my questions seem evident. I've checked the previous posts and haven't found an answer. I have data from an experiment in which we have registered physiological data from 30 subjects in 2 conditions (with 30 trials in each condition). These 30 trials are close in time, and we expected a higher correlation between closer trials, with the correlation decreasing as trials are further apart from each other. We are interested in the effect of condition, not in time effects. I think that the right way to analyze these data is to fit a multilevel model, in which TRIAL is a level 1 variable, SUBJECT is a level 2 variable, CONDITION is a fixed factor, and the DV is the physiological response (FR). The R command I'm using is:               lme(fixed= FR ~ CONDITION, data=mydata, random= ~ TIME | SUBJECT)      My (many) questions are both theoretical and practical:   * Which covariance structure `lme` does use by default? Is it a problem not to use the most appropriate covariance structure?    * I've read that an autoregressive covariance structure (AR1) refers to a constant variance at each time point and a weaker correlation as time points get further apart. My data only meets the second criterion. How can I know which covariance structure is right for my data? How important is it to the validity of the results?   * I'm only interested in the CONDITION effect, which is significant when I'm not including TIME in the model. I'm not interesting in TIME effects, nor I want to use the model to make a prediction, but only to check the significance of the CONDITION effects. Is it correct if I drop out TIME and fit a model only with CONDITION? Thank you for your help!


title: Compare the variances of several groups

text: I have five groups of data and I want to test if their variance is significantly different or not. The data are normally distributed and I do not have any repeated measurements in any groups. What sort of analysis should I use? Tukey's test compares the means, but how can I compare their variances? For example I have 3 groups A, B and C. In groups A, B and C I have 30, 50 and 20 observations (which are stored as double-precision numeric variables in R).


title: How to derive the probability mass function for Zero-inflated Poisson Distribution?

text: I am a maths student and I asked to derive the zero-inflated poisson distributon. I looked it up in the interent and I found the p.m.f. however I have no-idea why it looks the way it looks. Thank you so much for answers.


title: Interpreting main effect coefficient in different models

text: My interest lies in finding the "right" correlation between a continuous IV ($x$) and a continuous DV ($y$). At first I ran a simple linear regression: $$ y=a+b_1 x $$ However, lots of other factors influence $y$ besides $x$. One of them is a categorical variable $c$, with $n+1$ categories. Using dummy coding, I ran the regression:   $$ y = a + b_1x + b_2c_1 + b_3c_2 + \ldots + b_{n+1}c_n $$ $c$ and $x$ are not orthogonal, so I then added interaction effects:   $$ y = a + b_1x + b_2c_1 + b_3c_2 + \ldots + b_{n+1}c_n + b_{n+2}xc_1 + \ldots b_{2n+1}xc_n $$ The interaction effect is significant (although the coefficients for some categories*x are not), as is the main effect, in all models, and (most) of the category coefficients. What I noticed, however, is that the coefficient $b_1$ is somewhat different in each of the three models. What is the correct way to interpret the correlation of $x$ and $y$?


title: Using least squares to estimate mean

text: From _"Generalized Additive Models: an introduction with R" by Simon N. Wood, page 55, Exercise 1.9, Question 1_ : > 4 car journeys in London of length 1, 3, 4 and 5 kilometres took 0.1, 0.4, > 0.5 and 0.6 hours respectively. Find a least squares estimate of the mean > journey speed. _Hint: use a simple one parameter linear model, with journey time as the response (y) variable and journey length as the predictor (x) variable, think about the meaning of the parameter._ The answer given at the back of the book is: $\frac{51}{6.3}\approx8.1$ kilometres per hour. **Question: is that answer correct? It seems to me that it is stated in squared kilometres per hour. Why not $\frac{\sqrt{51}}{6.3}\approx1.13356$ kilometres per hour?** * * * My working: This is a no intercept linear model because at time zero no distance will have been covered. $y_i=\beta x_i+\epsilon_i$ Using ordinary least squares regression I can obtain an estimate of $\beta$: $\hat{\beta}=\frac{\sum_{i=1}^4{x_iy_i}}{\sum_{i=1}^4x_i^2}=\frac{(1)(0.1)+(3)(0.4)+(4)(0.5)+(5)(0.6)}{1^2+3^2+4^2+5^2}=\frac{6.3}{51}=0.1235294$ So the model is: $y_i=0.1235x_i$ Which means a one kilometre change in $x_i$ results in a 0.1235 time change in $y_i$. This is where I start to get confused. The total distance covered by the four cars is $\sum_{i=1}^4{x_i}=13$ kilometres. The total squared distances covered by the four cars is $\sum_{i=1}^4{x_i^2}=51$ kilometres$^2$. The total time taken by these four cars is $\sum_{i=1}^4{y_i}=1.6$ hours. As we are told to use least squares, we should then surly do: $\frac{\sqrt{51}}{1.6}\approx 4.46$ kilometres per hour I'm guessing I'm probably wrong would appreciate an explanation :) EDIT: Made correction based on @Whuber's comment


title: Likert scale inter rater agreement

text: What is the best inter rater agreement test for Likert scale type questions. As far I saw that cronbach alpha is for internal consistency, or it shows how good items are related to describe the main question. I want to measure inter- rater agreement. Is Cronbach Alpha the right metric for doing this?


title: Is it abnormal that out-of-sample fit is better than in-sample?

text: I'm using Eureqa, as machine learning tool to fit a formula to my data. I found out that the formula fits my test data better than my training data! Is this abnormal?


title: Constructing most powerful critical region basics

text: I've been trying to teach myself some of the fundamentals of statistics and would like to see if I'm in the right direction. **Here's the problem:** Construct the most powerful critical region of size $\alpha$ to test the null hypothesis $H_0: \theta=\theta_0$, where $\theta$ is the parameter of a binomial distribution $b(n,\theta)$ with a given value of $n$, against the alternative hypothesis $H_1: \theta=\theta_1<\theta_0$. Use $\theta_0=0.75$, $\theta_1=0.5$, $n=25$, and $\alpha=0.05$. **My work:** Since we have simple hypotheses, we can apply the Neyman-Pearson lemma. So, if $f(x;\theta)$ is the pmf we're working with, we have $$\frac{f(x;\theta_1)}{f(x;\theta_0)}= \frac{ {25\choose x}(1/2)^x (1/2)^{25-x}}{{25\choose x}(3/4)^x (1/4)^{25-x}}= 2^{25}(1/3)^x.$$ Now, this ratio is decreasing in $x$, which means we reject $H_0$ for small enough $x$. So, we want to find $c$ such that $$\mathbb{P}(X\leq c \hspace{1mm}; \theta=0.75)=0.05$$ and, by computation, I found $$\mathbb{P}(X\leq 14 \hspace{1mm}; \theta=0.75) \approx 0.03$$ and $$\mathbb{P}(X\leq 15 \hspace{1mm}; \theta=0.75) \approx 0.07.$$ Hence, by Neyman-Pearson, a test that rejects when $X\leq 14$ is a UMP level $\alpha(=0.05)$ test. So, the most powerful critical region of size $0.05$ is $0\leq X \leq 14$. I would very much appreciate any feedback/corrections/etc. Thanks!


title: Indicators of disease

text: I have this dataset with size (rows) and shape (columns) and in each cell #disease/#in cell . How would I answer the following questions? Is size an indicator of disease?   Is shape an indicator of disease?   Is there a size*shape effect?               Size          shapea    shapeb     shapec    %a    %b    %c     <4 mm       5 / 5771    0 / 26    1 / 522  0.1%  0.0%  0.2%     4-6 mm      9 /  739    1 / 11    2 / 149  1.2%  9.1%  1.3%     6-8 mm      3 /  496    0 /  7    1 / 126  0.6%  0.0%  0.8%     >=8 mm     33 /  276   18 / 26    8 / 146 12.0% 69.2%  5.5%          Total      50 / 7282   19 / 70   12 / 943  0.7% 27.1%  1.3%      Thanks, Luther88 **Edits** @suncoolsu in detail these are nodule data. size is the size of the nodule (the cutoffs are what have been published and accepted as standards). shape is the shape of the nodule. the numerator is the number of malignant nodules and the denominator is the number of nodules of given size and shape combinations. @onestop i do not have the raw data to do logistic regression. just this table. any further suggestions?


title: Normalizing data between 0 and 1

text: I am a novice when it comes to stats, so I apologize beforehand for the simplicity of my question. I trying to figure out the best way to normalize (this may be the wrong term) my data in so that the maximum value is 1 and the minimum value is 0. What are my options and what methods would you suggest as the best? My data will be non-negative and has no bound (though it may be bound in some situations, how would this make it different?). I am measuring content across many different dimensions and I want to be able to make comparisons in terms of how relevant a given piece of content is. Additionally, I want to display values across these dimensions that is explicable and easily understood. Thank you in advance for helping a newbie :)


title: Fitting 4-moment distribution with mixture gaussian

text: I know that Mclust does the fit on its own but I am trying to implement an optimization with the aim to generate a mixture of 2 gaussians with the combine moments as closed as possible to the moment of my returns' distribution. The objective is to `Min Abs((Mean Ret - MeanFit)/Mean Fit) + Abs((Std Ret -Stdev Fit)/Stdev) + Abs((Sk Ret-Sk fit)/Sk Fit) + Abs((Kurt Ret- Kurt Fit))` Taking into account that I fix the weight between the two gaussians at (0.2;0.8) I implement the below code in R:               distance <-function(parameter,x) {       u=mean(x)       s=sd(x)       sk=skewness(x)       kurt=kurtosis(x)       d1=dnorm(x,parameter[1],parameter[2])       d2=dnorm(x,parameter[3],parameter[4])       dfit=0.2d1+0.8d2       ufit=mean(dfit)       sdfit=sd(dfit)       skfit=skewness(fit)       kurtfit=kurtosis(fit)                 abs((u-ufit)/ufit)+abs(s-sdfit)/sdfit)+abs((sk-skfit)/skfit)+abs((kurt-kurtfit)/kurtfit))      }      Parameter<-c(0,0.01,0,0.01)  # starting point of the optimization      opp<-optim(parameter,distance,x=conv)         1. could anybody tell me whether it is the right approach ?    2. should I add some constraint like `ufit=0.2*mean(d1)+0.8*mean(d2)...` thank you very much in advance for your time and help. Sam


title: Sampled median's accuracy

text: I'm working on a problem where I need to calculate the median for a very large data set (for instance, 100M values) that has a log-normal distribution. Because of the data set's size, we were thinking about taking a sample (say, a random subset of 2000 values), and calculating its median. While this is much nicer from a computation perspective, I am very worried it will be inaccurate. What method could I use to determine how accurate this sampled median is?


title: Confidence interval when data is in range

text: I am trying to calculate the confidence interval when data is given in range. Example: Rather than telling their exact age in years, 9 sample students from a class tell their age as: > (23 +/-1) ; (22 +/-1) ; (20 +/-2 ) ; (25 +/-2) ; (24 +/-2) ; (25 +/-3) ; (22 > +/-3) ; (24 +/-4) ; (25 +/- 4 ) How do we calculate XY% confidence interval in such scenario.


title: Pseudo-$R^2$ in spatial regressions and Kelejian-Robinson test in R

text: I am currently trying to find out if there is a way to calculate a pseudo $R^2$ value from the output of maximum likelihood estimation of the spatial error model, and maximum likelihood estimation of the spatial lag model in R. I wonder also if there is a way to calculate the Kelejian-Robinson test to verify the spatial error dependence in R.


title: Validating the CART model in R

text: I have built the CART model, however I want to understand how we predict/validate the results with Validation data.               d=sort(sample(nrow(red_data),nrow(red_data)*.7))          #Selecting training sample     train_data=red_data[d,]     test_data=red_data[-d,]     nrow(train_data) #20133#     nrow(test_data)#8629#          rpart(formula = red_status ~ VALUE_CLASS + GAP_1 + cs_rr1_2,          data = train_data, method = "class", control = rpart.control(minsplit = 6))          print(cfit)     summary(cfit)          prp(cfit,type=2, extra=106, nn=TRUE,fallen.leaves=TRUE)      This is what I have done. Could some one help me on this


title: Importance Sampling to evaluate integral in R

text: I have asked the question here also. However, there might be something wrong with my theoretical understanding hence I'm asking here as it is more relevant. Kindly do not diss without looking first. I have referred to these links - here and here. The latter link is a question of mine, but the answers, as they turn out, aren't quite right. Hence I am creating asking a question with new attempts made to answer it.   * I have the following integral                          (x^-0.5)                                   ;in x=[0.01,1]      To solve this using Importance Sampling MC integration, one needs to select an **importance pdf** that is approximately the same as the **function plot** My R code to solve the same is this :               #function 1 - importance sampling     w <- function(x) dunif(x,0.01,1)/dbeta(x,0.7,1)     f <- function(x) x^(-0.5)     X <- rbeta(1000,0.7,1)     Y <- w(X)*f(X)     c(mean(Y),var(Y))      **True integral value** - 1.8 **Using the Importance Sampling code above** - 1.82 (where my importance PDF is **_Beta(0.7,1)_** which is quite alright so I'm assuming the code is correct.   * Now I have this integral  `[1+sinh(2x)ln(x)]^-1 ;in x=[0.3,8]` for which my code is :               #function 2      w <- function(x) dunif(x,0.01,1)/dnorm(x,0.5,0.25)     f <- function(x) (1+sinh(2*x)*log(x))^(-1)     X <- rnorm(1000,0.5,0.25)     Y <- w(X)*f(X)     Y <- Y[!is.na(Y)]     c(mean(Y),var(Y))      **True Integral Value** ~0.7014 **Value from executing above code** ~3.25 (where my importance PDF is **_normal(0.5,sd=0.25)_** * * * What am I doing wrong? 1) Take function to be evaluated as f(x). 2) Generate samples from Importance PDF g(x) that is truncated between the intervals. 3) Get the mean(f(x)/g(x)) which is the integral.


title: Validating a Data Set.. Which Test should I use

text: I have a data set(size around 1000) which was collected by a third party. I want to verify the correctness of this survey data using random sampling (by telephoning). Is this a correct approach to the verification? Also, which statistical test would be the best choice?


title: ML vs WLSMV: which is better for categorical data and why?

text: I was wondering which is a better estimator to use for categorical data: ML or WLSMV. I saw on a discussion on the Mplus website that they recommend WLSMV for categorical data but didn't explain why. Does anyone know specifically why ML doesn't work as well? Preferably, I am looking for a reference that compares these two estimation approaches, but have not been able to locate one after hours of searching. Thank you for sharing your knowledge and experience!


title: Modified K-means with unequal cluster variances

text: I wonder how I can modify the K-means algorithm so that the cluster volumes are not equal to each other. The K-means objective is to minimize within cluster sum of squares $\sum_{i=1}^{p} {\parallel \mathit{X}_i-\mathit{L}_{\mathit{Z}_i} \parallel}_2^2$, and this objective assumes that all cluster variances are the same. If we assume that the clusters are Gaussian with mean $\mathit{L}_{\mathit{Z}_i}$ and variance $\sigma_{\mathit{Z}_i}^2$ where $\mathit{Z}_i$ stands for the cluster assignment of data point $i$, then the objective for the cluster assignments become $\sum_{i=1}^{p} \frac {{\parallel \mathit{X}_i-\mathit{L}_{\mathit{Z}_i} \parallel}_2^2} {\sigma_{\mathit{Z}_i}^2}$. So, I tried modifying K-means such that $\mathit{Z}_i$ update is performed using this new update rule, and $\sigma_{\mathit{Z}_i}^2$ are also updated in each iteration. However, when I use this new modified K-means, almost all data points are assigned to the same cluster, which is weird. What might be the problem about that approach? I know EM can be used for this unequal-volume GMM purpose, but I want a simpler approach like K-means, and I am really curious about why what I tried is not feasible. Thanks!


title: Bayesian network to factor graph

text: I have a bayesian network with conditional probabilities as given by the diagram and I have converted it to factor graph. I just sort of read about factor graphs. Can someone be kind enough to let me know if I have done the factor graph representation correctly? I really appreciate your help. ![Graphical Model](http://i.stack.imgur.com/DwNCe.png) ![Factor Graph](http://i.stack.imgur.com/wrigJ.png) EDIT: After the comments, I thought about this again and came up with the following to express the joint distribution $$ P(w, \lambda, \phi | y) = P(y|w,\phi) \times P(w|\lambda) \times P(\lambda) \times P(\phi) $$ I hope it is correct. Note that $y$ is observed. Any suggestions would be greatly appreciated. ![enter image description here](http://i.stack.imgur.com/jKELT.png)


title: Comparing multiple gls' objects using anova.gls

text: I got a warning message when I was trying to do anova for two nlme::gls objects. Here is an example:               require(nlme)     set.seed(123)     y<-rnorm(100,10,2)     x1<-rnorm(100)     x2<-sample(1:5,100,T)     x3<-rt(100,20)     x4<-rbinom(100,1,0.3)     fit1<-gls(y~x1+x2+x3+x4,correlation=corAR1(form=~1|x4))     fit2<-gls(y~x3+x4,correlation=corAR1(form=~1|x4))     anova(fit1,fit2)          Model df      AIC      BIC    logLik   Test  L.Ratio p-value     fit1     1  7 421.5587 439.4359 -203.7794                             fit2     2  5 413.9855 426.8590 -201.9928 1 vs 2 3.573242  0.1675     Warning message:     In nlme::anova.lme(object = fit1, fit2) :       fitted objects with different fixed effects. REML comparisons are not meaningful.      Dose anybody know what the message means? Should I do the following rather than above?               > anova(fit1, L=c(-1, -1, 1, 1))     Denom. DF: 95       F-test for linear combination(s)     (Intercept)          x1          x2          x3               -1          -1           1           1        numDF  F-value p-value     1     1 281.9844  <.0001      Why the numDF is one here? It seems 2?


title: Wilson score for a range of possible values?

text: (Note: I did look at similarly worded questions and none seem to answer this one.) Given this scenario: People can rate 0 or more items on a scale of 0-X. I am trying to find the Wilson score interval (WSI) for each item with at least 1 rating. Does the WSI apply to such situations? If so, can Anyone point Me in the direction of an example of how this is applied in this scenario?


title: The difference between linear SVM and other linear classifiers?

text: The linear SVM in textbook takes form of maximizing $L_D = \sum_i{a_i} - \frac{1}{2}\sum_{i,j}{a_ia_jy_iy_jx_i^Tx_j}$ over $a_i$ where $a_i \geq 0$ and $\sum_i{a_iy_i} = 0$ Since $w = \sum_i{a_iy_ix_i}$, the classifier will take the form $\text{Sgn}(wx - b)$. Thus, it seems to solve linear SVM, I need to figure out $a_i$ with some gradient based methods. However, recently, I came across a paper which states that they try to minimize the following form: $L_P = \frac{1}{2}||w||^2+C\sum_i{\text{max}(0, 1-y_if_w(x_i))}$ and they claim $C$ is a constant. It seems to me that this form is quite different from primal form of $L_P$ in linear SVM because of the missing $a_i$. As far as the paper goes, it seems to me that they optimized on $w$ directly. I am puzzled here as if I missed something. Can you optimize $w$ directly on linear SVM? Why is that possible?


title: Brewer's method for sampling with unequal probabilities with n>2

text: I'm trying to use Brewer's method to sample 12 units out of a population of 73. I read on Brewer and Hanif's "Sampling with unequal probabilities" that the probability must be proportional to $$\frac{P_i (1 - P_i)}{(1-rP_i)}$$. I know that for $n = 2$ the first unit is selected with probability $$\frac{P_i (1 - P_i)}{D (1-2P_i)}$$ where $$D = \sum_{i=1}^N\frac{P_i(1-P_i)}{(1-2P_i)}$$ and the second unit is selected with probability $$\frac{P_j}{1-P_i}$$. I would like to know if this selection probabilities are the same for $n > 2$ or if I have to change the probabilities for the first unit according to the sample size.


title: How do I compare sample means in this experimental-control group study?

text: I am hoping to compare the sample mean of an experimental group to the sample mean of three control groups. Here is some context for the study: The experimental group consists of ~60 persons with a particular criminal record who were afforded various rehabilitative services. For each person, I have a count of their number of rearrests after they began their "treatment." The control groups are similarly sized (50-70 people) and share similar characteristics to the experimental group. For each person in the control group, I have a count of their number of arrests after the control group began their treatment (i.e. control group is monitored over the same time period). I now want to test if there is a statistically significant difference between the number of rearrests for the experimental group and the control groups. The distribution of arrests for the experimental group and the control groups is non-normal and skewed to the right. Most people either had zero rearrests or 1 rearrest while fewer and fewer people had >1 rearrest. I did however test for homogeneity of variances, which reflected that the variances are equal. **What test do I use in this case?** Also, what tests do I use under the following assumptions?   1. Normal Distribution = True and Homogeneity of Variance = True   2. Normal Distribution = True and Homogeneity of Variance = False   3. Normal Distribution = False and Homogeneity of Variance = True   4. Normal Distribution = False and Homogeneity of Variance = False


title: Singular information matrix error in lrm.fit in R

text: I am running an ordinal logistic regression in R and running into trouble when I include dummy variables. My model works great with my first set of predictors. Next I want to add dummy variables for each of the years represented in my dataset. I created the dummy variables with `car:recode` in this manner (one statement like this for each of the 11 years)               fsd$admityear2000 <- recode(fsd$ApplicationYear ,"2000=1;else=0")      The lrm model is specified as follows               library(Design)     ddist<- datadist(fsd)     options(datadist='ddist')     m4 <- lrm(Outcome ~ relGPA + mcAvgGPA + Interview_Z + WorkHistory_years + GMAT + UGI_Gourman + admityear1999 + admityear2000 + admityear2001 + admityear2002 + admityear2003 + admityear2004 + admityear2005 + admityear2006 + admityear2007 + admityear2008 + admityear2009, data=fsd)      (sorry for all of the other random variables, but I don't want to introduce confusion by changing my code) I get the error               singular information matrix in lrm.fit (rank= 22 ).  Offending variable(s):     admityear2009 admityear2000 admityear1999      Error in lrm(Outcome ~ relGPA + mcAvgGPA + Interview_Z + WorkHistory_years +  :        Unable to fit model using “lrm.fit”      I understand that including all options of a dummy variable over-defines the model, but I get the error whether I include all 11 years or just 10. I found a suggestion here to set the penalty parameter of `lrm` to a _small positive value_. Setting it to 1 or 5 changes the error such that it only names one of the variables as offending. The error doesn't go away even with `penalty=100`. I'm pretty new to R, but loving the freedom so far. Thanks for any help! **Responses and Lessons**   * Factors are awesome and I can't believe I didn't notice them earlier. Man that cleans up my code a lot. Thanks!   * My DV, 'Outcome' is indeed ordinal and after making it a factor(), I also made it ordered().   * The str() command is also awesome and this is what my data now looks like (with some of the non-relevant variables omitted) output:               str(fsd)          Outcome      : Ord.factor w/ 3 levels "0"<"1"<"2"     relGPA       : num       mcAvgGPA     : num      admitschool  : Factor w/ 4 levels "1","2","3","4"     appyear      : Factor w/ 11 levels "1999","2000",..        * both lrm() and polr() now run successfully, and they both deal with appyear by dropping some values of the factor. lrm() drops 1999, 2000, and 2001 while polr() just drops 1999 and 2000. lrm() gives no warnings while polr() says "design appears to be rank-deficient, so dropping some coefs." This is an improvement, but I still don't understand why more than one value needs to be dropped. xtabs shows that there isn't full seperation right? output:               xtabs(~fsd$appyear + fsd$Outcome)                fsd$Outcome     fsd$appyear    0    1    2            1999 1207  123  418            2000 1833  246  510            2001 1805  294  553            2002 1167  177  598            2003 4070  158 1076            2004 2803  106 1138            2005 3749  513 2141            2006 4429  519 2028            2007 6134  670 1947            2008 7446  662 1994            2009 4411   86 1118


title: Plotting a Gaussian in Python

text: I am trying to plot a histogram of my data, and I seem to be a little confused here. I am using matplotlib in Python. Here is the code from their website:               mu = 100 #mean     sigma = 15 #std deviation     x = mu + sigma * np.random.randn(10000)          # the histogram of the data     n, bins, patches = plt.hist(x, num_bins, normed=1, facecolor='green', alpha=0.5)      I am confused as to what the _x_ -axis should be for my use. I have calculated the standard deviation and the mean, but I am uncertain if I should replace the `np.random.randn(10000)` with the actual data.


title: Follow-up question: When should you center your data & when should you standardize?

text: I have a follow up question to MånsT's reply to the "When should you center your data & when should you standardize"-question. ( I cannot leave a comment as I am below the magic "50 reputation".) He says that > In addition to the remarks in the other answers, I'd like to point out that > the scale and location of the explanatory variables does not affect the > validity of the regression model in any way. > > Consider the model $y=\beta_0+\beta_1x_1+\beta_2x_2+…+ϵ $. > > The least squares estimators of $\beta_1,\beta_2$,… are not affected by > shifting. The reason is that these are the slopes of the fitting surface - > how much the surface changes if you change $x_1,x_2$,… one unit. This does > not depend on location. (The estimator of $\beta_0$, however, does.) > > By looking at the equations for the estimators you can see that scaling > $x_1$ with a factor $a$ scales $\beta_1$ by a factor 1/$a$. To see this, > note that > > $\hat{\beta_1}(x_1)=\frac{\sum_{i=1}^n > (x_{1,i}−\bar{x_1})(y_i−\bar{y})}{\sum_{i=1}^n(x_{1,i}-\bar{x_1})^2}.$ > > Thus > > $$ \hat{\beta_1} (ax_1)=\frac{\sum_{i=1}^n > (ax_{1,i}−ax_1)(y_i−\bar{y})}{\sum_{i=1}^2(ax_{1,i}−a\bar{x_1})^2}=\frac{a\sum_{i=1}^n > (x_{1,i}−\bar{x_1})(y_i−\bar{y})}{a^2\sum_{i=1}^n(x_{1,i}-\bar{x_1})^2}=\frac{\hat{\beta_1}(x_1)}{a}.$$ > > By looking at the corresponding formula for $\beta_2$ (for instance) it is > (hopefully) clear that this scaling doesn't affect the estimators of the > other slopes. It is clear to me that $\hat{\beta_1}$ is given by the provided expression above. However, this is the $\beta_1$ in the model: $y=\beta_0+\beta_1x.$, as MånsT has clarified by writing $\hat{\beta_1} (x_1)$. Will this argument still hold in the _general linear model_ , given by $y=\beta_0+\beta_1x_1+\beta_2x_2+…+ϵ $?? Since, in the general linear model, the $\beta$s are given by the normal equations: $\vec{\beta}=(X^t X)^{-1}X^t \vec{y}$ (full rank of $X^t X$ assumed). So: It is not clear to me that the $\beta_1$ in the model with only _one_ predictor will be the same as the corresponding $\beta_1$ in $\vec{\beta}$. According to http://data.princeton.edu/wws509/notes/c2.pdf, pages 21-24, there is a difference (the author makes the distinction between gross and net effects (gross: $\beta_1$ in one predictor model, net: $\beta_1$ in general linear model)). How will the scale-invariance argument be for the general linear model?


title: How to do meta-regression in SPSS?

text: I am trying to manage a meta-regression in SPSS17 using the effect size as the dependent variable. I want to explore if my independent variables affects the effect size. Some small practical questions:   1. What is the minimum number of studies necessary for a meta-regression?     * Some people suggest at least 10 studies are required. Why not 20 or 5 studies?      * Is the total sample size an important consideration?      * Why would 10 studies with 200 patients be enough, but 5 studies with 400 patients not be enough?   2. Can I enter all three regressors at once and report the global model, or do I have to enter one regressor at a time and report 3 models each one separately?      * How does the correlation between the independent variables affects this choice?      * How does the number of the studies affect the number of independent variables that I should enter simultaneously?   3. Does the independent variable have to be a scale variable?     * The dependent variable (=effect size) is, of course, scale. The independent variable must be also scale, or could be ordinal or nominal?   4. How can I weight my effect size for sample size?   5. What is the preferable level of significance?      * Is p<0.05 still acceptable for clinical research in such an analysis?


title: How to express error as a percentage?

text: I have data about how many unique users do a certain thing for each day of a month. I can average it, and i would like to display the variation in a intuitive format (such as % of something). Is there a standard way of doing this? I've found standard error, which is $\frac{\sigma}{\sqrt{n}}$, which is not particularly intuitive for users of the data. If anyone needs clarifications, please ask. I'm not too clear about this myself. EDIT: in response to the answers, it's for building out a analytics dashboard for use by the entire company (so many people probably don't understand standard deviation). In particular, we are doing A/B testing for various metrics over a period of say 1 month. We basically would average all the metric per day to give a number for that period, but there are variations day- to-day and would like to have some good way of expressing that.


title: ARDL, Lag Terms and Singularity

text: I am interested in fitting an ARDL model that has 4 lags for each explanatory variable. However, when I fitting the model in R. R says that coefficients are not defined because of singularities. Is this because I have too many lags ? Further, how should I choose the appropriate lag ? (all variables are stationary)


title: positive skewness in simulation results

text: I am using simulations to make a calculation. I generate many random numbers from a distribution for each input and then I take the mean and standard deviation of the outputs. I noticed that the mean output from the simulations is always slightly higher than the result would be had I not used the distribution. This is due to the distribution becoming more and more positive skewed at each operation. This example illustrates what happens even though I simply used lower and higher bounds instead of random numbers from a normal distribution: 10(+/- 2) * 10(+/- 2) 8 * 8 = 64 (36 lower) 10 * 10 = 100 12 * 12 = 144 (44 higher) The mean of the results is 102.7, not 100. In this example, would I say that 10(+/- 2) * 10(+/- 2) = 102.7? In my actual calculations, the result is fed back into the calculation in order to get an estimate for time 2. After time 10, the result is some 35% higher than it would be without using the distributions. Am I going about this correctly?


title: Change and anomaly detection

text: I recently wrote a program that graphed data points so that a user could scroll through them and find "interesting" parts of the data. Now I am looking at ways to make it even simpler by making a table of values deemed to be interesting. These are usually either the value switching to a new value, or spiking and returning to the old one. The main problem is that I have to take into account a lot of noise. Since I'm already saving hours of data trawling, making it conservative (i.e. is trigger happy) is fine. Algorithms and theories I'm looking at are k-nearest neighbour and Local Outlier factor. I am looking at using both, and then collecting them by saying errors will be unique within +-100 samples, giving errors detected by both a higher rating... Or something along those lines. What algorithms are there I should be looking at? And what would you recommend as the best algorithm? **Edit: More specific information:** The values are read in ~20ms intervals (but to say time is constant is ok), and are between +-10 (as single precision floats). The amount of noise is unknown, but is pretty consistent (i.e. same magnitude of noise throughout), and usually jumps between two values. A typical set of data is 120000+ records long (so records are read 1 at a time, and the whole set, or even it's size, is not known), changes can occur anywhere in that set, and at any frequency, though the change will usually be >2 (though that fact is not reliable) and aren't usually within 1000 records of each other (again, not reliable). Typically the change or spike is <2 records long, and there are few gradual changes (though the possibility of them cannot be ruled out. Changes and spikes are always obvious to a human, and since this is to point a human in the right direction it only needs to be accurate to +- 100 records. At the moment I am looking at something along the lines of: 1 - Get the following from 100 records the mean and standard deviation 2 - For each value test if it is within 2 standard deviations of the mean. 3 - if not then print the change and start from 1 again (ignoring 20 records) That would be simple to implement, and would show gradual changes as being changes. I could add more complex checks to find the type of exception (spike or change) and the 20 doesn't need to be a constant


title: How to make WinBUGS close when "Trap" occurs

text: I'm using `R2WinBUGS` to process several datasets in WinBUGS in "batch mode". Occasionally traps occur, mainly due to bad parameters (I guess), and I managed to avoid the whole R script to stop using `try`, and to retry automatically. Yet, when this happens, I have to manually close WinBUGS in order to let the R script go on. I'd rather find a way to make WinBUGS exit automatically as soon as a trap happens.


title: Difference between the two normal distributions

text: I have two random variables $X$ and $Y$ which follows Normal distribution , whose pdf's are given by $f(x)= \frac{1}{2 \sqrt{2 \pi} \sigma}[e^{\frac{-(x-1)^2}{2 \sigma^2}}+e^{\frac{-(x+1)^2}{2 \sigma^2}}]$ $f(y)= \frac{1}{ \sqrt{2 \pi} \sqrt{1+ \sigma^2}}[e^{\frac{-x^2}{2(1+\sigma^2)}}]$ Both Random Variables $X$ and $Y$ has Mean Zero and Variance $1+\sigma^2$. I am interested to evaluate the Higher moments for both $X$ and $Y$. Is the higher moments for both $X$ and $Y$ are equal. What is the difference between these two Random variables. Thank you.


What is the relationship between propensity scores and sufficient statistics? I just came across the idea.


title: Why is data grouped or binned in this way?

text: I am dealing with an existing algorithm where data is processed as follows: We have a large sample of data taken every minute, this is grouped into 'buckets' of data for every hour. (Each bucket contains a SUM(), MAX(), MIN() and the number of data points taken in the hour.) Now, ideally there should have been 60 samples in each bucket, but the algorithm takes the 60th minute data and adds it to current bucket AND the next bucket. I wanted to know if this is a standard processing technique? Note: The data is prepped for display in a graph later (so is it a smoothing technique?) Note 2: I have very little knowledge of Statistics.


title: Hypothesis Testing on Exponential distributions

text: Let $X_1, \dots, X_n$ be independent exponential $(\theta)$ random variables. Suppose we are interested in testing $H_0: \theta = \theta_0 = 1$ versus $H_A: \theta = \theta_1>1$. Consider two tests with rejection regions: Test 1: {$ (x_1, \dots, x_n) \in \mathbb{R}^n: \sum^n_{j=1} x_j \geq k_{\alpha}$}, Test 2: {$ (x_1, \dots, x_n) \in \mathbb{R}^n: X_{(1)} \geq h_{\alpha}$}. > (1) For $n = 10$ and $\alpha = 0.05$, calculate $k_{\alpha}$. $\bf{My \ thoughts:}$ I know I need to set if up so I find the distribution of $\sum^n_{j=1} x_j $. Not quite sure on setting it up to get a value so I have it equal to the $\alpha$ given. > (2) For $n = 10$ and $\alpha = 0.05$, calculate $h_{\alpha}$. $\bf{My \ thoughts:}$ Same idea as number (1), but again I am stuck on proper set-up. > (3) Using the values found for $k_{\alpha}$ and $h_{\alpha}$, calculate the > power of both tests with $\theta =2, 3, 4$. Based on the results of these > findings, which test should be used? $\bf{My \ thoughts:}$ I think this may require a bit more work but I will just need to set it up according to the cases of the values for $\theta$. Any help is greatly appreciated.


title: Using Q-Q Plot for determining distribution of two data sets

text: Suppose I have some baseline data $X_1$ and some post-baseline $X_2$. For example, the baseline data could be blood pressure. Some intervention is added (e.g. a blood pressure medication) and the post-baseline data is $X_2$ (blood pressure after medication). > **Question.** If I want to determine whether $X_1$ comes from an exponential > distribution and $X_2$ comes from an exponential distribution, does it > suffice to have a single Q-Q plot for the entire data set ($X_1 \cup X_2$)? > Or do we need to make separate Q-Q plots for baseline and post-baseline > data. I have very few observations in the baseline data so it doesn't seem > to make sense to make a Q-Q plot for $X_1$. **Added**. I am using a test that depends on distributional assumptions of $X_1$ and $X_2$. Can I assume that I cannot make distributional assumptions based on the qq-plot? This would imply that I would need to use a non- parametric version of the test.


title: How to set the optimal number of simulations

text: In geostatistical context it is common practice that for simulations of a variable of interest (e.g., grade of concentrations of metal in rock samples) the **number of simulations** at least needs to be **30**. I would ask what _the criteria_ is to choose a right trial number considering that each trial has an expense of large computation. Is that criteria generic enough to be applied on all simulations?


title: What statistical methods to use for analysis of a scenario questionnaire?

text: I've produced a questionnaire with 3 parts, demographics, attitude questions, and then a final section on particular scenarios. I'm looking to analysis people's response to the scenarios - as in what factors affected their choices. I've used Likert scales and also categorical - however, I've been using SPSS so have arranged Labels to Values so that data is nominal. I've used Chi-square tests in contingency tables to look what relationships between answers, and was thinking of using the Kuskal-Wallis test - however, I have thought about factor analysis. I'm not sure what is right or suits my data though! Any help would be greatly appreciated! R


title: Repeated observations on the same unit at the same time point, how to solve xtset repeated time values?

text: I’m interested in testing the impact of gender on personal selling. I have data on each salespeople’s amount of sales per customer for 22 months from Jan 2002 to Oct 2003 (amount of sales per customer $=y)$, each salespeople’s gender $(x_1)$, and a dummy variable (if the product is on sale, $x_2)$. Since gender is time invariant, I opt for the random effects model. Ideally, the model should look like: $y_{ijt}=b_0 + b_1x_{1ijt} + b_2x_{2ijt} + S_j + e_{ijt}$, where $i$ is customer $i, j$ is salespeople $j$, and $t$ is the selling time point, $S_j$ is the salespeople-specific random effects. Now when I specify `xtset` in Stata, I encounter the problem of repeated time values within panel since even on each single day there might be multiple observations for a salespeople. I’m not sure, but can I just specify `xtset salespeople_id` and run `xtreg y i.x1 i.x2 i.month, re` on the following model: $y_{ij}=b_0 + b_1x_{1ij} + b_2x_{2ij} + b\cdot \text{month}_{ij} + S_j + e_{ij}$? Thank you!


title: What is copula transformation

text: I have seen that copula transformation changes my sample space to the range of $[0 \; 1]^d$ where d is the number of dimensions. Can anyone explain me about copula transformation?


title: Compare two (or more?) groups of regression coefficients

text: Short version: I need to compare two small (~10) groups of numbers, usual setup for a non-paired t-test, or Mann-Whitney. But -- the numbers come _each_ with its own SE, and since the groups are small, this is presumably relevant. Question -- how do I take the individual SEs into account? Longer version: There are two groups of patients, treatment and placebo. The test to assess the influence of the treatment involves making many (~100) measurements for a period of time, for each patient. Then an exponential decay curve is fitted into these measurements; the relevant outcome is the time constant of the decay. Since it's estimated from noisy measurements, it has a SE. How do I take the individual SEs into account when comparing the two groups? Artificial example: is the mean of the (approximately estimated with given SEs) numbers in group A different from the mean of group B? $A=\{13\pm2,15\pm4,10\pm2\}$, $B=\{12\pm4,10\pm3,16\pm3,14\pm5\}$. Follow-up: What if there are more than two groups?


title: Why isn't bayesian statistics more popular for statistical process control?

text: My understanding of the bayesian vs frequentist debate is that frequentist statistics:   * is (or claims to be) objective   * or at least unbiased   * so different researchers, using different assumptions can still get quantitatively comparable results while bayesian statistics   * claims to make "better" predictions (i.e. lower expected loss), because it can use prior knowledge (among other reasons)   * needs fewer "ad hoc" choices, replacing them by prior/model choices that (at least in principle) have a real-world interpretation. Given that, I would have expected that bayesian statistics would be very popular in SPC: If I were a factory owner trying to control my process quality, I would primarily care about expected loss; If I could reduce that, because I have more/better prior knowledge than my competitors, even better. But practically everything I have read about SPC seems to be firmly frequentist (i.e. no prior distributions, point estimates of all parameters, many ad-hoc choices about sample size, p-values etc.) Why is that? I can see why frequentist statistics were a better choice in the 1960's, when SPC was done using pen and paper. But why hasn't anyone tried different methods since then?


title: variate transformation

text: logarithmic transformation which is used to stabilize the variance of the distribution of the sample variances , $S^2$. $$S^2\sim^{asy} N(\sigma^2,\frac{2\sigma^4}{n-1})$$ Is that $$log S^2\sim N(log\sigma^2,\frac{2}{n})$$ or, $$log S^2\sim N(log\sigma^2,\frac{2}{n-1})$$ I need online reference of logarithmic transformation.


title: Does this make sense when linear regression results are out of the -1/+1 bandwidth?

text: I've computed a regression on panel data using the `xtreg` command in stata. My coefficients on some of the independent variables (which are binary- variables) are out of the range of -1 +1, even though the dependent variable is log transformed. I find this rather odd, since a result of -2,3 would mean -230%, this doesn't seem possible with purely positive data. Anyone has any idea how to interpret this, or what mistake I could have made?