title: Compare several regression lines ANCOVA?

text: I´m about to make a comparison I had never performed so I ask you to excuse me if the question is too obvious. I´m trying to analyze the relation between number of movements and prey capture time on spiders under different temperature treatments. I have four regresion lines (corresponding to different temperatures) and I want to compare the slopes of these lines. I performed a preliminar covariance analysis in R, and it indicates there is no correlation between the data I´m analyzing. Nevertheless, reading some posts. I´ve found out that I made some violations, because data are very non- normal and the variances are non-homogeneous. After trying several transformations I could not get the conditions. So my questions are: 1) Do I need to know another assumption of the test? 2) If correct, does it exist a non parametric alternative to the ANCOVA or another way to compare the slopes? 3) Finally since I´m getting familiar with the test, I´d like to know why in some cases is it necessary to perform an ANOVA after making the ANCOVA? Just in case, I attach the preliminary analysis I performed in R. Thanks for any help or advice.               ancova <- lm(Tiempo~Movimientos*Grupo)     summary(ancova)               Coefficients:                              Estimate Std. Error t value Pr(>|t|)     (Intercept)          13.4444     1.0521  12.778 1.94e-13 ***     Movimientos          -1.0000     0.3051  -3.277  0.00272 **     GrupoG2              -5.8089     3.6488  -1.592  0.12223     GrupoG3              -3.7914     3.5522  -1.067  0.29463     GrupoG4               0.7345     1.4494   0.507  0.61615     Movimientos:GrupoG2   0.6356     0.4424   1.437  0.16148     Movimientos:GrupoG3   0.4082     0.5616   0.727  0.47321          Movimientos:GrupoG4   0.4184     0.3879   1.079  0.28959          anova(ancova)          Analysis of Variance Table          Response: Tiempo                            Df Sum Sq Mean Sq  F value    Pr(>F)     Movimientos        1 429.04  429.04 177.2294 6.988e-14 ***     Grupo              3  61.10   20.37   8.4125 0.0003554 ***     Movimientos:Grupo  3   5.33    1.78   0.7341 0.5401739     Residuals         29  70.20    2.42     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Coefficient marginal to interactions in linear regression

text: Consider this model: $$y = \beta_0 +\beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + \varepsilon$$ Somebody told me today that the coefficient for the main effect of $x_1$ (i.e., $\beta_1$) will be 'marginal to the interaction'. Could somebody explain what 'marginal to the interaction' means?


title: Which measure to use to classify some distributions

text: I have a lot of data sets corresponding to particular events, each event is a collection of some $N$ real values between 0 and 1. Some of these events don't interest me, and I want to filter them out. To this end I want to associate with each event a value that would let me order events from the most useful to least useful and institute an arbitrary cutoff. **Let's consider some histograms:** This is an example of a good event, that I definitely want to keep. ![The good](http://i.stack.imgur.com/fnSrC.png) This is an example of a bad event, that I want to potentially filter out. ![The bad](http://i.stack.imgur.com/PNUkt.png) This is an example of an obviously bad event, that I definitely want to filter out. ![The ugly](http://i.stack.imgur.com/7Zl38.png) What kind of measure should I use to classify these events? I've tried stabbing in the dark with kurtosis (gives best results), skewness, means and central moments.


title: Random forest computing time in R

text: I am using the party package in R with 10,000 rows and 34 features, and some factor features have more than 300 levels. The computing time is too long. (It has taken 3 hours so far and it hasn't finished yet.) I want to know what elements have a big effect on the computing time of a random forest. Is it having factors with too many levels? Are there any optimized methods to improve the RF computing time?


title: GBM relative influence

text: For a regression problem, I used gradient boosting machines and assessed RMSE. My dataset is comprised of 34 features and 10,000 records. Only 2 predictors were considered 'important' (importance for other predictors happened to be zero): they are both factor features with 300 levels or more, so lots of predictions on new data set gets the same result. When I delete these two features, I get nearly the same RMSE score even if only 3 or 4 predictors highlight some relative influence (again, with zero influence for other predictors). What could explain this result? Should I be concerned with levels of factor predictors?


title: Formula for p-value of Tukey HSD

text: What is the formula for calculating p-value of Tukey HSD http://onlinestatbook.com/calculators/tukeycdf.html


title: How can I calculate a rate between two "moving" data sets?

text: I own a software business that follows a common pattern: Interested people can download the software and try it free for 30 days, and if they like it, they can buy it at any time during those 30 days (or after). I'd like to calculate my "conversion rate," which I will define as the percentage of downloads that lead to sales within 30 days. I don't track individual users, so I have no way of knowing whether downloader X ever ended up buying. I do, however, have time series data for the number of downloads each day and the number of sales each day. I have this data going back several years. I would like to be able to do any of a few things with this information: 1. Compare this rate to a similar time period last year. 2. Determine whether the rate is generally increasing or decreasing. 3. Predict what the next X days sales might look like based on the last X days downloads. Given that information and those goals, is it possible to determine (or at least approximate), a useful conversion rate? If so, how is it done? (I'm a programmer who had one intro stats class in college, which was mostly about the probability of rolling dice. Answers that explain it like I'm 5 would be greatly appreciated)


title: How exactly does Chi-square feature selection work?

text: I know that for each feature-class pair, the value of the chi-square statistic is computed and compared against a threshold. I am a little confused though. If there are $m$ features and $k$ classes, how does one build the contingency table? How does one decide which features to keep and which ones to remove? Any clarification will be much appreciated. Thanks in advance


title: Is 2-sigma limit applicable as a measure of variance to a distribution?

text: There is a paper that uses something like the 2-sigma limit to check if a quantity is more-or-less constant. More specifically in this paper on Text Detection (Sec. 3.2 3rd para last line) this limit is used to check for constant stroke width i.e. something like if mean > 2*sigma. This relation is not explicitly stated in the paper (it writes about a paragraph to the same effect), so I tested it on text and works reasonably well in discriminating between text/non-text on the basis of constant stroke-width. The question is - why is the 2*sigma limit applicable as a measure of a quantity's constant nature? And in case this limit is satisfied, what conclusions can be drawn? Also interested in mathematically proving what is obtained from experimental results. Basically I am searching for a proper background on the 2-sigma limit: this one isn't so helpful.


title: Is there a way to model the error term in a linear regression with R?

text: I'm tryin to estimate a pretty basic regression. I have a dataset containing `x` and `y` and would like to esimate the following model in R               y = b0 + b1*x1 + b2*epsilon      where epsilon=error term. I'm trying the lm() function in R, but do not know how to include the error term in the regression. Does anyone know how to include the error term. Is there any good packages for R to do this? Thanks!


title: What is the best linear regression model to use when the shape of the data resembles a cubic distribution?

text: I have a data set with a distribution of one variable against the other resembling a cubic one (rises to some point and then falls to a steady level without a consequent rise). I know in which cases to use log-linear, log-lin, lin-log, and reciprocal or log reciprocal linear models, but I am not sure what to do here (I have checked all of the above and they not surprisingly turned out to be a bad fit). Is there any linear model that would help me in this case?


title: standardization to obtain unit variance

text: I've come across some papers in where certain forecast errors are standardized to have unit variance. Unfortunately that's the only information they provide and I have no idea on how to obtain/calculate their results. Assuming I have a vector of 3 year forecasts that i would like to standardized to have unit variance. I thus calculate the forecast errors such as actual value-forecasted value and get a vector of forecast errors. I've came up with the following: $x$ equals my forecast errors of the 3 year forecasts $x.stand=(x-mean(x))/std(x)$ is this correct or would this only be valid for one year forecasts?


title: Estimator (preferably unbiased) of $\ln (\text{E}[X])$

text: Given the distribution function of random variable $X$ I know how to estimate its mean. What would be an estimator (preferably unbiased) of $\ln(\text{E}[X])$ ?


title: Annual Probability of Occurrence

text: I am trying to determine, based on historical record, what the annual probability of occurrence would be for a specific hazard. Take tornadoes. Over a 59 year time span I have documented 45 tornadoes. Some years, across that 59 year span, saw no events occur which leaves me with 25 years out of 59 that had tornado occurrences. So how would I determine what my annual probability of an occurrence would be (ie: 10% chance of tornado a year).... I am not a mathematician, so please explain in a way that I can apply the formula to other historical events. I appreciate any help you can give.


title: Connections between groups in a hive plot

text: In a Hive Plot nodes are placed on lines and connections between groups are highlighted. This plot seems to lose the connections between nodes on the same line. It seems like an intentional loss of data to build a better plot, but is it possible to not lose this information in such a plot?


title: I study subpopulation so huge that chi-squares test insists it is not representative. How do I get confidence intervals for my population?

text: I want to get results about population A. But all I can study is subpopulation C of B, where B is a subpopulation of A (that is, C $\subset$ B $\subset$ A). Each member of A is classified into one of the groups G1, G2, G3 and G4. Selection from A into B is not random, but selection from B into C is random, although some members of B have better chances to be selected into C. To test if subpopulation C is representative of A, I constructed the "contingency" table ...                     A             B          C        Total?     G1     9817678      621718     90196    10,529,592     G2     1280567      191298     30799    1,502,664     G3     426856       79707      11000    517,563     G4     9817678      701425     87996    10,607,099     Total  21,342,779   1,594,148  219,991  23,156,918      ... ignoring my doubts about column "Total?". For example, grand total 23,156,918 is pretty meaningless, since there are only 21,342,779 physically distinct members. And, of course, for such huge populations the chi-squares test insists that B is not representative for A, and C is not representative for B.                d <- data.frame(        cbind(          A = c(9817678, 1280567, 426856, 9817678),          B = c(621718, 191298, 79707, 701425),          C = c(90196, 30799, 11000, 87996)),        row.names = c('G1', 'G2', 'G3', 'G4'))           chisq.test(d[,1:2])     #     #   Pearson's Chi-squared test     #     #data:  d[, 1:2]     #X-squared = 160929.2, df = 3, p-value < 2.2e-16           chisq.test(d[,2:3])     #     #   Pearson's Chi-squared test     #     #data:  d[, 2:3]     #X-squared = 1539.55, df = 3, p-value < 2.2e-16      Am I applying the wrong statistical test here? Now, if I would randomly sample 100 members from A, B and C, and get the "average" numbers for G1 - G4, the chi-squared will be pretty confident that C is representative of B:                d.r100 <- sapply(d, function(x) round(100*x/sum(x)))      d.r100     #      A  B  C     #[1,] 46 39 41     #[2,]  6 12 14     #[3,]  2  5  5     #[4,] 46 44 40           chisq.test(d.r100[,2:3])     #     #   Pearson's Chi-squared test     #     #data:  d.r100[, 2:3]     #X-squared = 0.3943, df = 3, p-value = 0.9414      It feels wrong that the more data I process, the worse my representativeness becomes. After all, if studying "subpopulation C" will fetch me the answer for "population A" plus-minus mere 2-3%, it will be fine. But how do I estimate my confidence intervals for "A" if all I can study is "C"?


title: Questions about large-scale multiple testing

text: 1. From Wikipedia > Traditional methods for multiple comparisons adjustments focus on correcting > for modest numbers of comparisons, often in an analysis of variance. A > different set of techniques have been developed for "large-scale multiple > testing", in which thousands or even greater numbers of tests are performed. > ... Particularly in the field of genetic association studies, **there has > been a serious problem with non-replication** — a result being strongly > statistically significant in one study but failing to be replicated in a > follow-up study. Such non-replication can have many causes, but it is widely > considered that failure to fully account for the consequences of making > multiple comparisons is one of the causes. I was wondering what the problem with non-replication is, and why it can be because of the failure to fully account for the consequences of making multiple comparisons?   2. > **For large-scale testing problems where the goal is to provide definitive results** , the familywise error rate remains the most accepted parameter for ascribing significance levels to statistical tests.  > > Alternatively, **if a study is viewed as exploratory, or if significant > results can be easily re-tested in an independent study** , control of the > false discovery rate (FDR) is often preferred. What does "the goal is to provide definitive results" mean, and why is FWER most accepted in such cases? What does "if a study is viewed as exploratory, or if significant results can be easily re-tested in an independent study" mean, and why is FDR preferred in such cases?   3. > The FDR, defined as the expected proportion of false positives among all significant tests, allows researchers to identify a set of "candidate positives", of which a high proportion are likely to be true. The false positives within the candidate set can then be identified in a follow-up study. Since FDR "allows researchers to identify a set of 'candidate positives'", which I understand as the estimate of false positive individual nulls, why do we need further that "the false positives within the candidate set can then be identified in a follow-up study"? Am I right that FWER also can "allow researchers to identify a set of 'candidate positives'", as the set of individual tests that rejects their nulls. Thanks and regards!


title: Is the likelihood ratio test for the global null considered to be an omnibus test

text: Continue from my last question about an omnibus test. When the (global) null is the intersection of several individual nulls, I was wondering if the likelihood ratio test for the global null is considered to be an omnibus test. For more details, I quote from Casella and Berger's Statistical Inference: ![enter image description here](http://i.stack.imgur.com/wnMj0.png) ... ![enter image description here](http://i.stack.imgur.com/eeUOZ.png) Thanks and regards!


title: Is there an interactive notebook tool for R as iPython Notebook?

text: I've been using RStudio and I like it, but for the type of work I do a Notebook environment (like iPython Notebook) is more useful than an IDE. Is there an iPython Notebook equivalent for R?


title: Probability of a subsequence occurring simultaneously in two longer sequences

text: Suppose I have two sequences of characters of length L and M respectively, with the characters chosen i.i.d. from an alphabet of A through H, each with probability p=1/8. I want to find the probability that there exists a subsequence of length N which occurs simultaneously in both the length L and M sequences. I don't actually care whether it occurs exactly once or at least once in each of the sequences. Whichever is easier will do. For example, let's say we have (L) AABBCCDDEEFFGGHH (M) GBCCDFH If N=4, I'm going to search through M and discover that BCCD also occurs in L (exactly once in both in this case). I want to know how likely that was to have happened (finding any length four sequence, not specifically BCCD). In the case where M=N, I thought the answer would be something like $$\left( \begin{array}{c} L-N+1\\ 1 \end{array} \right)*(p^N)^1*(1-p^N)^{L-N}$$ giving a lower bound, but of this, I'm not terribly confident. In any case, M will probably be significantly longer than N, so it would be good to get a more accurate answer.


title: How to model month to month effects in daily time series data?

text: I have two time series of daily data. One is `sign-ups` and the other `terminations` of subscriptions. I'd like to predict the latter using the information contained in both variables. Looking at the graph of these series it's obvious that terminations are correlated with multiples of the sign-ups the months before. That is, a spike in sign-ups on May 10th, will lead to an increase in terminations in June 10th, July 10th and August 10th and so on, although the effect wears off. I'm hoping to get a hint as to which models I might employ to model this specific problem. Any advice would be much appreciated.. So far, I've been thinking a VAR model, but I'm not sure how to include the monthly effect - use a really high order of lags or add a seasonal component somehow?


title: Differences between Dwass-Steel-Critchlow-Fligner and Mann-Whitney U-test for a post-hoc pairwise analysis

text: If you have one control and four groups of patients with a non-normal small sample size, there are differences in these tests for a post-hoc analysis after Kruskal-Wallis?


title: Can a Naive Bayes classificator "learn" variables which are not in the training set?

text: # Theory I have three data sets, let's call them A, B and C. These data sets contain the following variables:               A = { x, A1, A2, ..., An }     B = { B1, B2, ..., Bm }     C = { A1, A2, B1, B2 }      x is the dependent variable which I'd like to predict. As you can see, only the data set A contains x. A and B are disjoint. On the other hand, the variables of C are the union of a subset of A and a subset of B. **The goal is to predict x from the B variables.** To achieve this, I thought of creating a Naive Bayes classificator for A. As soon as the classificator has learned to classify items based on the A variables, I'd include the data from C. At this point, the classificator doesn't have any knowledge of the additional variables coming from C (B1, B2), so they don't have any influence on the probability distribution of x. This way, the classificator 'extends its knowledge' onto the variables (B1, B2), and thus will be able to predict x from the B variables, without the presence of the original A variables. # Example Suppose the classificator is an e-mail spam filter: The x variable is the categorical variable `(Spam | Not spam)`, the A variables are the words of the e-mail, and the B variables are other e-mail properties like it's origin (sender IP address).   * The A dataset is a collection of e-mails where we know which one is spam, but we don't know the origin IP adresses.   * The C dataset is another collection, where we don't know which e-mail is spam, but we know the e-mail text and the origins. **We want to be able to classify the e-mails based on their origin.** First, we'd create a Bayes classificator based on the A dataset. Now we're able to classify the mails based on their text. Then, we'd extend the variables by adding the C dataset, which would 'teach' the Bayes classificator what IP addresses are the sources of spam e-mail. After this step, the classificator could detect spam just by looking at the e-mail origins. # Question **Is that possible, feasible? Are there better ways?** I suppose it kind of violates the requirements of the Bayes classificator, since the variables are not independent (the e-mail origins correlate with the words in the e-mails). But still, it seems like it could work.


title: Can you develop an econometrics model for stress test purpose only focusing on 2008-2009 data?

text: I have become aware that a group at a large corporation is developing an econometrics model to forecast sales of their product. They are using this model solely to estimate sales in specified stress test economic scenarios where they are given what the economic environment will be like, including real GDP contraction, rising unemployment rate, etc... out to 2016. Because of the nature of those scenarios, they think the most proper way to construct this model is to focus solely on the 2008-2009 period capturing the main period of the recent financial crisis. They have monthly data, so that gives them 24 monthly data points. Given that GDP's frequency is really quarterly, on this one variable it gives them only 8 true datapoints. But, they extrapolate it into 24 month observations. For the record, if they chose to, they have good internal data going back to 2001 and up to the current period. But, as mentioned they decided to focus instead solely on the 2008-2009 period. I will also answer this question as I have built many such econometrics models. And, I invite others to debate and rebutt my answer... and to post your own better answer.


title: Intuitive explanation for density of transformed variable?

text: Suppose $X$ is a random variable with pdf $f_X(x)$. Then the random variable $Y=X^2$ has the pdf $f_Y(y)=\left\{\begin{array}{ll}\frac{1}{2\sqrt{y}}\left(f_X(\sqrt{y})+f_X(-\sqrt{y})\right) & y \ge 0 \\ 0 & y \lt 0\end{array}\right.$ I understand the calculus behind this. But I'm trying to think of a way to explain it to someone who doesn't know calculus. In particular, I'm trying to explain why the factor $\frac{1}{\sqrt{y}}$ appears out front. I'll take a stab at it: Suppose $X$ has a Gaussian distribution. Almost all the weight of its pdf is between the values, say, -3 and 3. But that maps to 0 to 9 for $Y$. So, the heavy weight in the pdf for $X$ has been extended across a wider range of values in the transformation to $Y$. Thus, for $f_Y(y)$ to be a true pdf the extra heavy weight must be downweighted by the multiplicative factor $\frac{1}{\sqrt{y}}$ How does that sound? If anyone can provide a better explanation of their own or link to one in a document or textbook I'd greatly appreciate it. I find this variable transformation example in several intro mathematical probability/stats books. But I never find an intuitive explanation with it :(


title: Solving for a pdf of a function of a continuous random variable. Justification and reason for the procedure

text: Short version: When solving for the pdf of a function of a continuous random variable(say, $Y=X^2$), why can't you just plug in inverse of that($\pm\sqrt{x}$) into the pdf of the RV? Why do you have to start from the cdf of the $Y$, substitute $Y$ with $X^2$, and differentiate that to get the pdf? The point of this question is regarding the intermediate step involved.   1. Why is it not allowed to simple plug-in to the pdf of RV?    2. why is necessary and also justified to do so involving the cdf? Long version : Let's say you have a random variable $X$ and a function of the random variable, $g(X)$ (such as $Y=X^2$). You are interested in finding the pdf of $Y$, i.e. $f_Y(y)$. Usually, to do this, you have to start with the cdf of $Y$ first, and differentiate that to get a pdf of it. For example, $$F_Y(y)=P(Y\le y)$$ and here you substitute $Y$ with $X^2$: $$P(Y=X^2 \le y)=2F_X (\sqrt{y})-1$$ and you you differentiate this to get pdf of $Y$ $$d(2F_X(\sqrt{y})-1)/dy=f_Y(y)$$ My question is : WHY? why not just solving for $$f_X(\pm\sqrt{y})$$ what is wrong with this approach? and why is it justified to use cdf to solve this? Please explain this in a layman's term since clearly I am a novice in statistics


title: Probability density function for a logarithmic spinner

text: I'm having trouble deriving the PDF for the following problem (from the book _Doing Bayesian analysis with R and BUGS_ ): > Consider a spinner of the kind often found with board games, but with a > log10 scale from 1 to 100 (covering 0 to 360 deg). The spinner is fair, so > any value is equally likely. What is the pdf at point $x$? I know that the $x$ value at angle $t$ is given by $$x = 100^\text{t/360}$$ Doing some transformation $$x = \exp^{(t/360) \times \log(100)}$$ $$pdf(x) = pm(x)/dx$$ I also know $$dt/dx = 360/(\log100 \times x)$$ ## Question > How to derive the answer below, using the information above? ### Book answer > $$pdf(x) = 1/(2 \times ln(10) \times x)$$ Thanks.


title: Estimating Parameters of truncated generalized pareto distribution by probability weighted moments

text: Is it possible to estimate parameters for a truncated Generalized Pareto Distribution (Wikipedia Article) by probability weighted moments? If not, why is that so? I am aware that it would be possible numerically by MLE or similar methods, but I am specifically interested in whether a PWM approach is feasible.


title: Configural frequency analysis

text: What do you think about configural frequency analysis? http://en.m.wikipedia.org/wiki/Configural_frequency_analysis I personally love the idea to extract types and antitypes in heterogenous samples in an explorative approach. In Germany, however, this statistical procedure seems baraly taught in University context. Why is this technique not common? Are there practical reasons ( e.g. no need) or is the technique out-dated or suffers from other problems?


title: Gibbs sampling for correlated random variables

text: # Short summary Suppose two latent variables of a hierarchical model are correlated. Let $1-\epsilon$ be the degree of correlation. As $\epsilon\rightarrow 0$ the variables become perfectly correlated and Gibbs sampling appears to only be possible if the two variables are considered as one variable (i.e. using block-Gibbs). Is this indeed the case or am I lacking understanding? Edit: Wikipedia discusses this problem in their section on "Failure modes". However, there is no reference or resolution to the problem. # Complete question with example Consider a population whose members have two binary attributes: age $X$ and health $Y$. We encode the variables as $$X=\begin{cases}0&\text{if young}\\1&\text{if old}\end{cases}$$ and $$Y=\begin{cases}0&\text{if healthy}\\1&\text{if sick}\end{cases}.$$ Let the joint probability distribution be given by $$ P\left(X=x\cap Y=y\right)\equiv p_{xy}=\frac{1}{2}\left(\begin{array}{cc}1-\epsilon&\epsilon\\\epsilon&1-\epsilon\end{array}\right) $$ such that $\sum_{xy}p_{xy}=1$. Let us expand the probability distribution such that $$ P\left(X=x\cap Y=y\right) = P\left(Y=y|X=x\right)P\left(X=x\right), $$ where $$ \begin{align} P\left(X=x\right)\equiv p_x &=\sum_y p_{xy}=\frac{1}{2}\quad\forall x\\ P\left(Y=y|X=x\right)\equiv p_{y|x}&=\frac{p_{xy}}{p_x}=\begin{cases}1-\epsilon&\text{if }x=y\\\epsilon&\text{if }x\neq y\end{cases} \end{align} $$ Suppose we make a noisy observation of the state of health an individual $\hat{Y}$ with error rate $\delta$ such that $P\left(\hat{y}=y|Y=y\right)=1-\delta$. We construct a hierarchical model as shown in the figure below. ![](http://i.stack.imgur.com/Ujls0.png) To use a Gibbs sampler, we construct the conditional distributions (assuming a flat prior on $X$) $$ P\left(X=x|Y=y\right)=\begin{cases}1-\epsilon&\text{if }x=y\\\epsilon&\text{if }x\neq y\end{cases}$$ $$ P\left(Y=y|X=x\cap\hat{Y}=\hat{y}\right)\propto\begin{cases}1-\epsilon&\text{if }y=x\\\epsilon&\text{if }y\neq x\end{cases}\times\begin{cases}1-\delta&\text{if }y=\hat{y}\\\delta&\text{if }y\neq \hat{y}\end{cases}. $$ However, in the limit $\epsilon\rightarrow 0$ knowledge of $x$ gives us full knowledge of $y$ and Gibbs sampling is no longer possible because the sampler gets stuck: If $X=0$, then $P\left(Y=y|X=0\cap\hat{Y}=\hat{y}\right)=\delta_{y0}\quad\forall\hat{y}$ and vice versa. What is the best way to deal with this situation?


title: How to model a Time Series which has different nonstationary segments by ARIMA?

text: I had a monthly river temperature (408 values, separated 360 for modeling). Then I deseasonalized and transformed it to a normal time series by a plotting position technique. Now I need to fit an ARMA model to the time series which I got. These plots show the time series which I want to fit a model, ACF and PACF plots. ![enter image description here](http://i.stack.imgur.com/jtHXn.jpg) According to the results I tried fitting ARMA(2,0,0) and ARIMA(1,0,1) using `arima(TS, c(2,0,0))` in R. but both models were very bad fit. For example, AR(2) results are like this (I just plotted a portion of results for a better view) ![enter image description here](http://i.stack.imgur.com/uyfT1.jpg) The only explanation that I can think of is this: the results of my `KPSS` test and `ADF` test is stationarity for whole data, but If I apply the test on the different portions of data (like last 50 or 70 values) it will result in nonstationarity. **So, Is my conclusion correct? or there is something else wrong with my modeling? If I am right how can I model a time series which has different form of nonstationarity in different parts of the data?** Thank you, I attached the excel file of my time series in this link: time series data (csv)


title: Is the following property for positive random variables fulfilled in general?

text: [I have cross-posted this from math.stackexchange: http://math.stackexchange.com/questions/476466/is-the-following-property-for- positive-random-variables-fulfilled-in-general ] Suppose we have a continuous random variable $X$, defined on the interval $[0, \infty)$, which has density $f(x)$ and a finite expectation and variance. I am wondering whether the following is true in general: Take any value $\alpha > 0$, and let $x_1 := \mathbb{E}(X|X \leq \alpha)$, while $x_2 := \mathbb{E}(X|X > \alpha)$. Will a new, discrete random variable $\tilde{X}$, where $\tilde{X} = x_1$ with probability $F(\alpha)$, and $\tilde{X} = x_2$ with probability $1-F(\alpha)$, always have a lower variance than $X$? I have checked that this property is true for the uniform distribution. Formally, it is straightforward to see that the property is fulfilled if for all $\alpha > 0$, $\int_{0}^{\infty} x^2 f(x)dx \geq \frac{\left(\int_0^\alpha x f(x) dx\right)^2}{F(\alpha)} + \frac{\left(\int_\alpha^\infty f(x) dx\right)^2}{1-F(\alpha)}$. Using Jensen's inequality, a sufficient condition for the above is that $\left(\int_{0}^{\infty} x f(x)dx\right)^2 \geq \frac{\left(\int_0^\alpha x f(x) dx\right)^2}{F(\alpha)} + \frac{\left(\int_\alpha^\infty f(x) dx\right)^2}{1-F(\alpha)}$. However, I am not really sure where to go from there. The integral in the LHS can still be split up into the range $0$ to $\alpha$ and $\alpha$ to $\infty$, but doesn't seem to help me much. Maybe the property is not true in general, after all? Any input is appreciated. Many thanks in advance!


title: Probability output from support vector machine (svm) with soft margin

text: Based on my very simple understanding of SVMs, it seems like a probabilistic output would be a very useful thing to have. Soft margin seems to part of the way toward accounting for noisy data, but still classifies it into one class or another without saying anything about how likely a given point is to be classified correctly. I've found a couple proposed ways to get this with a simple Google search, but am curious what is the state-of-the-art and what is the most straight forward way of getting this probabilistic output. Thanks


title: Normalizing Term Frequency for document clustering

text: I have a problem understanding the normalization of Term Frequency weight in document Vector Space Model for clustering. Let's say that for document d I have counted occurences of all terms. I remember reading somewhere some time ago that they should be divided by maximum term frequency for that specific document, or sum of all frequencies, but I cannot find the source. Is that correct? If so, I would really appreciate getting some source for this information - I need it for my thesis. My second question is about calculating TF-IDF weight - in this case, should the TF be normalized somehow, or should I take 'raw' frequencies?


title: Fitting an exponential model to data

text: I have 2 variables, both from class "numeric": `> head(y)   [1] 0.4651804 0.6185849 0.3766175 0.5489810 0.3695258 0.4002567`   `> head(x)   [1] 59.32820 68.46436 80.76974 132.90824 216.75995 153.25551` I plotted them, and now I would like to fit an exponential model to the data (and add it to the plot) but I cannot find any info on fitting models to multivariate data in R! Only to univariate data, can somebody help? I don't even know where to start... Thanks!


title: How to prove it's not just a chance finding when testing several variables

text: If I tested several similar variables, such as socioeconomic status measures, but only a few were significant, how would I prove that it's not just a chance finding? Besides having the literature to back it up, is there any sensitivity testing people usually use to prove it's not just chance? To what degree would I need to consider it may be a chance finding? I mean, if 2 out of 6 variables are significant, should I be worried? Or it would apply to situations with a much higher number of variables?


title: KS test Enquiry in R?

text: I have a data vector $p$ and I want to test if its true distribution is uniform[0,1]. So I ran a 2 similar KS tests, one with $y = p$ and the other with $y = ecdf(p)$.               ks.test(runif(length(p), min = 0 , max= 1), p)      Please forgive me for not being able to show you guys the data vector p, because its length is 100000. Nonetheless, all entries in vector $p$ takes values between 0 and 1, i.e. [0,1] (inclusive). The result of the above test is:                   Two-sample Kolmogorov-Smirnov test          data:  runif(length(p), min = 0, max = 1) and p     D = 0.0082, p-value = 0.2463     alternative hypothesis: two-sided      but it comes with a warning message at the end:               Warning message:     In ks.test(runif(length(p), min = 0, max = 1), p) :       cannot compute correct p-values with ties      The 2nd KS test that was ran _at the same time_ (which I thought could avoid getting warning messages):               print(ks.test(runif(length(p), min = 0 , max= 1), ecdf(p)))      Results:                   One-sample Kolmogorov-Smirnov test          data:  runif(length(p) + 1000, min = 0, max = 1)     D = 0.0093, p-value = 0.008757     alternative hypothesis: two-sided      I have looked at other posts, but I still don't understand the warning messages. **Question 1:** What are the "ties" referring to ? Now, I checked out R's help and it said $y$ can be either a data vector or cdf of data vector. In the latter case, it leads to a one sample KS test, which I have no idea how it differs from 2 samples KS test. **Question 2:** What is the difference between one-sample VS two-sample KS test ? I thought the 2 KS tests I ran above are very similar (if not identical), since I just want to if the data vector follows a continuous uniform distribution. And I thought the 2nd KS test (using $ecdf(y)$) would avoid getting warning messages. If the 2 KS tests are actually doing different things, please tell me and advise which one I should use ? **Question 3:** The 2 KS tests yield different p values, the first one > 0.05, while the second < 0.05. But if I plot a histogram for my data vector, it actually looks uniform to me, which seems to suggest that there is something subtle about the 2nd KS test that I dont know about (any ideas ?): **Question 4:** Please advise me how to handle the warning message from the 1st KS test ? Can I just ignored them ? Note: I wrote a function that ran those 2 KS tests above and plot the histogram. All are done at once. ![enter image description here](http://i.stack.imgur.com/ikmcq.png)


title: Count regression model Results

text: I am trying to draw some conclusions about the fitting of one model, but after looking at some examples in the internet I just can't get a hold of it, the interpretation of the results I mean. Since all of the example reach a different conclusion with results I see really similar. I have some count data with high variability so I tried fitting a negative binomial model to the data and got the following results > summary(monb)               Call:     glm.nb(formula = Counts ~ Hour + weekday, data = modtab, init.theta = 0.2910141397,          link = log)          Deviance Residuals:          Min       1Q   Median       3Q      Max          -0.8156  -0.7714  -0.6790  -0.5941   2.5479            Coefficients:                        Estimate Std. Error z value Pr(>|z|)       (Intercept)       -0.520449   0.284584  -1.829   0.0674 .     Hour              -0.007665   0.013999  -0.548   0.5840       weekdayDonnerstag  0.044863   0.343570   0.131   0.8961       weekdayFreitag    -0.742229   0.365964  -2.028   0.0425 *     weekdayMittwoch   -0.448842   0.381636  -1.176   0.2396       weekdayMontag     -0.493662   0.353680  -1.396   0.1628       weekdaySamstag    -0.006994   0.336181  -0.021   0.9834       weekdaySonntag    -0.235636   0.343460  -0.686   0.4927       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Negative Binomial(0.291) family taken to be 1)              Null deviance: 383.88  on 633  degrees of freedom         Residual deviance: 375.59  on 626  degrees of freedom         AIC: 1067.4          Number of Fisher Scoring iterations: 1                     Theta:  0.2910            Std. Err.:  0.0453                 2 x log-likelihood:  -1049.4140       I think that the model does not fit the data properly but I really can say any valid arguments of why. Can you say what are the reasons for a bad fitting or if I am wrong what are the evidences for a good fitting by looking at this summary? Thank you


title: What is the main motivation behind using linear regression to model linear data as opposed to finding the means at 2 values and inferring a line?

text: After doing some tests I know the answer to this. If the standard deviation is high enough, the means will not converge fast enough to give an acceptable error. Assume the distribution of residuals is uniform and normal. If you take $n$ samples from each of the 2 points, you can take the "realm of possibilities" as some confidence interval with bounds: $$ME = z^*\frac{\sigma}{\sqrt{n}}$$ For whatever z-score is acceptable. The maximum error is going to be 2 times this (the bottom of one interval to the top of the other). If the relationship is already very small, then you must have large $n$ anyways, but there also is the fact that $n$ is proportional to $\sigma^2$ so if the standard deviation is high, then n must be extremely high to get a good answer. Is there a more official way to say this, or a resource I could look up that gives a formal reason/proof/explanation? Also, could this possibly invalidate the results of matching if the underlying relationship is linear and the standard deviation is high? Or is matching merely looking for _any_ answer, as opposed to a precise one?


title: How to chose best coin flipping machine?

text: How to chose best coin flipping machine ? On job interview I was asked to solve the following problem. There is some amount of different coins, some are fair and some are not. Also there are two coin flipping machines. As a result of several coin flipping experiments a table reflecting machine outcomes was created. Also a person was asked to flip the same coins and his outcomes were added to the same table. Resulting table has the the following format:               Coin Id | Flip cnt | P hds | P tls | M1 hds | M1 tls | M2 hds | M2 tls      where: P hds - count of heads flipped by Person   P tls - count of tails flipped by Person   M1 hds - count of heads flipped by Machine 1   M1 tls - count of tails flipped by Machine 1 M2 hds - count of heads flipped by Machine 2   M2 tls - count of tails flipped by Machine 2 Question: How to determine machine that is better to use for flipping coins using data in this table? As I understood at the interview a "better" machine in this context is the one that has coin flipping closer to normal distribution. Any ideas how to solve this task?


title: Difference between multinomial logit regression and multinomial logistic regression

text: Is there a difference between the two models - multinomial logit regression and multinomial logistic regression?


title: How to test for multicollinearity among dummy explanatory variables?

text: I have a Box-Cox regression where the explanatory variables are almost all dummy variables. If I want to see if there is multicollinearity among them, what would be an appropriate test? Do variance inflation factor (VIF) tests work here? What about the Pearson correlation coefficient matrix?


title: In Brownian motion, how does the absolute value of a particle's distance from the center scale with respect to number of iterations?

text: The mean distance is probably zero due to averaging out. But what about the mean of the absolute distance?


title: How to create indexes from "open" variables

text: I am running an impact evaluation to investigate the impact of a specific program on several indicators x. These indicators have values between 0 and up to 10000, and a high standard deviation. For my purpose, I'd like to create one index for each indicator each that takes into account the initial level of the indicator x in the base year (year=0) and the change between the base year and the post-intervention year (year=1), for each observation (household). The index should be constructed from two components, the initial level wi and the change diffi, with wi and diffi each being normalized (0 to 1) and wi X diffi making the indicator. At least, thats how it was done in a previous paper, where the variables, however, had logical limits, i.e. months per year etc. A total loss of the initial value between year 0 and year 1 (i.e. diffi=-100%) should be 0 and a gain of 100% or more (.e. diffi >= 100%) should be 1. I cannot simply focus on the change itself, as initial levels determine the "impact" of the change on the household (outcome variables are measures for food security in a drought-affected developing country). Do you have any suggestions?


title: Box-Cox transformation for residuals in R

text: I have residuals for my model. They are simply measured-predicted. However, I notice that they do not follow a normal distribution. I want to make my residuals distribution normal so that I can obtain confidence intervals that are representative of the distribution. I have tried doing the following:               >library(MASS)     >boxcox(residuals)      Error: $ operator is invalid for atomic vectors      However, I get an error. After looking at the Box-Cox more closely I have to input a formula or fitted object. If I input the measured vs. predicted fit, will the Box-Cox transformation normalize the residuals? Is there another way to implement the Box-Cox transformation so that it only needs to look at the distribution of the data?


title: Are direct Monte Carlo and Markov Chain Monte Carlo equivalent?

text: I'm not quite sure where to start... I'm trying to wrap my head around two methods for drawing samples from a parameter distribution, given a forward model and some distribution of the observable. I'll call these methods _direct_ and _Markov Chain_ Monte Carlo and describe them below. But to be clear on the parameters of my problem, I've got some observable $x$, with some normally-distributed error $\sigma_x$. e.g. I read a paper that said $x=1.2\pm0.1$ or whatever, presuming Gaussian errors. Now, I have a forward model that takes a parameter $\theta$ and returns a value of $x$. I'll call this a forward model $f$ and write $f(\theta)=x$. I think it's important to know that this is a horrible, difficult, complicated function that is not easily inverted. What I'd like to do is look at the distribution of the parameter $\theta$. Here are my two ways of doing so. **Direct MC** First, I construct a realization of the observation. i.e. I draw from the observed distribution. In practice, I would do this by inventing a number $1.2 + 0.1\times$ (a random number drawn from $N(0,1)$). Then, I do my best to invert the model, in practice probably with a gradient descent method or something similar. This produces an optimal $\theta$, which is one draw from the parameter distribution. To build a sample, I do this for a suitably large number of different (random and independent) realizations of the observation. **Markov Chain MC** Basically, I run an MCMC method. That is, I build a chain of $\theta$s, starting from some initial value and comparing a new value, drawn from a proposal distribution. The forward model is used to compute the probability of the new value, in terms of the distribution of the observations. i.e. $x$ is computed from $\theta$ and evaluated in $N(x,\sigma_x)$ or whatever. The new value is accepted if more probable than the current value, or randomly accepted even if worse. Otherwise the current value is kept. However the new value is assigned, it is added to the sample (the "chain") and the process is repeated. It's more detailed than this, but I'm not sure it's necessary to explain MCMC completely and accurately for the purpose of this question. Anyway, if the chain is "burned-in" and run long enough, it samples the distribution. So, finally, my question is: _are the two samples reproduced in these ways the same?_ i.e. _do the two methods produce the same results?_ I intend to do a small numerical experiment when I have some time. My understanding was that the two are the same, but I'm not sure why and can't prove it. Last, but not least, please let me know if I can clarify the question, because my experience on Stats.SE is that it's quite hard to make sure everyone is on the same page, and that I'm not very good at it!


title: Where to get started in data analysis

text: I wish to break into the field of data science, more specifically, analysing data sets, extracting features and classifying them or identifying relationships. However, I have no formal math background and am trying to re-learn it from the ground up. I've started an introductory stats and probability course on Coursera, and am looking at some linear algebra via Khan academy but looking at the math I'm still quite lost. I was wondering if anyone knew of a good introduction to data analysis text that was aimed at non-mathematicians, a "black box" approach, if you will. That would allow me to learn some basic techniques while I progress in my math that I could be getting on with alongside my studies. I've had a play around with MATLAB, sci-kit learn, scipy, numpy for Python but all the texts I use are very math heavy. Is this unavoidable? Should I just knuckle down with the math first? What would your advice be?


title: Statistical package that works with sql database

text: At my firm we have been calculating fundamental stock factors each week for over a decade. We have the performance of each of these factors. I'd like to enhance the statistical analysis of these factors and improve the multi-factor optimization. Is there a somewhat beginner "popular" software that interfaces with a sql database that would help me to perform this analysis? I understand this is probably a very elementary question here, thank you for your time.


title: Anomaly Detection Algorithm versus "just doing some statistics"

text: I've proposed using an anomaly detection algorithm in a project. The algorithm would consist of choosing some features we think might be indicative of anomalous examples. Then using a training set of anomalous and non-anomalous examples to fit parameters to a Gaussian distribution. We would then use these parameters to create a probability function p where p(x) < epsilon for some epsilon would indicate an anomaly. How would I answer this challenge? > “We love this concept, but how does this differ from us just doing some > statistics on a set of results?” By just doing some statistics, I assume the challenger means detecting the anomalies manually in old data and coming up with a set of conditions that would indicate anomalies in future data.


title: How to get prediction intervals at mean & at max of covariate values in R

text: I'm still working on R problems from a book, and using my spending data. **Part 1:** I need to predict the amount that a male with average data for status, income and verbal would spend along a 95% CI. (NB, male is coded as 0 and female 1 in the data set.) First I used the linear model `newspend<-lm(sepnding ~ status + income + verbal + sex)` and from a previous question the sex coefficient is for female, so how would I get the male?               > qt(0.975,42)     [1] 2.018082     > c(-22.12-2.02*8.21,-22.12+2.02*8.21)     [1] -38.7042  -5.5358      _First Question:_ I don't know how what this means and how will I get a 95% CI for males if this is for females? **Part 2:** Then repeat the prediction for males with maximal values of status, income and verbal and determine which CI is wider. Second, I used the max. values for the data in a new equation.                    sex             status          income            Min.   :0.0000   Min.   :18.00   Min.   : 0.600        1st Qu.:0.0000   1st Qu.:28.00   1st Qu.: 2.000        Median :0.0000   Median :43.00   Median : 3.250        Mean   :0.4043   Mean   :45.23   Mean   : 4.642        3rd Qu.:1.0000   3rd Qu.:61.50   3rd Qu.: 6.210        Max.   :1.0000   Max.   :75.00   Max.   :15.000                 verbal          spending           Min.   : 1.00   Min.   :  0.0        1st Qu.: 6.00   1st Qu.:  1.1        Median : 7.00   Median :  6.0        Mean   : 6.66   Mean   : 19.3        3rd Qu.: 8.00   3rd Qu.: 19.4        Max.   :10.00   Max.   :156.0        `x0<-data.frame(status=75.00, income=15.00, verbal=10.00, sex=1.00)`   I'm not sure if this is how I get the maximal values for a male? (Or should the 2nd part be a new question since it's a lot of information?)


title: R/Bayes: Why is a wildly different prior not affecting the posterior in bayesglm{arm}?

text: I am trying to see how different priors affect the posterior estimates in `bayesglm{arm}`. But no matter what the prior is, why is the posterior hardly changing?               library(Zelig); data(turnout) #Using the "turnout" data from the "Zelig" package          library(arm) #Loading the "arm" package     model.1 = bayesglm(vote ~ race + age, family=binomial(link=logit),                        prior.df=1, prior.scale=2.5,                        data=turnout) #This is the default specificication with a weakly informative prior      summary(model.1) yields:               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) 0.042320   0.176567   0.240  0.81058         racewhite   0.641959   0.134136   4.786  1.7e-06 ***     age         0.011222   0.003046   3.683  0.00023 ***      Now compare this to:               model.2 = bayesglm(vote ~ race + age, family=binomial(link=logit),                        prior.df=1, prior.scale=2.5, prior.mean=c(-100,0), #Specifying a crazy prior                        data=turnout)      Summary(model.2) yields:               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) 0.039891   0.176720   0.226 0.821410         racewhite   0.645076   0.134465   4.797 1.61e-06 ***     age         0.011219   0.003047   3.682 0.000231 ***      What's going on here? Surely a prior which is so different would influence the maximum likelihood estimates enough to sharply change the posterior?


title: Creating a correlated prior

text: I would like to create a "weakly informative" prior distribution for a couple parameters. They both could theoretically take any value between 0 and 1, but I have reason to think that they should be negatively correlated. Is there any sort of standard joint distribution (a beta-analogue to the multivariate normal?) where I could easily specify such a distribution? In looking around, I've read about the Dirichlet distribution, which seems to be almost perfect, but inappropriate because of it's requirement that the sum = 1.


title: Starting coefficient vector for GLM

text: I would like to know how R chooses its starting coefficient vector for a GLM when its `start` argument is left blank and defaults to `NULL`. For my personal implementation of a GLM, I have simply initialized $ \boldsymbol \beta_0 $ to be all $1$s. However, while this generally is fine, it can cause the iterative algorithm to diverge. Basically, I am just looking for a simple algorithm/formula that takes into consideration the data points and family of the GLM to choose the original coefficient vector, $ \boldsymbol \beta_0 $.


title: Variance unbiased?

text: I'm having some trouble with a homework problem, and would appreciate help. You have an i.i.d. sample of n observations drawn from a distribution.If $\hat x$, the mean, is known, we can estimate the variance with: $\hat \theta_x^2 = 1/n \sum(X-\hat x)^2$. Show that $\hat \theta_x^2$ is and unbiased estimator of the variance $E[(X-\hat x)^2]$. I understand that the sample variance is unbiased if we divide it by n-1, but why is it also unbiased, according to the problem, if we divide by n?


title: Is it possible to predict the likelihood of an order of random events?

text: Sorry if this is a n00b question, I'm just trying to wrap my head around the problem. I am trying to refute conventional wisdom here, so any help is greatly appreciated. And now for the question: For example, I flip a coin 100 times and I end up with 50 consecutive heads followed by 50 consecutive tails. How would I calculate the probability or likelihood that the next time I replay the scenario that I will end up with the exact same result? Would the equation be different if we were using a ternary or quaternary event instead of a binary event? Conventional, non-statistical, wisdom would assume that coming to the same result on consecutive attempts would be incredibly remote, but that the odds of eventually getting an identical result would increase as more attempts are made, i.e. 50/50, 50/50 would be implausible, but 50/50, 49/51, 51/49 ... 50/50 would be plausible.


title: How do you show that two populations are statistically similar?

text: > **Possible Duplicate:**   >  Testing hypothesis of no group differences Suppose I have $k$ samples from 2 independent experiments (service times by 2 methods) and their means are similar. How do I statistically show that both methods have similar service times?


title: How to test for Pearson correlation when one variable is fixed on date?

text: I have 4 groups of different respondents, each group surveyed on four different dates (points of time). All respondents have answered a psychological questionnaire related to their perception of death (Perception of Death Questionnaire, PDQ). I also have a sample of popular newspapers that have circulated on and before those four points of time. Let's assume, that I can build an death-presentation index out of the newspapers and let's call it "DPI". I want to test whether PDQ and DPI are correlated with each-other. * * * _My current insight into the problem follows_ : **One case = one date** I can average all answers to PDQ given at each day, and get 4-cases sample with columns PDQ and DPI. AFAIK I can't use classical correlation significance even if I had more than 4 of days: effectively my sample will be time series, where cases are not mutually independent. Calculating Pearson's correlation coefficient will be OK, but for testing its significance I will need a different set of statistical tests (like VECM). **One case = one participant** I could treat the sample as it is - as a sample of respondents. Each respondent has day of participation, and so I can pair it with the associated value of DPI. Then I can calculate Pearson's correlation coefficient. Since I have the same number of respondents for each group, I believe that the correlation coefficient will have exactly the same value as in "One case = one date" case. But can I test for it significance?? My guts tell me, that from the same reasons as above, I cannot. But the more I think about it, the less I know why exactly. Please, help ;-)


title: How to escape symbolic value in R

text: I just loaded a csv file in R. When I ran the `summary` command for one of the columns, I got the following:               > Error: unexpected symbol in "summary k_low"      I'm pretty sure I know what the 'unexpected symbol' is: for observations in which we had no reliable data for the variable `k_low`, we simply entered a period, or a `.` Is there a way to escape these observations in R? My boss uses STATA, which apparently escapes periods automatically, but I don't have STATA on my personal computer (and would prefer to use R in any case). So basically: Is there a way to get R to bypass any 'symbols' and return a summary of only those observations for which numerical data was entered?


title: Generalised linear model fitted values

text: I have ran this model in R:               glm(alert ~ water.height + ssp*ssp.zone + log(count) + ssp*days,          family=quasibinomial, data=ScanSampling_sub_alert)        * `Alert` = proportion of a goose flock that were alert,    * `ssp` = subspecific identity of flock,    * `ssp.zone` = which part of site they were in,    * `log(count)` = log of flock size, and    * `days` = days from 27.9.13 This is a plot of the the predicted proportion of flock that is alert conditional on `days` and `ssp.zone`: ![enter image description here](http://i.stack.imgur.com/JpFL3.jpg) Two things are strike me as slightly strange:   1. The top two lines for predicted values are curved, even though the model doesn't include a polynomial term.    2. The confidence interval for the top two lines in uneven and the predicted values are not always in the centre of the confidence interval. Can anyone explain why I'm seeing a curved line for predictions and an uneven confidence interval? Is it something to do with using `family=quasibinomial`? Unfortunately I cannot recreate these patterns using any built-in R datasets.


title: Admissible decision rule is always Bayes for some proper prior

text: I think I am able to visualize this graphically, but can't think of how to write this out. Some other literature online suggests drawing a picture, but I was looking for a concrete proof. I started working through contradiction, but could not get anywhere. Any help would be great.


title: Proportion Composite score for DV in rANOVA- Need Help Interpreting Results

text: In my design, I have 3 binary items for each condition (4 conditions total) meant to measure the same thing. All participants were in all conditions (repeated measures, 12 scores each). I made a composite score for each participants for each condition (yes, I know, not advisable but this is the way I am doing it so work with me, please). I coded the items such that 0= small and 1= large. I then interpreted the composite score to mean the likelihood of choosing the large item. (With me so far?) I ran a rANOVA using my 2 IV's as within subjects factors. The output says that there are 2 main effects and an interaction (I can't post the image but I am basing this off of the table called "multivariate tests" with Wilks' Lambda, which is p<0.05 (as are all the other tests there). However, when I look at some of the levels (either in the main effects or in the interaction), I see that the confidence interval contains 0.5 (which to me means that we the level does not provide us any more info about how likely it is to choose "large.") Does the significance in the table I have shown here mean that the effects do exist? Should I interpret the effect by looking at the individual levels? (i.e. <0.5 means more likely to choose large) and see if the confidence intervals overlap? What does it mean if the table says a main effect with 2 levels is significant but one of the 2 levels has a confidence interval that contains 0.5? Thanks so much!


title: Compact letter displays with non-estimable cases

text: Compact letter displays are similar to grouping bars often used to display the results of a multiple-comparison procedure, but can be used more generally. You give a list of sample means (or other statistics) and associate one or more letters with each mean. Two means sharing the same letter are not significantly different according to the procedure used. Two means having no letters in common are significantly different. My question is: **What should one do if some of the means are not estimable?** (perhaps because of missing cells in the design) From one point of view, the non-estimable case(s) should get every letter used, because obviously you haven't found it to be statistically different than the others. So the display could look like this:                treat    lsmean      SE df  lower.CL upper.CL .group      A      6.666667 3.25747 18 -0.177024 13.51036  1          B      8.500000 3.25747 18  1.656309 15.34369  1          C     14.333333 3.25747 18  7.489643 21.17702  12         D            NA      NA NA        NA       NA  12         E     25.500000 3.25747 18 18.656309 32.34369   2         The other point of view would be that this unnecessarily clutters the display. So you should show this:                treat    lsmean      SE df  lower.CL upper.CL .group      A      6.666667 3.25747 18 -0.177024 13.51036  1          B      8.500000 3.25747 18  1.656309 15.34369  1          C     14.333333 3.25747 18  7.489643 21.17702  12         D            NA      NA NA        NA       NA      E     25.500000 3.25747 18 18.656309 32.34369   2         So what do you all think? Assign letters to the `NA` results, or not? ## footnote There are also cases where the differences are estimable, but the means themselves are not. So it is possible to get a display like this:                treat lsmean SE df lower.CL upper.CL .group      A         NA NA NA       NA       NA  1         B         NA NA NA       NA       NA  1         C         NA NA NA       NA       NA  12         D         NA NA NA       NA       NA  12         E         NA NA NA       NA       NA   2         Note that if scheme 1 was used in the above example, we cannot tell whether the differences involving `C` or `D` are estimable or not. Of course, in such a situation, we really would be well advised to display the actual estimates of the differences instead!


title: Interpretation of crossed random effect interactions in lme4

text: I'm considering a model in lme4 in which I am estimating random effects for two crossed factors, very similar to the Machines example in Bates' 2010 draft book (http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf; Chapter 4, section 4.1.2). Bates outlines two different ways to model the interaction between the crossed factors, and although I see the technical difference between them, I am not entirely confident that my substantive interpretation is correct. The first model is:               fm11 <- lmer(score ~ Machine + (Machine|Worker), Machines, REML=FALSE)      This model yields, for each worker, an intercept and a set of random effects for each machine:               ranef(fm11)       $Worker         (Intercept)    MachineB   MachineC       1   0.3188313   2.2294816  0.6119687       2   0.1784391  -0.9913178 -4.4293996        ... The other model is:               fm12 <- lmer(score ~ Machine + (1|Worker) + (1|Machine:Worker), Machines, REML=FALSE)      which yields:               ranef(fm12)         $`Machine:Worker`            (Intercept)       A:1  -0.7468217        A:2   1.5456840       A:3   1.7697364        ... as well as an intercept for each worker. Would it be correct to say that in the first model, the random effects are a function of Machine (that is, we need a different random effects variance parameter for each Machine, and some Machines might have more variance in their random effects than others), and in the second model, the random effects are a function of idiosyncratic combinations of Machine and Worker (that is, there is variability in the random effects for different combinations of Machine and Worker but that variability doesn't vary as a function of Machine; one general random effect variance parameter covers it all)? I suppose the conclusion I would like to draw is that if Model 1 fits better than Model 2, then the amount of random effects variance is different for different Machines (we need to specify which Machine we are talking about in order to model the effects properly), whereas if Model 2 fits as well as Model 1, then the random effects variance doesn't vary meaningfully by Machine (there is no added explanatory value in modeling separate random effects variance parameters for different Machines). Does that seem like a defensible interpretation, or am I reading too much into the existence of those extra parameters?


title: Alternative to MAPE when the data is not a time series

text: I have a data set where many of the actual values are zero, so I can't use MAPE. It's not a time series, so I can't use MASE ala our very own Rob Hyndman. Is there another alternative to MAPE that I could use?


title: TukeyHSD with factorial ANCOVA

text: I'm trying to run _post-hoc_ comparisons on a factorial ANCOVA model for a project of mine and I'm having some difficulty. First off, lets start with the experimental design:   1. Factors include: Focal plants infection status (2 levels), Competition (2 levels) and Damage (2 levels) crossed in a fully factorial design (i.e. 8 treatments)   2. Response variable is simply total biomass of focal plants.   3. Competitor plant biomass and Initial focal plant biomass are being included as continuous covariates. The model that I ran was as follows:               model1<-aov(TotalFB~Focal*Competition*Grazing+TotalIB+TotalBB, data=dat2)      I realize I cannot use the TukeyHSD function as it does not handle continuous covariates. I have played around with the ghlt from the multcomp package and cannot seem to get it to work correctly. I am specifically interested in the significant Focal*Grazing interaction and was wondering how I could see if there are significant differences within the levels of focal plant infection status but across the levels of grazing. Any help would be greatly appreciated. Feel free to ask me if you require additional information. Thanks


title: Acceptable to reverse score a principal component?

text: I have two datasets from similar psycholinguistic experiments. In both of them, information about the participant's reading and spelling ability was collected, then converted into standardized scores zRead and zSpell. The aim is to use these as covariates when examining the priming effects from the experiments. Because these predictors are correlated, and because theoretically it's useful to differentiate their shared variance from their unique effects, I'm looking into using principal components instead, i.e. generate two principal components from the zRead and zSpell scores and use these as the covariates in a linear mixed model. The problem is that although the principal components come out similarly in each experiment, with PC1 indexing the shared variance and PC2 differentiating reading and spelling ability, the sign of PC2 is different in each dataset, i.e.: Experiment 1 correlation table:                      zSpell zRead  PC1a   PC2a     zSpell  1.000 0.504 0.867 -0.498     zRead   0.504 1.000 0.867  0.498     PC1a    0.867 0.867 1.000  0.000     PC2a   -0.498 0.498 0.000  1.000      Experiment 2:                      zSpell  zRead  PC1a   PC2a     zSpell  1.000  0.485 0.862  0.508     zRead   0.485  1.000 0.862 -0.508     PC1a    0.862  0.862 1.000  0.000     PC2a    0.508 -0.508 0.000  1.000      The aim is to present both of the related datasets together, so explaining how PC2 represents a (slightly) different thing in each dataset might be confusing. Is it acceptable to reverse score individual components by multiplying by -1? If not, is it okay to reverse all the components at once if the interpretation would make more sense? I can't see this changing the correlational structure of the variables, but I'm not sure if there other reasons to avoid it.


title: How to use the datasets to practice machine learning algorithms?

text: I'm a beginner in ML, and I want to practice with some algorithms I learned from the book, after searching stats.SE, I found some places where I can get some datasets. But I have this stupid question, how to use them? For example, the Iris dataset from UCI contains 150 examples, 50 instances for each classes. OK, I thought I should take bulk of this dataset as **training** examples, spare rest of the dataset as **test** examples, right? And the problem is how should I divide this dataset? Shuffle them randomly and then draw the training examples randomly? I think different training set will affect the model learned from it, right? Please give me some advice on how to use the dataset, thanks.


title: Regression analysis?

text: I have scored pig carcasses for injuries at 2 points; immediately after slaughter (slaughter stage 1 [SS1]) and again after the carcases have been washed, scalded and dehaired (slaughter stage 2 [SS2]). I have used an ordinal scale to record damage - absent, mild, moderate and severe skin damage. I have reason to believe that mild skin lesions will be less visible at SS2 but severe skin lesions will actually be more visible at SS2. How can I test this? Do I need to do some sort of regression analysis? Thanks


title: Autoregressive prior distributions

text: I have a set of 50 (`P`) time series' and each one is on a slightly different but should be influenced by the same effects of time. To get the time effects I have constructed the following model random intercept GLM with a Poisson Link Function (below). Since I think observations from day to day are likley to be correlated, I put a prior on the coefficient for each day (`betastar[d]` on `beta[d]`) centered at the previous day. I haven't really seen this done before so I wanted to see if anyone sees any serious issues with this fairly simple model structure. Here is some jags/bugs code to explain:               model {          #prior on random intercepts     for (p in 1:P){         alpha[p] ~ dnorm(0,0.001)     }          for (d in 2:D){         betastar[d] ~ dnorm(beta[d-1],0.001)         beta[d] ~ dnorm(betastar[d],0.001)         }          beta[1] ~ dnorm(0,0.001)          #likelihood     for (i in 1:PDs){         y[i] ~ dpois(exp(alpha[lookupP[i]] + beta[lookupD[i]]))     }}


title: Exact central confidence interval for a correlation

text: I am working my way through (T.J. Diccio & B. Efron, "Bootstrap Confidence Intervals", _Statistical Science_ , 1996, _11_ (3), 189–228), and I'm stuck even before I get to the good stuff. In the introduction, there is a table of before/after treatment values called `cd4`, and the correlation $\hat\theta$ is given as 0.723. The paper then states, "We can find an exact interval for $\theta$ if we are willing to assume bivariate normality ... The exact central 90% interval is (0.47, 0.86)." Here I stipulate that I am a stats moron. When I key in the data and compute a confidence interval by hand in R, either using the Fisher transform or the _T_ distribution, I don't get that interval. How does one compute the interval they are talking about?               > cd4        subj base oneyr     1     1 2.12  2.47     2     2 4.35  4.61     3     3 3.39  5.26     4     4 2.51  3.02     5     5 4.04  6.36     6     6 5.10  5.93     7     7 3.77  3.93     8     8 3.35  4.09     9     9 4.10  4.88     10   10 3.35  3.81     11   11 4.15  4.74     12   12 3.56  3.29     13   13 3.39  5.55     14   14 1.88  2.82     15   15 2.56  4.23     16   16 2.96  3.23     17   17 2.49  2.56     18   18 3.03  4.31     19   19 2.66  4.37     20   20 3.00  2.40     > r = cor(cd4$base, cd4$oneyr)     > r     [1] 0.7231654      Fisher transform:               > se = 1/sqrt(17)     > se     [1] 0.2425356     > tanh(c(r-1.6448*se, r+1.6448*se))     [1] 0.3133382 0.8082940      _T_ distribution:               > sr = sqrt((1-r^2)/(20-2))     > sr     [1] 0.1627936     > tc = abs(qt(0.05/2, 18))     > c(r-tc*sr, r+tc*sr)     [1] 0.3811486 1.0651821


title: data mining/predictive modelling methods for small data sets

text: I have a small timeseries dataset that has 30 records with 6 predictor variable and 1 response variable. I would like to regress my time series response with 6 predictor variables. I have been using traditional statistical methods to do modeling. My question is are there any good data mining predictive modelling methods (example: neural networks/svm/decision tree) for smaller data sets as mentioned in the example above ?


title: How to calculated Confidence Interval for autocorrelated and lognormally distributed data?

text: My data is autocorrelated and is lognormally distributed, how can I calculate Confidence interval of that set of data?


title: Random forest - binary classification vs. regression?

text: I have a dataset that I'm trying to classify into 2 groups, A and B, using a random forest model. I know the true grouping and I'm trying to see how well I can model it using the other available variables. I've tried 2 different approaches that I thought would be equivalent, but which are actually giving me quite different results:   1. Reading in the grouping as a (non-numeric) factor in R, growing a classification forest, and taking the proportion of trees that vote for group A as my prediction.   2. Constructing an indicator variable for membership of group A, growing a regression forest, and taking the ensemble prediction as usual. The split between the 2 groups is roughly 90-10 A vs. B. I'm growing 240 trees from ~200k observations of the same variables. I've left most of the settings at the defaults for the R randomForest package, but to keep the processing time down to a manageable level I've increased the node size to 200. The results are as follows:   1. In the vast majority of cases, all 240 trees vote for A. The average predicted chance of any one observation being in A is about 99.9%. Worse still, not a single member of group B gets a majority of votes for group B!   2. I get a wide range of predictions, with the mean prediction lying close to the observed mean of ~90%. How can two apparently similar methods give such different results? As for how I ended up trying this - I was initially trying to classify my dataset into a larger number of groups, of which B was one, but I noticed that B was being classified almost 100% incorrectly. The other groups are all much better behaved, even though most of them make up a far smaller proportion of my data.


title: Is there overfitting in this modellng approach

text: I recently was told that the process I followed (component of a MS Thesis) could be seen as over-fitting. I am looking to get a better understanding of this and see if others agree. **The objective** of this part of the paper is to   * Compare performance of Gradient Boosted Regression Trees against Random Forests on a data set.   * Look at the performance of the final model chosen (either GBM or RF). The `gbm` and `randomForest` packages in R are being used, along with   `caret`. **The process** followed was as follows:   * Preliminary pre-processing of the data (e.g. plugging missing values of the nominal predictors with a distinct category called "Missing"). The target variable was not looked at in regard to any pre-processing (which was very minimal).   * Create a grid of values for the meta-parameters of each algorithm (e.g. number of iterations for GBM).   * Create 25 random splits of the data set (65% training and 35% test). **Repeat 25 times the following for GBM** (Each time utilizing one of the random train/test splits. Each time, which training and test set are "current" of course changes - this is repeated leave-group-out cross validation):   * Use 5-fold cross validation to find the "optimal" parameter settings of the algorithm over the grid search. Nothing from prior runs used at all in the current run.   * Once determined, fit a model to the full "current" training set and predict the "current" test set. Set aside the performance measure of this run. Once 25 performance measures (actually a domain specific measure, but think of it as accuracy) are obtained in this fashion, follow the exact same process, using the exact same independent samples of train and test, for RF (same process, just with different grid search of course). Now,I have 25 measures of performance from the then "current" test sets for GBM and RF. I compare them using a Wilcoxon Signed Rank Test and also a permutation test. I found GBM to be superior. I also claimed that the distribution of the performance measure from these 25 runs for GBM is the expected performance of the final GBM classifier. What I did not do, was to pull out a random test set from the very beginning and set it aside to be compared to the final GBM model built from all the training data. I contend that what I did was actually much better as I repeated the split data / tune model / test on hold out process 25 times versus only once. Is there over-fitting here? Since the 25 runs were used to select GBM versus RF does that mean that the performance measures acquired from the process can not be used as performance estimation for the full model? **EDIT** In response to Wayne's comment, here is what was done during each of the 25 runs:   1. The data sampled for the ith training set (i=1,..,25) was broken up into 5 equally sized groups. A model was fit using 4 of the 5 groups, setting the GBM parameters (e.g. number of iterations) equal to the values in the jth (j=1,..,18) grid.    2. The performance on the 5th group was calculated using this model.   3. Steps 1 and 2 were repeated 4 more times (regular old k-fold CV with k=5). The performance was averaged from the 5 sub-runs and this made up the expected performance of GBM with that certain set of parameter values.   4. Steps 1 -3 were repeated for the 17 other "rows" in the grid. Once completed, the best parameter values from the above exercise were determined and a GBM was fit using these parameter values and the full ith training set. Its performance was estimated on the ith test set. Once this whole process was done 25 times, there were 25 measures of performance available for GBM. Then they were gathered for RF in the exact same way. After comparing and choosing GBM, I looked at those 25 performance measures and took the mean and Stnd Error to determine a confidence interval for a GBM model on this data.


title: Nonparametric repeated measures: all versus control test

text: https://www.youtube.com/watch?v=gI0xQcLIplA In the video link I have attached, a particular Friedman's test is implemented for a comparison of several treatments against a control. I am not sure if this implementation is well known. My knowledge is that Friedman compares several treatments in an all-to-all fashion. Also, I do not know if this tests maintains the Type I error constant. This proposed test would be a non-parametric version of a Dunnett's test, with repeated measures. Particularly, I wanted to know:   * Are there open source references about this test? the video has a website reference (http://appliednonparametricstatistic.blogspot.com), and there is a cited book there called "Applied Nonparametric Statistics" by Wayne A. Daniel http://www.amazon.com/Applied-Nonparametric-Statistics-Wayne-Daniel/dp/0534381944. I guess that more explanations about publications could be there. Any guiding content with more references would be useful.   * Is the Type I error controlled by using this test on a comparison of a higher number of treatments against a control? Regards, Ignacio


title: Influencing scores with level of uncertainty when sample is small

text: I am doing a thesis research on fraud in a review dataset and have a hard time solving the following problem. I created a feature that calculates the disagreement between two types of users (new and frequent reviewers). I arbitrarily chose to require at least 2 reviews from each type of user before I can say something about how they disagree with a score for a restaurant. Because the required amount of 2 reviews is not scientifically chosen I desire a solution that calculates a score that is influenced by the size of the samples. So if both types of users have 5+ reviews, I want the score to reflect that it is more trustworthy because of the larger samples. Vice versa, if the score is calculated on two small samples it should reflect that it is less "useful". My goal is thus to end up with one score that is already influenced by the level of trustworthiness. If you have any suggestions please try to explain it in as simple terms as possible, if not obvious already, I am by far not an expert in this field. Could not find anything on the forum but I really don't know what to search for.


title: Confidence interval to estimate the quality of a set of observations by matching - nonmatching ratio of second observations

text: I'll ask my question as if it was a school assignment. John makes 400 observations. These observations can be of 5 types: A, B, C, D, or E.                   observation 1:   A         observation 2:   E         observation 3:   A         observation 4:   D         observation 5:   C         observation 6:   A         observation 7:   B         ....         observation 399: A         observation 400: E      Marie re-makes 45 randomly chosen of the 400 observations. 38 of them match John's observations and 7 of them do not match John's observations. In other words, 15% of all Marie's observations do not match John's observations. What is the 95% confidence interval of this percentage ? Let's assume the answer is: 5% - 25% (or 15% ± 10%) Does it mean that if 40% of John's observations are "A", I'm actually sure at 95% that there is between 40+40*0.1 (=44%) and 40-40*0.1 (=36%) of observations that really are of type A?


title: Many multiple linear models with predictors from a larger set

text: I wish to perform multiple linear regression of 5 or 6 predictors at a time. I have total no. of 50 predictors. I want to find regression coefficients for all possible combinations of predictors from the master dataset of 50 predictors taking 5 and 6 predictors at a time. $$Y = a.X_1 + b.X_2 + c.X_3 + d.X_4 + e.X_5 + f$$ $$Y = a.X_1 + b.X_2 + c.X_3 + d.X_4 + e.X_5 + f.X_6 + g$$ I also want to find out the correlation coefficient between original Y and predicted Y after multiple linear regression. How can I perform this using matlab ? My purpose is to find correlation coefficient value greater than 0.85 for between Y data and different combinations of X data (master set 50). can you suggest any alternative to multiple linear regression by which I can achieve such r value ? I am a beginner in both statistics and MATLAB. Hence, I am randomly trying all possible combinations.


title: Modeling conditionally independent observations using logistic regressions

text: I'm interested in modeling the probability of successfully arriving at a spawning site for an individual $i$, given two impediments that are conditioned on one another. I know whether an individual made it pass hurdle 1 ($y_{1,i}$) and hurdle 2 ($y_{2,i}$), and several other measurements that I want to use to model said probabilities with. I'm interested on the effect of these variables on the two (assumed conditionally independent) hurdle probabilities. $$ Pr(Y|y_{2,i}=1,y_{1,i}=1,X) = Pr(y_{2,i}=1|y_{1,i}=1,X)Pr(y_{1,i}=1|X) $$ Is is appropriate to use two separate, independent logistic regressions to model this system (assuming that the processes describing the two hurdles are different)?


title: In a meta-analysis, how should one handle non-significant studies containing no raw data?

text: Let's say that I'm conducting a meta-analysis, looking at the performance of group A and group B with respect to a certain construct. Now, some of the studies that I'll come across will report that no statistical differences could be found between the two groups but no exact test statistics and/or raw data will be presented. In a meta-analysis, how should I handle such studies? Basically, I see three different alternatives here:   1. Include them all and assign to each one of them an effect size of 0.   2. Throw them all out.   3. Do some kind of power analysis for each one of them or set a threshold at a certain number of participants. Include all which should have been able to reach statistical significance and assign to each one of them an effect size of 0. Throw the rest out. I can see merits with all the different options. Option one is fairly conservative and you'll only risk making a type II error. Option two raises the risk for making a type I error, but it also avoids having your results ruined because of a bunch of underpowered studies. Option three seems like the middle road between option one and option two, but a lot of assumptions and/or pure guesses will have to be made (What effect size should you base your power analyses on? What number of participants should you demand from each study for it to pass?), probably making the final result less reliable and more subjective.


title: Maximum likelihood of two Poisson distributions

text: Suppose I have two sets of samples $\{a_i\}$ for $i=1, \ldots, n_a$, and $\{b_j\}$ for $j=1,\ldots, n_b$, drawn respectively from two Poisson random variables $A \sim\textrm{Poiss}(\lambda_a)$ and $B \sim \textrm{Poiss}(\lambda_b)$ with the constraint $$ \frac{\lambda_a}{c_a}=\frac{\lambda_b}{c_b}=\mu, $$ where $c_{a}$ and $c_{b}$ are known constants. What is the analytical expression for the maximum likelihood of $\mu$?


title: Find minimum number of batteries to last a week with 95% probability given mean and sd of lifetime

text: A TA needs a very important program on his computer but his power cable is missing. Luckily, he happens to have 100 batteries lying around and a contraption that can store batteries to use one at a time. From his understanding, the average lifetime of a battery will last 2 hours for his computer with a standard deviation of 30 minutes. How many batteries will he need to place in the contraption to keep his computer continuously running for the next 7 days (168 hours) with at least a 95% probability? I'm trying to approach this question using confidence intervals, but i'm confused what the sample proportion would be, would it be 84 batteries? cause u need 84 to run the computer for 168 hours, assuming each battery lasts exactly 2 hours. Would my standard error just be the standard deviation given? Or is there another way to approach this question


title: How to compare frequency distributions between 3 groups

text: I have characterized cells into 6 categories and want to compare the percent distribution of the different categories at 3 time points. Once I figure that out, I actually have additional treatment groups as well. What is the appropriate statistical test for comparing the frequency distribution between the 3 timepoints? Also, I am not a programmer, but hoping to do this in GraphPad.


title: What's the sample size required to confirm that two treatments are identical?

text: Goal: Determine if a machine is defective Say I am going to give the same treatment to a group of people. I have a machine that will randomly split them into $2$ groups before giving the same treatment ($50\%$ chance of being in either group). Treatment has a known success rate of about $20\%$. So, chance of a successful treatment in either group is $0.2\cdot 0.5 = 10\%$. I'm trying to test the machine for an unknown bias in how equally it administers treatment and how randomly it splits the people. So if I have $100$ people. I would expect $50$ people to end up in each group and $10$ in each group to show successful treatment. If I test fewer people, there is a higher chance of success rates' being different for the two samples due to chance and sample sizes being different. But with more people, the sizes of the two groups and the success rates for the two groups should approach the expected numbers **if** the machine is working properly. How many people do I need to test to see if the success rates are different enough that I can conclude the machine is defective with $90\%$ certainty? Should I use power analysis to determine the sample size needed to rule out effects greater than say $1\%$ relative difference? It seems a power analysis would dictate an infinite sample size to rule out an effect close to zero. Is there an assumption in power analysis that is violated when used in this situation? Would I adapt some method like this?: http://en.wikipedia.org/wiki/Checking_whether_a_coin_is_fair#Examples Thank you.


title: Paired t-test of medians

text: I need to run an inferential test on some large data where the individual data points have a heavily skewed distribution. I'm considering doing a paired t-test across a number of days comparing the results of the test group to the control group. But if I take daily means, I'm still subject to outliers influencing those values. What if I take daily medians and run the paired t-test in the same way? Is there anything "wrong" with this approach? Something tells me there is, primarily because I don't know how medians "behave". Does the CLT apply to medians also? Should they be normally distributed?


title: How to wrap up individual One sample t-test results in a single one, in SPSS?

text: I'm evaluating a questionnaire about job satisfaction, based on Likert scale. As for any questionnaire, there are some correlated questions assessing one factor. Since I’m using SPSS, I know that I can evaluate the population’s average of “each question” with One Sample t-test; However, what I need is an overall answer. As you can see in the table, the t-test result of each question is determined, but how can I make the final conclusion? How much will the corresponding overall average be? (Test value is 3) my one sample t-test result table


title: Error while mapping SYMBOLS to ENTREZID

text: I am getting a strange error converting Gene Symbols to Entrez ID. Here is my code:               testData = read.delim("IL_CellVar.txt",head=T,row.names = 2)     testData[1:5,1:3]      ClustID Genes.Symbol ChrLoc NM_001034168.1 4 Ank2 chrNA:-1--1 NM_013795.4 4 Atp5l chrNA:-1--1 NM_018770 4 Igsf4a chrNA:-1--1 NM_146150.2 4 Nrd1 chrNA:-1--1 NM_134065.3 4 Epdr1 chrNA:-1--1               clustNum = 5     filteredClust = testData[testData$ClustID == clustNum,]     any(is.na(filteredClust$Genes.Symbol))      [1] FALSE               selectedEntrezIds <- unlist(mget(filteredClust$Genes.Symbol,org.Mm.egSYMBOL2EG))      Error in unlist(mget(filteredClust$Genes.Symbol, org.Mm.egSYMBOL2EG)) : error in evaluating the argument 'x' in selecting a method for function 'unlist': Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs Another approach fails too:               selectedEntrezIds = select(org.Mm.eg.db,filteredClust$Genes.Symbol, "ENTREZID")      Error in .select(x, keys, columns, keytype = extraArgs[["kt"]], jointype = jointype) : 'keys' must be a character vector Just for the sake or error, removing 'NA', doesn't help:               a <- filteredClust$Genes.Symbol[!is.na(filteredClust$Genes.Symbol)]     selectedEntrezIds <- unlist(mget(a,org.Mm.egSYMBOL2EG))      Error in unlist(mget(a, org.Mm.egSYMBOL2EG)) : error in evaluating the argument 'x' in selecting a method for function 'unlist': Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs I am not sure why I am getting this error as the master file from which gene symbols were extracted for testData gives no problem while converting to EntrezID. Would apprecite help on this. Thanks AK


title: Publication bias - trim and fill - macro/syntax

text: I'm looking for an SPSS or Excel macro/syntax that will allow me to compute Duval and Tweedie’s trim and fill procedure (see paper here: http://onlinelibrary.wiley.com/doi/10.1111/j.0006-341X.2000.00455.x/abstract). I have a relatively large number of variables which would make this very time consuming to do by hand. Although I have found some macros for R I'm looking to see if something is available in SPSS or Excel.


title: Rules of thumb for PERMANOVA

text: I'm currently assessing the best option for looking at some community composition data, and determining how they differ among locations. Ultimately, as I'm simply looking to see if there is a difference depending on site (4 level factor), I was looking at Permanovas. Does anybody have any good references which explains the assumptions, such as maximum number of response variables (e.g. how large can the species composition be)? Many thanks.


title: Overfitting times series data

text: I am trying to use Robert Shiller's stock market data to estimate future 5 year returns. My first naive attempt used M5P in Weka (decision tree regression) with 60% of the data for training and 40% for validation. This gave a model with 120 rules and an excellent fit, but it is clear that the validation set is "nearly twinning" the training set, and M5P is overfitting. How do I choose a model in these circumstances? Sampling the data would reduce the model complexity, but what is the rationale for that and how do I determine sample size? Intuition tells me to just fit a linear relationship between CAPE, date and return, but I'd like to operate on a slightly more formal basis. **Edit 1:** There has to be a difference between the training and test set. What are the metrics for this difference?


title: What is the PDF of $[(X-a)^2 + (Y-b)^2]^{1/2}$ where $X$ and $Y$ are two non-standard normal random variables?

text: I have to conduct an experiment getting data from a system. These data are the estimated values, provided by the system, of a true value that we know beforehand. I then compare the estimated values with the true value to analyse how accurate the estimations of the system are. In this context, we have a set of estimated values and a true value. Each of the estimated values and the true value are represented in terms of $(x,y)$ coordinates. Assume that the coordinate of the true value is $(a,b)$ where $a$ and $b$ are constants. Also assume that the $x$-coordinate and $y$-coordinate of the estimated values are two normally distributed random variables with different means and variances, that is $X \sim \mathcal N(m_1, v_1)$ and $Y \sim \mathcal N(m_2, v_2)$. My question is: Is there any parametric method to calculate the probability of the random variable $$Z= \sqrt{ (X-a)^2 + (Y-b)^2 }$$ @update: Thank you very much for your help. After analyzing the estimated values, I have learned that the covariance between X and Y is not zero, which is: $$cov(X,Y) \neq 0$$


title: automatic assign class name based on text

text: My question is , I have a set of plain text , i want to create category based on the text. Eg: i have written something about Soup recepie then the algorithm must create a category called Food. After it creates, features are extracted and stored in food category. Next time the text that contain something of that feature will allot that particular category. But the thing is how to automatically first generate category for the text?


title: Is it OK to increase validation checks and decrease min gradient while training neural network?

text: My input vector is a 130*85 matrix and my target vector is 130*26 matrix. I am using the below parameters for training the network with 60 hidden nodes.                      net.trainParam.max_fail = 50;            net.trainParam.min_grad=1e-10;            net.trainParam.show=10;            net.trainParam.lr=0.01;            net.trainParam.epochs=1000;            net.trainParam.goal=0.001;`      As you can see I have set the max_fail to 50 and min_grade to 1e-10. While the default values will be 6 and 1e-5 respectively. But with the default values, the training stops early with out reaching the performance goal. By setting these parameters,the training is stopped when the performance goal is reached. So. Is it OK to change these parameters?


title: Random Forest - How to handle overfitting

text: I have a computer science background but am trying to teach myself data science by solving problems on the internet. I have been working on this problem for the last couple of weeks (approx 900 rows and 10 features). I was initially using logistic regression but now I have switched to random forests. When I run my random forest model on my training data I get really high values for auc (> 99%). However when I run the same model on the test data the results are not so good (Accuracy of approx 77%). This leads me to believe that I am over fitting the test data. What are the best practices regarding preventing over fitting in random forests? I am using r and rstudio as my development environment. I am using the `randomForest` package and have accepted defaults for all parameters


title: Expectation Maximization (EM) Method - All Constructs or One Construct at a time?

text: I want to ask whether I can run Expectation Maximization (EM) method in SPSS to replace ALL missing values of ALL constructs at one run, or I have to do for each construct separately? Thank you