title: Statistical Test for Analyzing Practice Test Data

text: A problem I've been toying around with: Company A helps a group of students to prepare for a standardized test. A perfect score on the test is 100, but most of Company A's students receive scores between 60-90. To help their students prepare to take the real standardized test, Company A offers a series of 10 practice tests on consecutive weekends. Company A would like to assess whether any of the tests in their practice test line (some of which have been created by the company itself) are, on average, too hard or too easy. Company A has a data set consisting of hundreds of students, all of whom have taken the full series and who eventually take the real test. Since Company A helps students prepare for the test, they expect that, given tests of equal difficulty, students' ability will improve over the 10 tests. Company A does not wish to assume that the improvement is necessarily linear. Here are my questions:   1. What sampling procedure should the company use? Company A wants to sample as few students as possible because accessing old data is time intensive.   2. What statistical tests should the company apply to the data to see whether there is a statistically significant difference in difficulty between tests?   3. How should Company A account for the fact that they expect students' scores to improve naturally over the course of the 10 tests? (This is a slightly concealed version of a problem I ran into at work, and it seemed like a tricky statistical problem. But maybe I'm just not familiar with relevant tools.)


title: Variables for post-stratification weights?

text: What justifies the usage of a variable for post-stratification? I am working with a constituent survey of a non-profit's constituent with 2500 responses out of a much larger sample and even larger population. I have _many_ variables about the target population, which are all active constituents. In literature I've read, it's common to use demographic variables (age, gender, and race, for example), but in my experience with this data, demographics have relatively high data quality errors and weak correlation to non-response error, while behavioral data (for example, donation history) are recorded reliably and correlate better to non-response. I assume demographics are common because many surveys try to get a nationally representative sample, and the government publishes demographic information for this population. Because I have them, is there anything wrong with using the behavioral variables instead of, or in addition to, the demographics? Is there a practical empirical method to choose variables? If the suggestion is to use behavioral variables in addition to demographics, how would I detect or prevent overfitting when raking weights with many variables?


title: frequentist coverage, and besides?

text: One day I gave a $95\%$-confidence interval to a requester who know nothing about statistics. He asked me what does it mean. Roughly, I answered "The population parameter is inside the interval $95\%$ of times". Then he asked "And when it is outside, can it be far from the interval ?". I said no but I was a little confused. What do we know about this ? It would be interesting for instance to know the distribution of the distance from $\theta_0$ to the confidence interval conditionally to the event on which $\theta_0$ is not inside the interval. Moreover could we identify some characteristics of the data sample under which the true parameter is more likely to be to the left or to the right of the confidence interval when it is outside ? (for instance the skewness of the data sample in the case of the confidence interval about a Gaussian mean).


title: Proving Causality with t-test/regression

text: Earlier today I was discussing statistical analysis software with a colleague of mine. My colleague had primarily used SPSS in previous work for performing t-tests, anovas, manovas, and other statistical tests similar to those. My colleague also mentioned that SPSS does not handle regressions well. I have no idea whether or not SPSS handles regressions well, but I did respond that a t-test can be formulated as a regression, so SPSS must not be all that bad with regressions. We got to talking about t-tests, regression, and causality, and it came up that "you cannot prove causality with regression, while t-tests are able to prove causality." I've always thought that causality could be established given an appropriate experimental design, regardless of whether you use a t-test or regression to estimate something or perform a hypothesis test. Is it possible to establish a causal relationship using a t-test and not regression? Is it possible to establish a causal relationship without knowing about the experimental design?


title: Is there such a thing as a multidimensional kappa test?

text: I am involved with designing a new diagnosis method for a routine clinical procedure. I want to compare the agreement between operators to the agreement between methods. Is there such a thing as a multidimensional Kappa test, and if so how do I perform it? EDIT: I could carry out a number of individual Kappa tests, finding the agreement between methods for each operator, as well as the agreement between operators for each method, but there may be some interaction between operators and methods which would be useful to know.


title: Average of predicted values with logistic regression

text: I have a large unbalanced dataset (the target has ~1500x more 0's than 1's) on which I train a logistic regression algorithm to predict the probability of success (Not a binary outcome but a real number between 0 and 1). Some additional details:   * about 100,000,000 lines in my dataset   * I'm using Online Logistic Regression with a Stochastic Gradient Descent.   * I apply a decay factor to my learning rate The average value of the target variable on the training set is about 0.0005. the average of the target variable in the test set is almost exactly the same as in the training set (~0.0005). However, when I test my model on the test set, the overall average of my predictions is a bit out (~0.001). Should I expect from logistic regression to centre the results around the mean? I know that the main purpose of the intercept term is to do so, but what level of difference should I expect? Is the difference in observed vs. predicted mean negligible here?


title: Alternatives to the Baron-Kenny approach to modeling mediation

text: I'm about to open the door to a very thorny issue in the social sciences. How does one correctly model and test hypotheses about mediating variables using observational data? I'm familiar with the Baron-Kenny approach to mediation (see previous answer here), and also with structural equation modeling. However, I've heard both approaches disparaged by more quantitative-minded social scientists than myself -- especially when one is using observational rather than experimental data. So, let's say that I'm trying to resolve the following: **Y** is a behavioral outcome. Both **X** and **Z** are observed characteristics of subjects that cannot be manipulated by an experiment. **X** is an attitude (something that can be changed over the long term) and **Z** is an unchangeable characteristic such as age, race, etc. I hypothesize that **X** is a mediating variable, thus **Z** affects **Y** through the pathway of **X**. While it's reasonable to suggest that **Z** is in some way correlated with **X** , my theory argues that it has no impact on **Y** ( **other than through **X** ). How would one best test these hypotheses using best practices in current research?


title: How do you explain 'co-integration' to determine spurious regression to a fairly new time series student?

text: I would appreciate it if someone could explain the concept of spurious regression intuitively / without it being too technical. So far the student understands that spurious regression is basically when two processes have no relation with each other, but would still conclude that there is one using statistical significance.


title: Inferring speed from non-I.I.D. GPS data

text: I have a collection of GPS log of many vehicles. Each row has timestamp, location, speed and vehicle number. I'd like to infer to distribution of speed, or at the minimum, the portion of time they are stationary v.s. when they are moving. The vendor do not make any claim on the independence of the data. In fact, our theory is that it logs more often when the vehicle is moving than when it is stationary. So the collection of data maybe correlated to speed that we are trying to measure. Is there any valid approach to make inference from this data? Does bootstrapping offer any help for this problem? **Edit:** Here is one view of the data. I plot the speed against the location. Basically I treat them as Bernoulli distribution of red v.s. blue. Just counting the data point it comes out to be 28% blue, which is reasonable. But I want to make a more statistical valid estimate from this. Speed/Distance plot It is a small surprise to see the apparent vertical band. I added jitter on the y-axis but the x-axis is suppose to be continuous. Thanks!!


title: Visualizing probability of event over time, based on epoch time

text: I have a list of epoch/unix times at which an event happened. My hypothesis is that there are certain times during the week when this event might happen more frequently. How can I visualize/determine which times are most probable for this event to happen? I'd like to know both time and day of week. I'm working with excel.               1353951487     1353950051     1353948443     1353919982     1353882762     1353724050     1353693717     1353505109     1353474394     1353415614     1353415429     1353397809     1353383231     1353375862     1353342807     1353298391     1353298644     1353295194     1353250737     etc...      * * * Since epoch time starts on a Thursday morning at 0:00 UTC, I subtracted three days from each time, to shift all times to seconds since Sunday Jan 4th 1970. (since the values listed are epoch seconds) I then found the time modulo (seconds in a week) to find the remainder of seconds which took place after the most recent Sunday. For the times listed above, here are the seconds since the most recent Sunday at 0:00 UTC.               Time        Seconds since most recent Sunday 0:00 UTC     1353951487  149887     1353950051  148451     1353948443  406043     1353919982  377582     1353882762  340362     1353724050  181650     1353693717  151317     1353505109  567509     1353474394  536794     1353415614  478014     1353415429  477829     1353397809  460209     1353383231  445631     1353375862  438262     1353342807  405207     1353298391  360791     1353298644  361044     1353295194  357594     1353250737  313137


title: ARIMA model interpretation

text: I have a question about ARIMA models. Let's say I have a time series $Y_t$ that I would like to forecast and an $\text{ARIMA}(2,2)$ model seems like a good way to conduct the forecasting exercise. $$ \Delta Y_t = \alpha_1 \Delta Y_{t-1} + \alpha_2 \Delta Y_{t-2} + \nu_{t} + \theta_1 \nu_{t-1} + \theta_2 \nu_{t-2} $$ Now the lagged $Y$'s imply that my series today is influenced by prior events. This makes sense. But what is the interpretation of the errors? My prior residual (how off I was in my calculation) is influencing the value of my series today? How are the lagged residuals calculated in this regression as it is the product / remainder of the regression?


title: Conditional heteroscedasticity vs iid hypothesis - McLeod-Li test

text: I am using the function McLeod.Li.test in the TSA package. According to the help file in R, this function “Performs the McLeod-Li test for **conditional heteroscedasticity** (ARCH).” However, it has been mentioned on page 50 of the book “Introduction to Time Series and Forecasting, Brockwell and Davis, 2nd edition” that the McLeod-Li test “can be used as a further test for the **iid hypothesis** ”. Question: Is it true to say that the McLeod-Li test can be used to check the “conditional heteroscedasticity” as well as the “iid assumption”? If not, which one is the correct null hypothesis?


title: Non-fair die - College Probability

text: How many times must you roll a non-fair die to be at least 84% sure that the sample probability will be within 3% from the actual probability. Since the die is not-fair, we do not know p. My question is, can I assume it to be around 1/6, or not? And so far, I found the Z-scores for 84% of the probability from the mean to be + & - 1.4. But, I do not know how to move forward with an unknown probability.


title: Tests for normed vector

text: I have many vectors $a_i\in\mathbb{R}^p:||a||=1,1\leq i \leq n$. I would like to test whether the $a_i$'s are randomly spread out on $S^{p-1}$ (the $p$ variate unit circle). Can anyone point to a test for that? Thanks in advance,


title: What's the difference between statistics and informatics?

text: We always say that statistics is just dealing with data. But we also know that informatics is also getting knowledge from data analysis. For example, bioinformatics people can totally go without biostatistics. I want to know what is the essential difference between statistics and informatics.


title: How should I account for small samples in a larger data set? Should I remove them?

text: We're trying to develop a simple tool that will help our teams optimise their clients' Facebook posting strategies. In our experience, time of post can have a big impact on audience response; but those times will differ from audience to audience based on their behaviours (can they access Facebook from work? What times of day are they online? When are they at a loose end, and when are other activities taking up their concentration?) We expect responsiveness to vary throughout the day for these (and other known & unknown) reasons. Furthermore, of course, time of day isn't the only reason that response rates vary; the content of the post is one very significant reason. I don't know if this background is useful; whether it makes the following question clearer. The following R script (I lack sufficient reputation to post charts) contains data that cover posting activity and audience response by hour of post. You'll notice that the 6 am point (for which we have two observations) massively outperforms the rest of the day. We often see these very high response rates for hours that have fewer observations.               post.hour <- c(0:23)     posts <- c(0,0,0,0,0,0,2,15,16,10,17,13,29,21,23,18,29,24,34,42,51,48,49,17)     response <- c(0,0,0,0,0,0,5282,8627,6080,2716,2831,3258,6291,7756,4008,4614,11838,2611,10527,14706,5416,10970,19098,9505)     mean.response <- response/posts          d <- data.frame(post.hour,posts,response,mean.response)          library(ggplot2)          response.chart <- ggplot(d,aes(post.hour,mean.response)) +      geom_point() + geom_line() +      ylab("mean response") +      opts (title="mean audience response by post hour")     response.chart      I've tried manually removing hours that have fewer than a certain number of observations: This seems unsatisfactory, and (furthermore) is hard to automate across widely differing data sets.   1. Should I remove these hours? If so, what statistical tools should I use to perform this repeatedly and automatically across multiple data sets?    2. If the answer is "no" -- and I suspect that it might be -- what should I do to account for these tiny sample sizes? Conscious that this is possibly Statistics 101 stuff. It should be clear from this that I have not - in fact - had any statistical training. I'd greatly appreciate you taking this into account if and when you respond to this question. I'm not even sure that the hourly data represent samples...


title: What is heavy hitter analysis?

text: I have some data on the number of times each of my machines turned off (due to an error) in a particular time period. There are about 6 different classes of machines being used to construct a total population of 50 machines. I wanted to analyze the stability of the 6 classes of machine relative to each other. An acquaintance suggested that I could do some heavy-hitter analysis during a brief chat to determine if some machines are shutting down more than others. Can someone tell me how to systematically perform this analysis or if there is a formal name for this kind of analysis?


title: Constructing an appropriate aov command in R

text: I have some interlaboratory data where we have tested recovery weights of 3 different batches of product in three labs, with varying numbers of operators in each lab (see below). I want to know whether the labs are all getting similar results. Note that we expect the three different products to have different results. My `aov` line is as below. Is this correct ? EDIT: I want to know if I have specified the model correctly for what I am trying to achieve. Apologies for not making this clear. I.e. can I draw from this that all three labs are getting similar results ?               wt<-read.csv(header=T, text="     Date,Lab,Operator,Sample,Weight     05/11/13,L1,Op1,Prod1,0.42     07/11/13,L1,Op2,Prod1,0.58     08/11/13,L1,Op3,Prod1,0.95     30/10/13,L2,Op5,Prod1,1.77     30/10/13,L2,Op4,Prod1,1.02     31/10/13,L2,Op5,Prod1,1.67     01/11/13,L2,Op4,Prod1,0.73     05/11/13,L1,Op1,Prod2,0.06     07/11/13,L1,Op2,Prod2,0.18     08/11/13,L1,Op3,Prod2,0.42     30/10/13,L2,Op5,Prod2,0.66     30/10/13,L2,Op4,Prod2,0.34     31/10/13,L2,Op5,Prod2,0.24     01/11/13,L2,Op4,Prod2,0.73     05/11/13,L1,Op1,Prod3,0.12     07/11/13,L1,Op2,Prod3,0.14     08/11/13,L1,Op3,Prod3,0.3     30/10/13,L2,Op5,Prod3,0.8     30/10/13,L2,Op4,Prod3,0.04     31/10/13,L2,Op5,Prod3,0.2     01/11/13,L2,Op4,Prod3,0.08     01/01/14,L3,Op6,Prod3,0.16     01/01/14,L3,Op6,Prod1,1.23     01/01/14,L3,Op6,Prod2,0.10")          summary(aov(Weight ~ (Sample * Lab) + Error(Operator/(Sample*Lab)),data=wt))      (note in real life I have also done outlier tests and removed some samples).


title: Sanity check: how low can a p-value go?

text: I am using a rank sum test to compare the median of two samples (n=120000) and have found that they are significantly different with: p=1.12E-207. Should I be suspicious of such a small p-value or should I attribute it to the high statistical power associated with having a very large sample? Is there any such thing as a suspiciously low p-value? Feedback appreciated.


title: How to derive the Fisher information in the laplace approximation of a generalized linear mixed model?

text: I am currently using the Laplace approximation to fit some geostatistical models for binomial data. Regarding parameters estimation I do not have any problem. I can easily implement the Laplace approximation and get my estimates. But when it comes to obtain the standard errors I am puzzled about how I can derive the hessian of my objective function, hence the Fisher information. I tried to derive this analytically but it seems very laborious and tedious. Is there any reference about this or someone has a solution? Thanks for all your answers in advance. EG


title: Using Cronbach's alpha to test reliability of a questionnaire

text: I administer a questionnaire to managers asking them to what extent they agree / disagree that they perform certain functions (listed in the questionnaire) on a scale of 1-5. The sample functions are   * Setting Work Objectives    * Providing Resources    * Communicating    * Organizing    * Controlling Should I use Cronbach's alpha to test the reliability of this questionnaire?


title: Clustering method that can use graph links, discrete and continuous properties?

text: I have an un-weighted, directed graph that clusters ok using MCL or other graph clustering algorithms. However, I also have additional discrete and continuous properties of the nodes being clustered that I would like to use to get more easily interpretable clusters. What are some clustering algorithms or tools that can work with all three kinds of data: graph, discrete, continuous?


title: How to predict how much data to collect?

text: I'm currently in the planning stages of an experiment that I hope will detect a 0.155% signal variation (relative magnitude) within a reasonable time frame (less than 6 months ideally). I've calculated the rate of (usable) data will be around 68 events per day, though it should be stressed this is a random variable. Now I'm trying to work out **how many days will I need to run the detector for to see the variation with a confidence level of $3\sigma$?** Some other details that may (or may not) be relevant include: the variation in the signal is expected to be sinusoidal with a period of 0.5 days. For this reason I reduced my useful event rate to 34 (Ie half) as clearly there is no variation to see when the sinusoidal signal is at or close to the mean value. I've been googling for a method to predict the size of a data set necessary to see such a small signal variation but have come up with nothing. I would be extremely grateful for any hints / tips anyone could offer. **EDIT:** With the greatest thanks to @shabbychef I feel this question has now been satisfactorally answered on physics.SE here.


title: How do you call this dynamic sample-size selection strategy?

text: Imagine that you want to assess the compressibility of a large document very fast. You could randomly pick a subsequence, try to compress it. This can serve as a prediction for the overall compressibility of the document. But how big should your sample be? We have come up with the following strategy:   1. Pick an arbitrary (small) sample size. Measure the compressibility.    2. Next, double the sample size and measure the compressibility again. If there is little change (say less than 10%), then conclude that you have reliably determined the compressibility of the document. If not, double the sample size again, and so on. We are quite certain that this is not a new strategy, and we are wondering whether it is related to some well-known strategy used by statisticians. ("Compression" here is just an example. Basically, we are interested in a metric that has no known nice mathematical properties, so that it is not possible to determine analytically what could be a good sample size. We have no choice but to fall back on such heuristics.)


title: Exponentially decaying variable

text: I have a variable which decays very quickly (exponentially). I want to transform it in a way it maintains "exponentially decaying" property but which will allow me to control the "speed" of a decay.               original <- c(100, 80, 70, 65, 64, 63, 63, 63, 63, 62, 62, 61);     plot(original);     target <- c(100, 90, 70, 60, 65, 64, 64, 63, 63, 63, 63);      target is decaying much slower. I want to be able to controll for that. How do I do that?


title: What R packages are available for nonparametric regression of two predictors?

text: I would like to fit a nonparametric regression model with two predictors. What R packages are available? I found 'loess' function in 'stats' package and 'gam' function in 'gam' package. Is there anything else?


title: Understanding output stepAIC

text: I am using the `stepAIC` function in R to do a bi-directional (forward and backward) stepwise regression. I do not understand what each return value from the function means. The output is:                         Df     Sum of Sq    RSS       AIC     <none>                        350.71   -5406.0     - aaa      1     0.283        350.99   -5405.9     - bbb      1     0.339        351.05   -5405.4     - ccc      1     0.982        351.69   -5400.5     - ddd      1     0.989        351.70   -5400.5      Question Are the values listed under `Df`, `Sum of Sq`, `RSS`, and `AIC` the values for a model where only one variable would be considered as the independent variable (i.e. y ~ aaa, y ~ bbb, etc.)? **NOTE: Question also in data science exchange:** data science link


title: Linear Mixed Model: Specification of matrices

text: I would like to conduct a power analysis for a linear mixed model with fixed effects for Treatment (Two levels) and Time (four time points: pre, mid, post treatment, 3 months post treatment) and correlated random effects for children (intercepts and slopes). I am using the following code in R using the `lmer` function:               library(lme4)     expdat <- expand.grid(kid = factor(1:500), Time = factor(1:4), Treat = c("XTx", "BAU"))     expdat$obs <- factor(seq(nrow(expdat)))          set.seed(101)     nsim <- 20     beta <- c(100, -7, 8, 15, 20, 0, 0, 0)     theta <- c(15.000000, 7.500000, 7.500000, 7.500000, 12.990381, 4.330127,                 4.330127, 12.247449, 3.061862, 11.858541)     ss <- simulate(~Treat*Time + (1+Time | kid), nsim = nsim, family = gaussian,                     weights = rep(25, nrow(expdat)), newdata = expdat,      newparams = list(theta = theta, beta = beta, sigma = 1))     expdat$Outcome <- ss[, 1]     fit1 <- lmer(Outcome ~ Treat*Time + (1+Time | kid), data = expdat)      This is a repeated measurement design where each person is measured four times. To simulate the data assuming random intercepts and random slopes and their correlations, shouldn't I specify the covariance matrix for the random effects? How would one do that in R?


title: Logistic regression: maximum likelihood vs misclassification

text: Traditionally the fitting of the logistic regression function is explained using maximum likelihood. Could one fit the logistic regression function as well based on either least-squares or based on misclassification error, or is this not possible / has drawbacks?


title: Converting data from classical scaling in MDS

text: I implemented the classical scaling and obtained the relative coordinates in Java.   I want to translate these relative positions, but I am not getting how I should be doing this. This is just like:               T (x, y) =      (x y) (0.981 −0.196            0.196 0.981) + (0.016 − 0.835)       There is a 2 × 2 matrix inside. I want to implement it in Java, and it is a part of MDS-MAP. Does anybody have an idea about this?


title: Toastmasters speech slot availability: what to measure?

text: I belong to a new Toastmasters club that is growing quickly. We meet twice a month and each meeting has four opportunities (slots) for members to give prepared speeches. As more members vie for the available slots, we're looking at ways to respond to the demand. I figure the first step is to measure the availability of the slots. I've automated the process of recording all available slots, polling each day. Now I'm trying to come up with a way to present the availability. My first pass is a graph showing number of days until the next free slot. But this needs refinement, and I'd appreciate your help. One challenge is that there are effectively more free slots than it appears, since speakers sometimes give up their slot in the week before the meeting (schedule conflict or not prepared to speak). That would make my simple graph look weird.   * today: next slot in 49 days    * tomorrow: next slot in 8 days I've also looked at measuring total number of slots available over some period… but that doesn't feel like enough guidance. I've tried finding some prior art on meeting room scheduling, which would seem similar, but no luck so far. Thanks for any ideas you have.


title: Mode of Joint Posterior - Maximization Problems

text: I have a problem whereby I get two different answers if I try to maximize a function. let $ \begin{bmatrix} Y_{o}\\ Y_{a} \end{bmatrix}|\phi\sim N (0,\phi^{-1}A^{-1}) $ $\pi(\phi)=\frac{1}{\phi}$, where $Y_{o} \in \mathbb{R}^{n_{o}}$, $Y_{a} \in \mathbb{R}^{n_{a}}$ I am trying to find the mode of the joint posterior $f(Y_{a},\phi|Y_{o})$ **Method 1** $f(Y_{a},\phi|Y_{o})=f(Y_{a}|\phi,Y_{o})f(\phi|Y_{o})$. $f(Y_{a}|\phi,Y_{o})$ is normal and its mode is independent of $\phi$, so $Y_{a}^{*}=-A_{aa}^{-1}A_{ao}Y_{o}$. $\phi^{*}$ is then the mode of $f(\phi|Y_{o})$ , $\phi^{*}=\frac{n_{o}-2}{Y_{o}^{T}A_{oo.a}Y_{o}}$ **Method 2** Maximizing the log of the joint. $\log f(Y_{a},\phi|Y_{o})=$ $-\frac{\phi}{2} Y_{o}^{T}A_{oo.a}Y_{o} -\frac{\phi}{2}(Y_{a}+A_{aa}^{-1}A_{ao}Y_{o})^{T}A_{aa}(Y_{a}+A_{aa}^{-1}A_{ao}Y_{o}) +\frac{n_{o}+n_{a}-1}{2}\log \phi$ yields $Y_{a}^{*}=-A_{aa}^{-1}A_{ao}Y_{o}$ and $\phi^{*}=\frac{n_{o}+n_{a}-2}{Y_{o}^{T}A_{oo.a}Y_{o}}$ **Modal $\phi^{*}$ differ by $n_{a}$** I don't know why they differ by $n_{a}$, where is my error?


title: Bayes Interval Estimator using Neyman-Pearson

text: This problem is 9.56 from Casella-Berger. Let $X \sim f(x|\theta)$ and suppose that we want to estimate $\theta$ with an interval estimator $C$ using the loss $$L(\theta, A) = bLength(A(X)) - I(\theta \in A(X))$$ If $\theta$ has the prior pdf $\pi(\theta)$, show that the Bayes rule is given by $$ C^{\pi} = \left\{ \theta: \pi(\theta|x) \geq b\right\}$$ The hint given is to write $Length(C) = \int_{C} 1 d\theta$ and use the Neyman-Pearson lemma. After working on this, I believe I have shown the desired result but did not apply the Neyman-Pearson lemma. I am hoping to see this problem solved using Neyman-Pearson, as that may be the point of the problem. Thanks for any help.


title: How to calculate the ratio of "data used in analysis" and "saved data for validation" from a sample set?

text: Provided a sample size S that I plan on using to forecast data. What are some of the ways to subdivide the data so that I use some of it to establish a model, and the remainder data to validate the model? I know there is no black and white answer to this but it would be interesting to know some "rules of thumb" or usually used ratios. I know back at university, one of our professors used to say model on 60% and validate on 40%.


title: How to represent an unbounded variable as number between 0 and 1

text: I want to represent a variable as a number between 0 and 1. The variable is a non-negative integer with no inherent bound. I map 0 to 0 but what can I map to 1 or numbers between 0 and 1? I could use the history of that variable to provide the limits. This would mean I have to restate old statistics if the maximum increases. Do I have to do this or are there other tricks I should know about?


I have computed percentage change from time1 to time2 for several variables. Can I predict percentage change in earnings from percentage change in produced and percentage changed in price? When I ran a model with actual data and dummy coded time (time1=1, time2=0), the dummy variable was not statistically significant. But there are large changes.


title: Welch's t-test gives worse p-value for more extreme difference

text: Here are four different sets of numbers: A = {95.47, 87.90, 99.00}   B = {79.2, 75.3, 66.3}   C = {38.4, 40.4, 32.8}   D = {1.8, 1.2, 1.1} Using a two-sample t-test without assuming equal variances, I compare B, C, and D to A and get the following p-values: 0.015827 (A vs B)   0.000283 (A vs C)   0.001190 (A vs D) I find it strange that the p-value from the A-D test is worse than the A-C test: the difference between the means is clearly much bigger AND the variance of D is much lower than the variance of C. Intuitively (at least for my intuition), both these facts should drive the p-value lower. Could someone explain if this is a desired or expected behavior of the t-test or whether it has to do more with my particular data set (extreme low sample size perhaps?). Is the t-test inappropriate for this particular set of data? From a purely computational point of view, the reason for a worse p-value seems to be the degrees of freedom, which in the A-D comparison is 2.018 while it is 3.566 in the A-C comparison. But surely, if you just saw those numbers, wouldn't you think that there is stronger evidence for rejecting the null hypothesis in the A-D case compared to A-C? Some might suggest that this is not a problem here since all p-values are quite low anyway. My problem is that these 3 tests are part of a suite of tests that I am performing. After correcting for multiple testing, the A-D comparison doesn't make the cut, while the A-C comparison does. Imagine plotting those numbers (say bar-plots with error bars as biologists often do) and trying to justify why C is significantly different from A but D is not... well, I can't. **Update: why this is really important** Let me clarify why this observation could have a great impact on interpreting past studies. In bioinfomatics, I have seen the t-test be applied to small sample sizes on a large scale (think differential gene expression of hundreds or thousands of genes, or the effect of many different drugs on a cell line, using only 3-5 replicates). The usual procedure is to do many t-tests (one for each gene or drug) followed by multiple testing correction, usually FDR. Given the above observation of the behaviour of Welch's t-test, this means that some of the very best cases are being systematically filtered out. Although most people will look at the actual data for the comparisons at the top of their list (the ones with best p-values), I don't know of anyone who will look through the list of all comparisons where the null hypothesis wasn't rejected.


title: Range features with SVM

text: I have some data that I'd like to apply classification to:               Label    Start1, End1,   Start2, End2,    Start3, End3         0       range 1         range 2          range 3         1 ...         2 ...         1 ...         ...      SVM seems like a good approach for this particular type of data, but I'm not very confident with selecting features. The natural features are the ranges above, expressed as a start and end value, and classification seeks to find the label for particular value that falls within a range (any range). Should there be 6 features, or would it perhaps be better to use Start and Length? Thanks


title: How does the mice imputation function work?

text: I was wondering if anyone had experience using the mice function, as described in mice: Multivariate Imputation by Chained Equations in R (JSS 2011 45(3))? I have a dataset with a number of variables, each with varying degrees of missing data. My primary question is: say I use Bayesian linear regression to impute missing data, does `mice` automatically use predictor variables from most significant to least significant to impute? Also, is it common to perhaps average all the imputed datasets?


title: Correlated brownian interpolation

text: I would like to generate conditional correlated random variables. I have a correlation matrix between normal variables, and these variables are modeled through SDEs. What are the algorithms to generate such random values? One example I found is; http://www.mathworks.com/help/toolbox/econ/interpolate.html


title: Nadaraya Watson model

text: I am working through Bishop's book PRML. I have been stuck on question 6.18, for a while and would appreciate some guidance. I feel like I am fundamentally not understanding something. ![enter image description here](http://i.stack.imgur.com/MgKDe.png)


title: Convergence of a genetic algorithm

text: Does anyone know of any method for deciding when a genetic algorithm is done? In MCMC (e.g, BUGS), several chains are started at different, random points. When they all look the same, it is done. Has this approach ever been tried with GA? Any other ideas? thanks


title: R from command line

text: I have a question. Here it is guys. > How do you use R from the command line?


title: How to analyze data with restricted range proportions

text: I have a dataset where one variable is a proportion with a restricted range (say 4%). I would like to use this variable to predict a death rate. Running a simple correlation I'm guessing would underestimate the relationship (?). I also have raw counts for both variables. What analysis do I do here? I'm stuck. Binomial regression? Below is an excerpt of my early attempts.... Firearm denial proportion and firearm suicide rate are positively correlated (r = .25, p < .01). There is a range restriction problem with using denials as a proportion (R = 4.13), but that should only result in underestimating the correlation. Thanks!


title: Choosing correct C and g parameters for libsvm

text: libsvm 3.18     Features: 10      I have used following, parameter range:               C-parameters: 2^-3  , 2^-2  .. 2^11,     g-parameters: 2^-11 , 2^-10 .. 2^-3,      with `10-fold` cross validation(-v 10). Maximum `accuracy(81.8646%)` at `c-value 2048` and `g-value 0.125`. When I used default c and g: 80.9531% accuracy is noticed. Question 1: If high values of C and g parameters caused any over-fitting and increased accuracy ? Question 2: Should I change the grid size for rescanning ? Question 3: Can I proceed with the c and g of 2048 and 0.125 ? Thanks in advance.


title: post hoc test other than Tukey of factorial Anova data in R

text: I am looking for help on post-hoc test of my group data (treatment and stage and interaction) after 2 way anova in R software. The data shown below is a example only actually my data shows significant result after ANOVA test. When I tried post hoc with TukeyHSD it gives pair wise comparisons which is a long list. But I am looking for other test where I can compare treatment mean alone, stage mean alone as well as stage*treatment with letter after mean showing significant or non-significant in group.Could there be other post-hoc text which fulfil my need? My data looks like Treatment Stage Chlorophyll Salt Green 0.2 Salt Green 0.3 Salt Green 0.4 Salt Pink 0.5 Salt Pink 0.3 Salt Pink 0.2 Salt Red 0.5 Salt Red 0.6 Salt Red 0.7 Nitrogen Green 0.4 Nitrogen Green 0.6 Nitrogen Green 0.9 Nitrogen Pink 0.2 Nitrogen Pink 0.3 Nitrogen Pink 0.5 Nitrogen Red 0.4 Nitrogen Red 0.2 Nitrogen Red 0.3 Control Green 0.5 Control Green 0.6 Control Green 0.8 Control Pink 0.5 Control Pink 0.4 Control Pink 0.6 Control Red 0.2 Control Red 0.3 Control Red 0.1


title: Smoothing strategies for features assuming values from countably infinite domains

text: I am in the midst of programming a simple Naive Bayes classifier as an exercise. It is supposed to perform word-sense disambiguation on natural language phrases, e.g. predicting the correct meaning of the ambiguous word _bass_ in a sentence like _I love the taste of bass_. At the moment I am struggling with handling discrete features potentially assuming one of infinitely many values (i.e. features assuming values from countably infinite domains). An example for such a feature might be the length of the preceding word. In the example above, the feature value would be 2, as the word preceding _bass_ has this many characters. The root of my problem is that in some languages (maybe not in English, but certainly in other languages as for instance German) one can form new words of arbitrary length by concatenating an arbitrary number of already existing words (the linguistic term for this process is _compounding_ ). For example, the three German words _Fußball_ (football), _Schuh_ (shoe) and _Schnürsenkel_ (laces) make up the valid compound _Fußballschuhschnürsenkel_. That is, the domain of the feature _"Length of the preceding word"_ is the set of natural numbers (including 0 for observations where the target word is the first word in the sentence), i.e. it is countably infinite. If I decided to use this feature, the classifier during training would have to determine the distribution for the probability $P(L_p = l_p|W = w)$. That is to say, given a target word $w$, how probable is it that the length of the previous word is equal to $l_p$. Afterwards, when classifying input sentences or phrases, it can happen that due to data sparseness no probability for a certain $(L_p, W)$ pair could be determined during training. For example, it might happen that in an input sentence the target word ( _bass_ in the above example) is preceded by a very long word (e.g. 30 characters), but that in the training data no such cases were present, i.e. the probability $P(L_p = 30 | W = bass)$ would be zero. Facing the issue of zero values described above, it is necessary to avoid them by using an appropriate kind of smoothing. Obviously, not all established smoothing methods are applicable in such a case. For instance, the widely used add-one smoothing which transforms a probability $n_{F=v}/N$ into $(n_{F=v}+1)/(N+B)$ (where $n_{F=v}$ is the number of observed situations in which feature $F$ assumes the value $v$, $N$ is the total amount of observations and $B$ is the number of values that feature $F$ can assume) would fail in this case as $B$ would be $\infty$ (since the number of values $F$ can assume is unbounded because the preceding word can have any length). The considerations above, finally, lead me to the following question: What smoothing methods are there to deal with features assuming values from countably infinite domains?


title: R: subset data frame by logical vector derived from column groups

text: I have a data frame like this:                 period     x   y db     perc     1 2013-08-26 4 166 nh 2.409639     2 2013-09-02 5 222 nh 2.252252     3 2013-09-09 3 223 nh 1.345291     4 2013-09-16 9 198 nh 4.545455     5 2013-09-23 3 213 nh 1.408451     6 2013-09-30 5 226 nh 2.212389     ...      There are many observations per db. I wish to make a subset of the dataframe excluding any rows where the sum of the x values for that db == 0, and have derived a logical vector `removal_candidates` like this:               xsums <- by(drawdata$x, drawdata$db, sum)     removal_candidates <- xsums == 0     ...             ne    nf    nh    ni    nj    nl    nm    nn    no    np    nq    nr    nu      FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE      In this example I'd like to make a copy of drawdata without the rows for **nf** but naive attempts such as `drawdata[drawdata$db == removal_candidates]` doesn't work.


title: How to use the hier.part package for hierarchical partitioning?

text: Hello after struggling with using R for the last couple of days I was hoping someone could help me with a statistical analysis I am completing for an environmental science honours project. Using R statistics is not something we have been taught and I am worried that I may have bitten of more then I can chew, however my whole project is based around the **hierarchical partitioning method and the exhaustive search multiple regression analysis method.** The hier.part package was installed along with gtools. I have converted my dataset to a .csv file with seven independent variables and one dependant variable with around 400 replicates (my intention is to do this analysis on eight datasets in total with different amounts of replicates and another dependant variable, but I am starting with this one). The dependant variable is GPP, the independent variables are, NDVI, Temperature, Precipitation, Solar Radiation, Nutrient Availability and Soil Available Water Capacity. Secondly I imported the .csv file into R using the script               GPPANDDRIVER <- read.table("C:\\etc, header=T, sep=",")      This works fine and I can edit the table using               edit(GPPANDDRIVER)      After looking at the `hier.part` package documentation available here it seems like I need to define Y which in the script below is the dependent variable and define `scan` which is the independent variables (mentioned before).               hier.part(y, xcan, family = "gaussian", gof = "RMSPE", barplot = TRUE)      I was defining the dependant `y` vector as               y <- as.vector(GPPANDDRIVER["GPP"])      This also works fine and I have my y vector. However I am not sure how to load independent variables onto the xcan dataframe part of the script. I have tried typing in two scripts but they have not worked.               xcan <- as.vector(GPPANDDRIVER[-GPP])     ## AND     xcan <- data.frame(GPPANDDRIVER[-GPP])      If anyone could help me find the right script for representing my independant variables as xcan that would be greatly appreciated. Also once defined if I entered in the hier.part script mentioned above would R then show me results of the analysis after processing? I will be moving onto to the regression analysis after this if anyone can shed some light on this first problem.               *information on hier.part arguments.*          **Arguments**          y a vector containing the dependent variables          xcan a dataframe containing the n independent variables          family family argument of glm          gof Goodness-of-fit measure. Currently "RMSPE", Root-mean-square ’prediction’          error, "logLik", Log-Likelihood or "Rsqu", R-squared          print.vars if FALSE, the function returns a vector of goodness-of-fit measures. If TRUE, a data frame is returned with first column listing variable combinations and the     second column listing goodness-of-fit measures.


title: Bootstrap, Monte Carlo

text: I have been set the following question as part of homework: > Design and implement a simulation study to examine the performance of the > bootstrap for obtaining 95% confidence intervals on the mean of a univariate > sample of data. Your implementation can be in R or SAS. > > Aspects of performance you may want to look at are confidence interval > coverage (i.e., what proportion of times the confidence interval contains > the true mean) and Monte Carlo variation (i.e., how much do the upper and > lower confidence limits vary between simulations)' Does anyone know how to go about the Monte Carlo variation aspect of this? I can't seem to even work out an algorithm or anything. Is it to do with Monte Carlo integration? Thanks!


title: Variable Clustering Analysis

text: I have a data set that consists of 143 variables (~11000 observations), and I wish to do variable clustering to reduce the dimension. I am using `hclustvar` function in R from `ClustOfVar` package. After performing the hierarchical clustering using `tree <\- hclustvar(X)`, I wish to use the `stability` function to find out the optimal number of clusters. The following is the result: ![enter image description here](http://i.stack.imgur.com/GjTCc.png) Any help is highly appreciated Regards


title: How to assess differences between three groups in a pretest - posttest design using SPSS?

text: I have a dataset that includes one independent variable with 3 levels (group 1, group 2, and group 3) and a repeated dependant variable (pre test score, post test score). I am using the statistics package SPSS and when I took the method I thought I would use "repeated measures anova" there didn't seem to be a way of splitting the independent variable into the 3 groups. In the data view of the spreadsheet each row represents a participant and so they are grouped by including a numerical value to each group 1 = group 1, 2 = group 2, and 3 = group 3. I need to compare the mean changes/improvements over the 2 time points (pre, and post) across the groups to assess if there is a significant improvement in any of the groups. Do you know the most appropriate approach to do this? The only option I am considering is creating a new spreadsheet where each row represents a participant and each column represents the difference between pre and post trials, and then entering this into a one way independent anova.


title: Which algorithmic approach is best for this problem

text: I have a problem where I get data from 3 different sources **s1, s2, s3** , and I have the target ( **value** ). I might have missing values from some of the predictor variables at certain rows. This is the structure of my data set:               s1 s2 s3 | value     --------- ------     34 23 34 | 34     *  13 12 | 12     76 12 75 | 75     ...      At prediction time, I will get data from some of the sources, and I need to predict a value for the target. I am confused about how to start. How should I best go about formalizing this problem--is this a classification or a regression problem? Can I formulate it as a recommender system?


title: Overall Significance vs. Individual Variable Significance in Mutliple Regression

text: I'm running a statistical analysis for work in which I'm trying to determine which (if any) key economic indicators influence our sales. Here is the summary data that I'm getting when I run a multiple regression: ![enter image description here](http://i.stack.imgur.com/GQseA.jpg) It appears to me that the overall regression has strong significance, as does the HousingStarts variable. The R-squared and Adjusted R-squared values also appear to indicate strong correlations. My question is on how I should interpret the ConsumerSpend variable, which has a somewhat high p-value. Does that indicate a weak correlation to the overall regression and therefore I should disregard it, or does the overall p-value trump the individual variable p-value? I should note that I ran separate linear regressions for the two variables, and both of the p-values were well below .001 (although the R-squared and adjusted R-squared for the ConsumerSpend regression was .2285 and .2163, respectively). Sorry if this is a basic question, but I took one stats class in business school, and I'm trying to apply what I learned!


title: Estimation of regression with autocorrelated errors

text: In a book it is written that, In regression work we typically assume that the observational errors are pairwise uncorrelated. But in most time series data , the successive residuals have tendency to correlated with themselves.If we find autocorrelated errors, we need to modify the regression procedure to remove the effect of autocorrelated errors. Typically this is done by an appropriate transformation of the variables used in the regression estimation procedure. Consider a multiple regression model with autocorrelated errors $$Y_t=\beta_0+\beta_1X_{1t}+\beta_2X_{2t}\ldots+\beta_pX_{pt}+\epsilon_t\ldots(1)$$ and the same regression model at time ${t-1}$ is $$Y_{t-1}=\beta_0+\beta_1X_{1,{t-1}}+\beta_2X_{2,{t-1}}\ldots+\beta_pX_{p,{t-1}}+\epsilon_{t-1}\ldots(2)$$ Multiplying both sides of this equation by correlation coefficient $\rho$,the correlation between adjacent errors gives $$\rho Y_{t-1}=\rho \beta_0+\rho \beta_1X_{1,{t-1}}+\rho \beta_2X_{2,{t-1}}\ldots+\rho \beta_pX_{p,{t-1}}+\rho \epsilon_{t-1}$$ Then we subtract this equation from the first equation to obtain $Y_t-\rho Y_{t-1}=\beta_0(1-\rho)+\beta_1(X_{1t}-\rho X_{1,{t-1}})+\beta_2(X_{2t}-\rho X_{2,{t-1}})\ldots+\beta_p(X_{pt}-\rho X_{p,{t-1}})+u_t\ldots(3)$ where $u_t=\epsilon_t-\rho \epsilon_{t-1}$ has uniform variance and is not autocorrelated. ### $\bullet$ My first question is : Why do they multiply $\rho $ in model $(2)$? Why not in model $(1)$? ### $\bullet$ My second question is :Are the data of independent variable $X$ in model$(1)$ same as the data of independent variable $X$ in model $(2)$ ? If this data depends on time, then they seem to be not equal data. If so then my parameters $\beta_0,\beta_1,\ldots\beta_p$ are not same for the model $(1)$ and $(2)$. Then how can i write model $(3)$? More specifically, how can i take $\beta_0,\beta_1,\ldots\beta_p$ as common factor ? ### $\bullet$ My last question is : how does $u_t=\epsilon_t-\rho \epsilon_{t-1}$ have uniform variance and why is not autocorrelated?


title: Datasets for market segmentation

text: What are the gold standard datasets to test market segmentation algorithms with? I'd like to try a few algorithms on known datasets for comparison before I try my own dataset.


title: Multicollinearity over intervals

text: You have multicollinearity when you have 2 variables $(X_1, X_2)$ that have a relationship, $X_1=a+X_2$ where $a$ is constant. My question is: is there still a multicollinearity issue if $a$ is not constant, $X_1=X_3+X_2$? For example your testing children and you have their birthday and the date of the test, so `date of test` $-$ `birthday`$=$`age at test`. Could you use age at test and birthday in a regression model? I ask because there would be a perfect correlation for one test interval, but as you add test intervals the correlation decreases.


title: Distribution of distance from center of sample group

text: We have a bivariate normal process where $X, Y \sim N(0, \sigma)$, with no covariance. (For convenience we can assert that $\sigma = 1$, or that we have a good estimate for its value.) **What is the distribution of the random variable $R(n) = \sqrt{\overline{x_i}^2 + \overline{y_i}^2}$ -- i.e., the Euclidean distance between the sample center of a _n_ points and the true center (at the origin)?** Note that as defined:   1. $R(n) \ge 0$   2. $E[R(n)] \to 0$ monotonically as $n \to \infty$   3. The Rayleigh distribution gives us $R(1) = \sigma \sqrt{\pi/2} \approx \sigma 1.25$ Furthermore, based on a Monte Carlo simulation for $n \in [2, 25]$ with $\sigma = 1$:   1. Variance decreases monotonically as _n_ increases   2. Skewness appears constant across _n_ at 0.63   3. Kurtosis appears constant across _n_ at about 0.24 (This question is a simplified version of this slightly more complicated one that seems to have gotten derailed in complications.)


title: Hypothesis testing exogenous variables among simultaneous equations

text: Consider the model $q_t^d\;=\;\alpha_0+\alpha_1p_t+\alpha_2y_t+e_t\\ q_t^s\;=\;\beta_0+\beta_1p_t+\beta_3z_t+u_t\\ q_t^d\;=\;q_t^s$ where $p_t$ is endogenous and $y_t$ and $z_t$ are exogenous. Suggest how to test $H_0\;:\;\alpha_2=\beta_3=0 \\ H_A\;:\;\alpha_2\neq0,\:\beta_3\neq0$ I combined both equations, wrote everything in terms of $p_t$, and combine all the variables into a matrix. $p_t\;=\;X\gamma +\epsilon$ where $X=(1|y|z)$ and then was trying to set up a test from there. I don't know if this makes sense. Some hints would be very helpful. ALSO: if someone could go over the rank and order conditions for identification, that would be really helpful. I know rank condition is necessary and sufficient while order condition is necessary--but I don't know why this is the case.


title: Cross Validation for mixed models?

text: My colleague and I are fitting a range of linear and nonlinear mixed effect models in R. We are asked to perform cross-validation on the fitted models so that one can verify that the effects observed are relatively generalizable. This is normally a trivial task, but in our case, we have to split the whole data into a training part and a testing part (for CV purposes) that share no common levels. For example, The training data may be based on Groups 1,2,3,4; The fitted model is then cross-validated on Group 5. So this creates a problem since the group-based random effects estimated on the training data do not apply to the testing data. Thus, we cannot CV the model. Is there a relatively straightforward solution to this? Or has anyone written a package to tackle this problem yet? Any hint is welcome! Thanks!


title: Advice on classifier input correlation

text: I have a lot of data on previous race history and I'm trying to predict a percentage chance of winning the next race using Regression, kNN, and SVM learning algorithms. Say a race has 5 runners, and each runner has a previous best course time of, say $T_i$ (seconds). I've also introduced an additional input for RANK of previous best course time of the 5 runners with value 0 to $1 - \frac{T_i-T_{min}}{T_{max}-T_{min}}$ My question is: does introducing both the absolute best course time and rank best course time cause any problems? I understand that these inputs are likely to correlate but if someone runs a world record time they are more likely to win easily but this will get lost using the rank input only which would assign them a rank of 1.


title: Confidence Interval ON a P-value

text: I am using the coin package in R to do a permutation test and it gives output that looks like this:               > pvalue(pt)     [1] 0.1479615     99 percent confidence interval:     0.1450811 0.1508754       Why would one want a confidence interval for a p-value? How should it be interpreted, and how should it be used?


title: To rotate or not to rotate post-PCA and pre-cluster analysis

text: Questions in respect to rotation post-PCA have been answered before -> its all in the hands of the researcher... Same answer to the question if rotation (orthogonal or not) makes sense before plugging the components into a cluster analysis?


title: Reasons for autocorrelation in time-series residuals

text: Why are residuals usually autocorrelated in time-series data? Could it stem from the autocorrelation of the response variable? Is the reason that in some cases the differencing (i.e., the differences between adjacent values of the response variable) is used?


title: Help on subsetting data frames using multiple logical operators in R

text:
In R, I am trying to write a function to subset and exclude observations in a data frame based on three variables. My data looks something like this:               data.frame':   43 obs. of  8 variables:      $ V1: chr  "ENSG00000008438" "ENSG00000048462" "ENSG00000006075" "ENSG00000049130" ...
          $ V2: chr  "ENST00000008938" "ENST00000053243" "ENST00000225245" "ENST00000228280" ...      $ V3: chr  "ENSP00000008938" "ENSP00000053243" "ENSP00000225245" "ENSP00000228280" ...
          $ V4: chr  "PGLYRP1" "TNFRSF17" "CCL3" "KITLG" ...      $ V5: chr  "19" "16" "17" "12" ...
          $ V6: chr  "q13.32" "p13.13" "q12" "q21.32" ...      $ V7: int  46522415 12058964 34415602 88886566 8276226 150285143 29138654 76424442 136871919 6664568 ...
          $ V8: int  46526304 12061925 34417515 88974238 8291203 150294844 29190208 76448092 136875735 6670599 ...     >      What I am trying to do is this:               input = function(x) {             pruned1 <- subset(x[-which(x$V5 == 1 & x$V7 > 113818477 & x$V8 < 114658477), ])             pruned2 <- subset(pruned1[-which(pruned1$V5 == 1 & pruned1$V7 > 192461456 & pruned1$V8 < 192549912), ])             # and so on             }      So I want to take my original data frame, exlude observations that fulfil all three of the subsetting conditions, make a new object of the remaining observations (pruned1) and exclude any observations that fulfil a new set of condidtions and so on. The problem I am having seems to be the third condition in the first subsetting, as this returns:               'data.frame':   1 obs. of  8 variables:      $ V1: chr NA
          $ V2: chr NA      $ V3: chr NA
          $ V4: chr NA      $ V5: chr NA
          $ V6: chr NA      $ V7: int NA
          $ V8: int NA      While only using the first two conditions returns:               'data.frame':   38 obs. of  8 variables:      $ V1: chr  "ENSG00000008438" "ENSG00000048462" "ENSG00000006075" "ENSG00000049130" ...
          $ V2: chr  "ENST00000008938" "ENST00000053243" "ENST00000225245" "ENST00000228280" ...      $ V3: chr  "ENSP00000008938" "ENSP00000053243" "ENSP00000225245" "ENSP00000228280" ...
          $ V4: chr  "PGLYRP1" "TNFRSF17" "CCL3" "KITLG" ...      $ V5: chr  "19" "16" "17" "12" ...
          $ V6: chr  "q13.32" "p13.13" "q12" "q21.32" ...      $ V7: int  46522415 12058964 34415602 88886566 8276226 150285143 29138654 76424442 136871919 6664568 ...
          $ V8: int  46526304 12061925 34417515 88974238 8291203 150294844 29190208 76448092 136875735 6670599 ...     >      Doing this manually is not an options, since I have a large number of input files that I would like so prune in the same manner. Thanks in advance for any help!


title: monotonicity of power function

text: Following is Example 4.1.4 from Bickel and Doksum's mathematical statistics book. Given an iid sample $X_1, ..., X_n$. each $X_i$ has a Bernoulli distribution with unknown parameter $\theta$. Given a constant $\theta_0$, we want to test $H: \theta=\theta_0$ versus $K: \theta >\theta_0$. The test statistic is $S_n = \sum_{i=1}^n X_i$, and the testing procedure $I_{S_n \geq k}$. If $\theta_0 =0.3$ and $n=10$ and $k=6$, the power function of $\theta$ is a increasing function. See Figure 4.1.1. My question is to go beyond this example to more general case. In a general hypothesis testing task or some more general one than the example, is there some property of the power function that is similar to or generalizes monotonicity of the power function in the example? Thanks and regards!


title: Trendline formula in R

text: I have created a simple scatterplot to analyse the correlation between two variables, x and y. Using the abline command, I added a trendline to this data. I now want to know the formula for this trendline. Is there a command that will give this formula as one is able to do in a simple Excel scatterplot?


title: Specification of error terms in aov() in R

text: For a two-way repeated measures design, we can specify the model using `aov` in the following fashion:               aov(dv ~ iv1 + iv2 + Error(subject/(iv1*iv2)), data=dataset)      Take the sample data from `personality-project.org` (see code below), the error term for this particular model is `Error(Subject/(Task*Valence))`.               datafilename="http://personality-project.org/r/datasets/R.appendix4.data"     data.ex4=read.table(datafilename,header=T)   #read the data into a table     data.ex4                                      #show the data     aov.ex4=aov(Recall~(Task*Valence)+Error(Subject/(Task*Valence)),data.ex4 )      What if I specify the error term instead as `Error(Subject)`. **What kind of design is this?** From what I recall, whether you include a factor - e.g., whether `Task` and `Valence` are included - in the error term depends on whether a factor is a "within-subjects" factor. If, say, this particular recall study also included "sex" as a predictor, then we would not include it in the error term because it is a between-subjects factor. But I am not sure what it means when we don't include any of the predictors in the error term. **Is it basically an intercept-only model much like the one below specified using`lmer`**?               lmer(Recall~Task*Valence + (1|Subject), data=data.ex4)


title: Model Identification

text: Could somebody please explain why this model is "just identified" ![enter image description here](http://i.stack.imgur.com/yRCmk.png) As I see it, there are 5 * 4 / 2 = 10 variances/covariances, 4 observed means, giving 14 available degrees of freedom 5 DF are used on the 5 paths 2 intercepts are estimated 2 error variances are estimated 2 endogenous variances are estimated 2 exogenous variances are estimated Giving 13 used DF. What am I missing ?


title: PyMC for nonparametric clustering: Dirichlet process to estimate Gaussian mixture's parameters fails to cluster

text: ## Problem setup One of the first toy problems I wanted to apply PyMC to is nonparametric clustering: given some data, model it as a Gaussian mixture, and learn the number of clusters and each cluster's mean and covariance. Most of what I know about this method comes from video lectures by Michael Jordan and Yee Whye Teh, circa 2007 (before sparsity became the rage), and the last couple of days' reading Dr Fonnesbeck's and E. Chen's tutorials [fn1], [fn2]. But the problem is well-studied and has some reliable implementations [fn3]. In this toy problem, I generate ten draws from a one-dimensional Gaussian $\mathcal{N}(\mu=0, \sigma=1)$ and forty draws from $\mathcal{N}(\mu=4, \sigma=2)$. As you can see below, I didn't shuffle the draws, to make it easy to tell which samples came from which mixture component. ![Gaussian mixture model data](http://i.stack.imgur.com/JFwFZ.png) I model each data sample $y_i \sim \mathcal{N}(\mu_{z_i}, \sigma_{z_i})$, for $i=1,...,50$ and where $z_i$ indicates the cluster for this $i$th data point: $z_i \in [1,...,N_{DP}]$. $N_{DP}$ here is the length of the truncated Dirichlet process used: for me, $N_{DP}=50$. Expanding the Dirichlet process infrastructure, each $z_i$ cluster ID is a draw from a categorical random variable, whose probability mass function is given by the stick-breaking construct: $z_i \sim Categorical(p)$ with $p \sim Stick(\alpha)$ for a concentration parameter $\alpha$. Stick-breaking constructs the $N_{DP}$-long vector $p$, which must sum to 1, by first obtaining $N_{DP}$ i.i.d. Beta-distributed draws that depend on $\alpha$, see [fn1]. And since I'd like to the data to inform my ignorance of $\alpha$, I follow [fn1] and assume $\alpha \sim Uniform(0.3, 100)$. This specifies how each data sample's cluster ID is generated. Each of the $N_{DP}$ clusters has an associated mean and standard deviation, $\mu_{z_i}$ and $\sigma_{z_i}$. Then, $\mu_{z_i} \sim \mathcal{N}(\mu=0, \sigma=50)$ and $\sigma_{z_i} \sim Uniform(0, 100)$. (I was previously following [fn1] unthinkingly and placing a hyperprior on $\mu_{z_i}$, that is, $\mu_{z_i} \sim \mathcal{N}(\mu_0, \sigma_0)$ with $\mu_0$ itself a draw from a fixed-parameter normal distribution, and $\sigma_0$ from a uniform. But per http://stats.stackexchange.com/a/71932/31187, my data doesn't support this kind of hierarchical hyperprior.) In summary, my model is: $y_i \sim \mathcal{N}(\mu_{z_i}, \sigma_{z_i})$ where $i$ runs from 1 to 50 (the number of data samples). $z_i \sim Categorical(p)$ and can take on values between 0 and $N_{DP}-1=49$; $p \sim Stick(\alpha)$, a $N_{DP}$-long vector; and $\alpha \sim Uniform(0.3, 100)$, a scalar. (I now slightly regret making the number of data samples equal to the truncated length of the Dirichlet prior, but I hope it's clear.) $\mu_{z_i} \sim \mathcal{N}(\mu=0, \sigma=50)$ and $\sigma_{z_i} \sim Uniform(0, 100)$. There's $N_{DP}$ of these means and standard deviations (one for each of the $N_{DP}$ possible clusters.) Here's the graphical model: the names are variable names, see the code section below. ![Graph](http://i.stack.imgur.com/CtwMs.png) ## Problem statement Despite several tweaks and failed fixes, the parameters learned are not at all similar to the true values that generated the data. Currently, I'm initializing most of the random variables to fixed values. The mean and standard deviation variables are initialized to their expected values (i.e., 0 for the normal ones, the middle of their support for the uniform ones). I initialize all $z_i$ cluster IDs to 0. And I initialize the concentration parameter $\alpha=5$. With such initializations, 100'000 MCMC iterations simply can't find a second cluster. The first element of $p$ is close to 1, and nearly all draws of $\mu_{z_i}$ for all data samples $i$ are the same, around 3.5. I show every 100th draw here for the first twenty data samples, i.e., $\mu_{z_i}$ for $i=1,...,20$: ![Means with zero-initialized cluster ID](http://i.stack.imgur.com/JqXCg.png) Recalling that the first ten data samples were from one mode and the rest were from the other, the above result clearly fails to capture that. If I allow random initialization of the cluster IDs, then I obtain more than one cluster but the cluster means all wander around the same 3.5 level: ![Means with random-initialized cluster ID](http://i.stack.imgur.com/VPU6w.png) This suggests to me that it's the usual problem with MCMC, that it can't reach another mode of the posterior from the one it's at: recall that these different results happen after just changing the _initialization_ of the cluster IDs $z_i$, not their priors or anything else. Am I making any modeling mistakes? Similar question: http://stackoverflow.com/q/19114790/500207 wants to use a Dirichlet distribution and fit a 3-element Gaussian mixture and is running into somewhat similar problems. Should I consider setting up a fully conjugate model and using Gibbs sampling for this kind of clustering? (I implemented a Gibbs sampler for the parametric Dirichlet distribution case, except using a fixed concentration $\alpha$, back in the day and it worked well, so expect PyMC to be able to solve at least that problem handily.) ## Appendix: code               import pymc     import numpy as np          ### Data generation          # Means and standard deviations of the Gaussian mixture model. The inference     # engine doesn't know these.     means = [0, 4.0]     stdevs = [1, 2.0]          # Rather than randomizing between the mixands, just specify how many     # to draw from each. This makes it really easy to know which draws     # came from which mixands (the first N1 from the first, the rest from     # the secon). The inference engine doesn't know about N1 and N2, only Ndata     N1 = 10     N2 = 40     Ndata = N1+N2          # Seed both the data generator RNG  as well as the global seed (for PyMC)     RNGseed = 123     np.random.seed(RNGseed)          def generate_data(draws_per_mixand):         """Draw samples from a two-element Gaussian mixture reproducibly.              Input sequence indicates the number of draws from each mixand. Resulting         draws are concantenated together.              """         RNG = np.random.RandomState(RNGseed)         values = np.hstack([RNG.normal(means[i], stdevs[i], ndraws)                             for (i,ndraws) in enumerate(draws_per_mixand)])         return values          observed_data = generate_data([N1, N2])               ### PyMC model setup, step 1: the Dirichlet process and stick-breaking          # Truncation level of the Dirichlet process     Ndp = 50          # "alpha", or the concentration of the stick-breaking construction. There exists     # some interplay between choice of Ndp and concentration: a high concentration     # value implies many clusters, in turn implying low values for the leading     # elements of the probability mass function built by stick-breaking. Since we     # enforce the resulting PMF to sum to one, the probability of the last cluster     # might be then be set artificially high. This may interfere with the Dirichlet     # process' clustering ability.     #     # An example: if Ndp===4, and concentration high enough, stick-breaking might     # yield p===[.1, .1, .1, .7], which isn't desireable. You want to initialize     # concentration so that the last element of the PMF is less than or not much     # more than the a few of the previous ones. So you'd want to initialize at a     # smaller concentration to get something more like, say, p===[.35, .3, .25, .1].     #     # A thought: maybe we can avoid this interdependency by, rather than setting the     # final value of the PMF vector, scale the entire PMF vector to sum to 1? FIXME,     # TODO.     concinit = 5.0     conclo = 0.3     conchi = 100.0     concentration = pymc.Uniform('concentration', lower=conclo, upper=conchi,                                  value=concinit)          # The stick-breaking construction: requires Ndp beta draws dependent on the     # concentration, before the probability mass function is actually constructed.     betas = pymc.Beta('betas', alpha=1, beta=concentration, size=Ndp)          @pymc.deterministic     def pmf(betas=betas):         "Construct a probability mass function for the truncated Dirichlet process"         # prod = lambda x: np.exp(np.sum(np.log(x))) # Slow but more accurate(?)         prod = np.prod         value = map(lambda (i,u): u * prod(1.0 - betas[:i]), enumerate(betas))         value[-1] = 1.0 - sum(value[:-1]) # force value to sum to 1         return value          # The cluster assignments: each data point's estimated cluster ID.     # Remove idinit to allow clusterid to be randomly initialized:     idinit = np.zeros(Ndata, dtype=np.int64)     clusterid = pymc.Categorical('clusterid', p=pmf, size=Ndata, value=idinit)          ### PyMC model setup, step 2: clusters' means and stdevs          # An individual data sample is drawn from a Gaussian, whose mean and stdev is     # what we're seeking.          # Hyperprior on clusters' means     mu0_mean = 0.0     mu0_std = 50.0     mu0_prec = 1.0/mu0_std**2     mu0_init = np.zeros(Ndp)     clustermean = pymc.Normal('clustermean', mu=mu0_mean, tau=mu0_prec,                               size=Ndp, value=mu0_init)          # The cluster's stdev     clustersig_lo = 0.0     clustersig_hi = 100.0     clustersig_init = 50*np.ones(Ndp) # Again, don't really care?     clustersig = pymc.Uniform('clustersig', lower=clustersig_lo,                               upper=clustersig_hi, size=Ndp, value=clustersig_init)     clusterprec = clustersig ** -2          ### PyMC model setup, step 3: data          # So now we have means and stdevs for each of the Ndp clusters. We also have a     # probability mass function over all clusters, and a cluster ID indicating which     # cluster a particular data sample belongs to.          @pymc.deterministic     def data_cluster_mean(clusterid=clusterid, clustermean=clustermean):         "Converts Ndata cluster IDs and Ndp cluster means to Ndata means."         return clustermean[clusterid]          @pymc.deterministic     def data_cluster_prec(clusterid=clusterid, clusterprec=clusterprec):         "Converts Ndata cluster IDs and Ndp cluster precs to Ndata precs."         return clusterprec[clusterid]          data = pymc.Normal('data', mu=data_cluster_mean, tau=data_cluster_prec,                        observed=True, value=observed_data)      ## References   1. fn1: http://nbviewer.ipython.org/urls/raw.github.com/fonnesbeck/Bios366/master/notebooks/Section5_2-Dirichlet-Processes.ipynb   2. fn2: http://blog.echen.me/2012/03/20/infinite-mixture-models-with-nonparametric-bayes-and-the-dirichlet-process/   3. fn3: http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm.html#example-mixture-plot-gmm-py


title: Fitting curved lines to scatterplots in R

text: I have created a scatterplot. I want to see the trend of how the number of fishing cat scats vary with an increase in perimeter of water body. Hence my response variable is `number of fishing cat scats` and predictor variable is `perimeter of water body`. After creating a scatterplot, I gave an `abline` (`lm`) command, which has given me a linear regression line, which doesn't exactly portray the relationship between `number of fishing cat scats` and `perimeter of water body`. What command do I need to use to generate a curve which fits the data? I am using R as my statistical software.


title: Time-scale in cox porportinal hazards model with multiple time-dependent variables

text: I'm doing survival analysis and I have three time-dependent variables, with different time scales each. I've already done univariate analysis and I want to proceed with multivariate analysis. My questions is how can I choose or deal with different time-scales? I'm using R, survival package.


title: Probability of exceedance and reliability of a sample range estimation

text: Consider that $P$ is the water pressure coming out of a valve $A$. Let $P_{dif}$ be the difference between the maximum and the minimum pressure of valve $A$: $$P_{dif}≔P_{max}-P_{min}$$ Now, what I want to do is estimate $P_{dif}$. In order to do that, I take a number of water pressure samples from valve $A$. Let $S$ be a set of $3$ measured samples: $$S = {5,7,1}$$ That is, $S$ contains 3 random samples of $P$, therefore, by placing S in ascending order, I can estimate $P_{dif}$ like that: $$\hat{P}_{dif} = S_{(n:3)}-S_{(n:1)} = 7-1 = 6$$ Questions:   * What is the probability of exceeding this estimation ($\hat{P}_{dif}$). That is, what is the probability that the population parameter ($P_{dif}$) will exceed the estimation ($\hat{P}_{dif}$).   * Reliability of the estimation ($\hat{P}_{dif}$): this refers to the probability that the estimation is wrong (it may be possible to find this by deriving the confidence interval of the sample range, but I don't know how). Any help towards this will be greatly appreciated.


title: Book for non parametric statistics

text: What would be a good book for non-parametric statistics. Not just the introduction but advanced level. Also I am looking at something I can use for learn and not for reference. In particular I am looking for a book that can contain basics behind non- parameteric methods, non-parameter inference, methods to evaluate non parametrics, e.g., KS test, $t$ test, etc. , bootstrapping ....


title: Matrix decomposition

text: I have a symmetric matrix $V$ of order 10. I want to decompose $V$ in such a way that $$V=SS'$$ with $S$ being non triangular. The matrix $S$ has some restrictions that 45 cells have ‘0’ values and 55 other remaining cells are supposed to be solved. I am looking for solving this using sas. Here is the structure of $S$ matrix (all 1's are to be replaced by variables; 55 variables and 45 zeros) $$ \begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & \\ 1 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & \\ 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & \\ 1 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & \\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \\ \end{bmatrix}$$


title: Is the matrix dimension important for performing a valid PCA?

text: If $X$ is a $m × n$ matrix, where $m$ is the number of measurement types (variables) and $n$ is the number of samples, would it be correct to perform a PCA on a matrix that has $m \geq n$ ? If not, please provide some arguments why would this be a problem. I remember having heard that doing such an analysis would be invalid, but the Wikipedia page for PCA doesn't mention a low $n/m$ ratio as being a potential limitation for using the method. Please note that I am a biologist and aim at a more practical answer (if possible).


title: Chi-square plot to check for multivariate normality: Should I use the sample covariance or "normal" covariance?

text: I'd like to check if my observations are normally distributed by using a chi- square plot. To calculate the generalized squared distances (squared Mahalanobis distances) I need the covariance matrix of my observations. Is it correct that I have to use the sample covariance instead of the "normal" covariance $$E[ (X - \bar{X})(Y - \bar{Y}) ]$$ in all cases where I don't know the actual E[X] / E[X] and have to estimate it with the respective mean?


title: Ornstein-Uhlenbeck vs. Random Normal?

text: I was asked to implement an Ornstein-Uhlenbeck process in one of my simulations. I have coded the process to visualize the results and I was wondering, if my first value is at the mean, why bother using an O-U process? I thought the advantage of an O-U process is that it is mean-reverting, but if I am starting at the mean and just want to fluctuate about the mean, would a random Normal generator work equally well? Do I gain anything with the O-U model?


title: Econometrics: Sargan test

text: Here are 3 questions about econometrics and `R` codes. Test the endogeneity of the variable EDUC:               etudes1 <- lm(EDUC ~ EXPER+EXPERSQ+SMSA+SOUTH+FATHEDUC+MOTHEDUC, data=Dataset)          Call:     lm(formula = LWAGE ~ EDUC + EXPER + EXPERSQ + SMSA + SOUTH + residuals(etudes1))          Residuals:          Min       1Q   Median       3Q      Max      -1.73299 -0.21985  0.02344  0.24196  1.35879           Coefficients:                          Estimate Std. Error t value Pr(>|t|)         (Intercept)         3.9002276  0.1998339  19.517  < 2e-16 ***        EDUC                0.1214847  0.0116061  10.467  < 2e-16 ***       EXPER               0.1073102  0.0092801  11.563  < 2e-16 ***        EXPERSQ            -0.0024888  0.0004016  -6.197 6.86e-10 ***        SMSA                0.1355027  0.0193295   7.010 3.15e-12 ***        SOUTH              -0.1311756  0.0178376  -7.354 2.70e-13 ***        residuals(etudes1) -0.0463500  0.0123884  -3.741 0.000188 ***     ---          Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.3811 on 2213 degrees of freedom     Multiple R-squared: 0.251,  Adjusted R-squared: 0.2489      F-statistic: 123.6 on 6 and 2213 DF,  p-value: < 2.2e-16         1. The "residuals" variable is significant (p-value = 0.000188), which means that we reject $H_{0}$: the variable EDUC is **not** EXOGENOUS. That means that we must use the 2SLS (instead of the OLS) to estimate the efficiency of the education (EDUC). Is it right? **The Sargan Test**               Call:     lm(formula = residuals(iv1) ~ EXPER + EXPERSQ + SMSA + SOUTH + FATHEDUC + MOTHEDUC)          Residuals:          Min       1Q   Median       3Q      Max      -1.90114 -0.23667  0.02866  0.24715  1.43550           Coefficients:                    Estimate Std. Error t value Pr(>|t|)       (Intercept) -7.948e-03  5.472e-02  -0.145    0.885         EXPER       -1.996e-04  8.019e-03  -0.025    0.980        EXPERSQ      4.986e-06  4.113e-04   0.012    0.990        SMSA         1.486e-03  1.907e-02   0.078    0.938         SOUTH       -1.191e-04  1.771e-02  -0.007    0.995         FATHEDUC    -3.054e-03  2.998e-03  -1.019    0.309         MOTHEDUC     3.669e-03  3.548e-03   1.034    0.301          Residual standard error: 0.3907 on 2213 degrees of freedom     Multiple R-squared: 0.000607,   Adjusted R-squared: -0.002103      F-statistic: 0.224 on 6 and 2213 DF,  p-value: 0.9691           jour_sargan <- summary(sargan)       r2iv   <- jour_sargan$r.square     # $     r2iv                               # [1] 0.0006070252         obs    <- jour_sargan$r.square*500         obs                                # [1] 0.3035126        pvalue <- 1-pchisq(obs, df=1)         pvalue                             # [1] 0.5816887        1. Do I have to look at the 2 p-values: FATHEDUC (p-value = 0.309) and MOTHEDUC (p-value = 0.301)? Or do I have to look at the final p-value of the Sargan Test (p=0.5816887) to take my final decision?   2. We don't reject $H_{0}$, so the instrumental variables (FATHEDUC and MOTHEDUC) are **not** EXOGENOUS. Is it right?


title: Having trouble understanding cross-validation results from scikit-learn

text: Actually, my question may just be about cross-validation in general. Here's what I'm doing: I'm trying to come up with a model using scikit-learn to learn on some data I've got. I've decided to use an SVM, using various kernels, to do the modelling. I've got about 50,000 data points from which to extract features. In an effort to make sure that my model is not over- or under- fitting, I've decided to run all of my models through cross-validation using scikit-learn's `cross_validation` functionality. I'm setting aside 40% of my training data for cross-validation, and so training on 60%. I do this iteratively until I come up with a set of features and a model that gives me a cross-validation score of about 0.96. Great! Here's the problem - when I use this model to predict results for my test data, I only get a score of about 0.79! I don't understand that result. My question is, am I misunderstanding the cross validation score? Shouldn't I be able to expect similar results for my test data when using the model cross-validated to 0.96? I even used the GridSearchCV to come up with the best parameters to use for the SVM kernel. I also made sure to train on the full set of training data when training my model before running predict. This is my first real attempt to use machine learning for a cool project, and I'm totally confused on my expectations. Any guidance would be greatly appreciated! Thanks, guys and gals!


title: How to calculate overall error from a series of dependent experiments?

text: Salmon were collected and induced to spawn. The percent of the fish collected who spawned successfully was 63%, so P(S)=0.63. A sample of the eggs were then fertilized and the percent of eggs that were successfully fertilized was 70%, so P(F|S)=0.70. The fertilized eggs were allowed to hatch and the percent of fertilized eggs that produced normal larvae was 80%, so P(L|FS)=0.80. The probability that the entire process is completed successfully is P(L)=(0.63)(0.70)(0.80)=0.3528. The sample standard deviation for fertilization is reported as +- 28, and sample standard deviation for larvae is reported as +- 22 (both seem large). The standard deviation for spawning data is not reported. Is there any way to estimate the error in the probability that the entire process is completed successfully? The data at each stage is distributed binomially, so what is the distribution of the product? Is it possible to calculate the standard deviation of the final product using only this information?


title: Distribution of p-values from multiple experiments

text: My experimental setup is the following:   1. compute an optimal solution $S$ to problem $A$ on instance $I$,   2. evaluate test statistic $f(S)$,   3. compute $f(S_i)$ for $S_i$ $(0 \leq i \leq n)$ drawn from a uniform distribution of feasible solutions to $A$, and   4. compute a p-value for S. The null hypothesis is that $f$ is not correlated with the optimality of $A$. The alternative is that $f$ is. The p-value for $S$ is then $l/n$ where $l$ is the number of $S_i$s such that $f(S_i) > f(S)$. If we repeat this test for many instances $I_j$, we get a distribution of p-values. I would like to determine if we can reject the null hypothesis over all tests. I've looked at Fisher's method, but I have some pvalues that are so small, that the result is 0. If I make a histogram of the pvalues I get the following: ![enter image description here](http://i.stack.imgur.com/jwuuo.png) There is a clear bias towards very small p-values, but what worries me is the uniform looking distribution of the other p-values. Does our data reject the null hypothesis? EDIT: What if the left-most bar indicated less tests with such small pvalues. How would we evaluate the following distribution of p-values? ![enter image description here](http://i.stack.imgur.com/klf0W.png)


title: Missing values - for cases of "I do not know" or "Does not apply"

text: If I have a questionnaire with items on a Likert scale (1 to 5) and a sixth option "I do not know/Does not apply", is it okay to consider a "I do not know/Does not apply" response as a missing value and use an imputation method, just like I would consider and impute the missing values in the case when respondent to not answer at all for certain items?


title: Probability on Z-score

text: I have a data of 30 values that represent returns on investment and which have a mean of 7 and standard deviation of 14. We are talking here about normal distribution. Now I am being required to do a following: > If you are given a choice between the top 6% of returns and a return of 25%, > which option would you choose? Between what two values of returns > (symmetrically distributed around the mean) will 32.6% of all possible > returns contained? I was thinking to sum top 6% of the values of a given data, and divide it by number of it, and to find out z score of 25, and to compare those two results. Is my approach good? I am a bit confused about how to find between what two values of returns data I will get 32.6%?


title: Help: Am I using the correct test to detect differences from zero (multivariate ANOVA repeated mesures)

text: 01 August 2014 11:17 Hi, I have 5 measures repeated 3 times under 3 different scenarios (conditions). The measures are survey questions and participants use a visual analogue scale to mark a response from 0.0 to 15.0. The same 5 measures (questions) are repeated for each of the 3 scenarios (conditions) making 15 measures in total. I have one group: The Participants. The control group is a 'virtual control group' where every measure (15) is set to '0.00'. The survey questions are randomised and asked at the same time. I would like to understand: Differences from 0.00: a). If there is a significant difference from 0.00 in any of the 5 measures (in all conditions). b). If there is a sig diff from 0.00 in any one of the 3 conditions. c). Understand any interactions between conditions and measures Differences from each other: c). Understand if there is a sig difference of the 5 measures between themselves across all 3 conditions. d). Understand if there is a sig difference of the 5 measures in any of 3 conditions. I'm using GLM in SPSS (repeated measures) in have found the following ways: Conditions(3 levels), Measures (5 levels) as factors with repeated contrasts or, Conditions (3 levels) as a factors and Measures (5) as measures with repeated contrasts I am doing this with the participant dataset (n=30) and a dataset that includes the participants and the control group (n=60) where I have added control participants with their measures set to '0.00' and differentiated the 2 groups by adding a group variable set at 1 or 0. My questions revolve around homogeneity of variance: I'm unsure that I can run a repeated measures MANOVA as the variance in the control group is zero whilst I have non-normal participant data with large variance. Any thoughts - on how I should be doing this? Thank you


title: Correlation of errors between regression equations?

text: I now that one can use Seemingly Unrelated Regressions when for an observation i, there are two equations whose errors are connected. What is the intuition here among errors of difference equations have some sort of relationship? In the SUR case, is it simply because unobservable factors (say, for example, general attitude towards spending money) that explain things like consumption for food in observation 1 will also explain consumption for clothing, for the same observation 1? Also, how does one call correlation between errors in this case? I understand it's not serial correlation, nor spatial correlation in this case.


title: How does pooling and resampling affect variance of sample mean?

text: Suppose I have $N$ independent random variables $X_n$. I draw a sample of predetermined size $K_n$ from each of them. Denote the average of each sample $\bar{\hat{X}}_n$, and the total number of observations $M = \sum_{n=1}^{N} K_n$. Next, I use all the $M$ sample observations $\hat{X}_{n,k}$, with no regard as to which r.v. they came from, to construct the empirical distribution $\hat{F}$. Next, I create $N$ new i.i.d. random variables $Y_n$ distributed as $\hat{F}$, and draw from the them i.i.d. samples of the same sizes $K_n$. Denote the average of each sample $\bar{\hat{Y}}_n$. Let $T$ be some commonly used sample statistic that represents the dispersion; e.g., variance (or interquartile distance). Finally, I compare $T_X = E[T(\bar{\hat{X}}_1, …, \bar{\hat{X}}_N)]$ and $T_Y = E[T(\bar{\hat{Y}}_1, …, \bar{\hat{Y}}_N)]$. Under what (preferably intuitive) conditions would I find $T_X = T_Y$, $T_X > T_Y$, or $T_X < T_Y$? Motivation: I've seen bootstrap used in social sciences to analyze a large group of independent random variables. Specifically, the researcher would examine the variation in the sample statistics across these variables. If it seems that the variation is bigger than "expected", the researcher would draw some conclusions about the original r.v. But I could never figure out precisely what null hypothesis this approach tests, so I'm also unsure about what conclusions would be justified.


title: How to compare measurement data resulting from automatic processing vs. visual inspection?

text: First set: 12, 23, 24, 14, 29, 30, 31, 16. Second set: 12, 22, 24, 14, 30, 30, 31, 16. These data are obtained using a software (first set) and visual observation (second set) of the same process. I do not know if it matters it is expressed in minutes per 2 hours interval. Thanks.


title: Analysing periodicity in sparse time series

text: I am interested in looking for periodicities in a several day long recording of electrical activity. The traces present a very steady baseline over which, from time to time, some short events (300-500 ms) appear (hence the _sparse_ in the title, although I am not sure it is the right term to use). Now, the nature of the data prevents me to analyse the trace as a whole (~20 days recorded at 10KHz, my computer will definitely not handle that), so I just wrote a series of routines to find the interesting events, and export their date/time (also their shape, but that is another story) to a text file. So, I end up with a series of times to analyse. Just by looking at the raw trace what I can see is:   1. A ~2 hours periodicity in the events   2. The events appear in groups. That is, I have a series of 20-30 events in the course of a couple of minutes, then nothing for ~2 hours, then other 20-30 events and so on. Note that the number of events in a group is variable, can be 30 in one and 5 in the next one.   3. The periodicity is fairly obvious, but it's not perfect: events happen roughly at 2 hours interval, and the exact time may vary from day to day.   4. There may be superimposed periodicities. In particular, I see events at a specific time of the day, which may or may not fit in the abovementioned 2 hours periodicity. What can I do to statistically determine the periodicity? What I have tried to do is to make a frequency histogram of the events, binning every 10 minutes, which empirically seemed like a sensible bin size, and then looking at the ACF or the FFT of the counts, but only in a few cases something pops out with the ACF (nothing with FFT). So, how would you analyse this type of problem? **Bonus question** : in certain cases I have missing data (1 or 2 days missing for technical reasons). How would I account for that? PS: I am using R, but any non-R solution will do as well! **EDIT** Here is some sample data to play with: http://dl.dropbox.com/u/11676289/exampletimes.txt Here is a plot of the times: ![Plot of the times of the events](http://i.stack.imgur.com/Oqy9z.jpg) and their histogram, in 10 minutes breaks ![histogram](http://i.stack.imgur.com/KvYAh.jpg)


title: Plot by cdplot in R doesn't fit with data

text: I am trying to understand how cdplot in R behaves, but I am missing something. When I copy/paste the following example given in the documentation for cdplot:               > fail <- factor(c(2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1,                        1, 2, 1, 1, 1, 1, 1),                levels = 1:2, labels = c("no", "yes"))          > temperature <- c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,                        70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81)     > cdplot(fail ~ temperature)      this is the plot that I get: ![Conditional Density Plot for fail given temperature](http://s18.postimage.org/88qqnrjqh/cdplot.jpg) Based on this graph, my expectation is that $P(\text{fail} = 1 | \text{temperature} = 67) \approx 0.77$ But when I look at the data:               > which(temperature == 67)     6 7 8          > which(temperature == 67 & fail == 1)     integer(0)      which together imply that $P(\text{fail} = 1 | \text{temperature} = 67) = 0$ (at least based on this data). What am I missing here? Why is there such a huge difference between what cdplot gives and a basic hand calculation? Also, for some variables, I get plots of this shape: cdplot gone wild What has happened in such a plot?


title: Help: Repeated Measures Multivariate Analysis in SPSS

text: I have 5 measures repeated 3 times under 3 different scenarios (conditions). The measures are survey questions and participants use a visual analogue scale to mark a response from 0.0 to 15.0. Thus the same 5 measures (questions) are repeated for each of the 3 scenarios (conditions) making 15 measures in total. I have one group: The Participants. The control group is a 'virtual control group' where every measure (15) is set to '0.00. The survey questions are randomised and asked at the same time. I would like to: Differences from 0.00 a). Understand if there is a significant difference from 0.00 in any of the 5 measures (in all conditions). b). Understand is there is a sig diff from 0.00 in any one of the 3 conditions. c). Understand any interactions. Differences from each other: c). Understand if there is a sig difference of the 5 measures between themselves across all 3 conditions. d). Understand if there is a sig difference of the 5 measures in any of 3 conditions. I'm using GLM in SPSS (repeated measures) in have found the following ways: Conditions(3 levels), Measures (5 levels) as factors with repeated contrasts or, Conditions (3 levels) as a factors and Measures (5) as measures with repeated contrasts I am doing this with the participant dataset (n=30) and a dataset that includes the participants (n=30) and a


title: Using multinomial logit to rank observations

text: I'm trying to devise a model to rank a number of links according to their popularity. The links refer to upcoming events and job offers, and ideally I'd like to have a reasonably simple model that, given the values of some covariates such as the number of clicks on the link and the number of days remaining till the event, would provide me with a value I could use to rank the observations according to their popularity. What is actually needed is a list of the top-ten events. My initial, naive, approach was to ask the committent to rank a small subset of the observations into three ordinal categories, namely "High", "Average" and "Low" popularity. I then reasoned that having something like a training dataset I could try fitting a multinomial logistic regression, and then I could sort the observations according to the predicted probability of falling into the "High" popularity category; i.e., given a subset of covariates $X$ the closer the predicted probability $P(``High"|X)$ gets to 1 the closer to the top that observation gets. What do you think about this approach? This seems a problem suited for machine learning methods, about which sadly I know next to nothing. Is there a more exact term for what I'm trying to do? I'd be grateful for any pointer to the relevant literature, altough keep in mind the "solution" has to be quite simple, a single equation would be ideal. Thank you for your time!


title: Tendency of Early stop boundary to maximum in a neural network optimization structure

text: I have an optimization algorithm that optimize numbers of neurons, early stopping (stopping after X) and best fractures of an binary classification problem with (`patternnet` and `feedforwardnet` with threshold of 0.5). Besides it I use 5-fold cross validation in every iteration of optimization algorithm. My optimization algorithm (for instance genetic algorithm) select the number of layers and neurons in every layer (maximum number of layers is 2). The best outputs of optimization algorithm based on best classification accuracy results maximum early stopping that I had set for early stopping limitations. Why? (For instance when I set early stopping in range [6 20] searching boundary I have 20 as best output for early stopping and when I set [6 30] I have 30.) Totally what is your idea about this structure? Do you think it is good to optimize Early stop criteria or I must change this optimization structure. You can refer to this link about another problem that I had in my system : Performance comparison of “patternnet” and “newff” for binary classification in MATLAB R2014a Thanks.