title: Drawing a conclusion from data without a correlation to other supporting data

text: I am looking for the word/phrase that is used to say that a (potential) erroneous conclusion has been drawn from data without having a correlation to other supporting data points. For example, if someone has the dataset of the number of red cars sold in the US per month as `Jan: 10, Feb: 10, Mar: 9, Apr: 8, May: 7, etc.` and draws the conclusion of `"Red cars are going out of style"` without correlating this data to total car sales for the same time period, or the number of other color cars sold, etc. This conclusion might be correct if there was data showing that car sales as a whole have not declined, but this conclusion could be false if car sales as a whole have declined, etc, etc. Is this called "Sampling Bias" or something else?


title: Neighborhood in Gaussian graphical models

text: A gaussian random vector $X$ can be represented using a graph where two nodes $a$ and $b$ are connected $\Leftrightarrow X_a$ is dependent on $X_b$ given all the remaining random variables. I have two questions - 1) It is stated that the neighborhood of $a$ is the smallest subset such that given this subset, $X_a$ is conditionally independent of all the remaining variables. Is this easy to prove or does the proof depend on proving a series of theorems? I'm wondering whether to find and read the proof or just assume this fact. 2) It is generally assumed that the neighborhood of any node is small in size i.e. a node is not connected to all other nodes. Why is this a reasonable assumption in practice? Thanks!


title: If the distribution is bimodal, does p-value mean anything?

text: P-value is defined the probability of obtaining a test-statistic at least as extreme as what is observed, assuming null-hypothesis is true. In other words, $$P( X \ge t | H_0 )$$ But what if the test-statistic is bimodal in distribution? does p-value mean anything in this context? For example, I am going to simulate some bimodal data in R:               set.seed(0)     # Generate bi-modal distribution     bimodal <- c(rnorm(n=100,mean=25,sd=3),rnorm(n=100,mean=100,sd=5))      hist(bimodal, breaks=100)      ![enter image description here](http://i.stack.imgur.com/spCCE.png) And let's assume we observe a test statistic value of 60. And here we know from the picture _this value is very unlikely_. So ideally, I would want a statistic procedure that I use(say, pvalue) to reveal this. But if we compute for the p-value as defined, we get a pretty high-p value               observed <- 60          # Get Pvalue     sum(bimodal[bimodal >= 60])/sum(bimodal)     [1] 0.7991993      If I did not know the distribution, I would conclude that what I observed is simply by random chance. But we know this is not true. I guess the question I have is Why, when computing p-value, do we compute the probability for the values "at least as extreme as" the observed? And if I encounter a situation like the one I simulated above, what is the alternative solution?


title: Markov chain and process

text: This is not homework. Just practicing for an upcoming exam. Question is taken from a web pdf : http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter11.pdf Consider the following process. We have two coins, one of which is fair, and the other of which has heads on both sides. We give these two coins to our friend, who chooses one of them at random (each with probability 1/2). During the rest of the process, she uses only the coin that she chose. She now proceeds to toss the coin many times, reporting the results. We consider this process to consist solely of what she reports to us. (a) Given that she reports a head on the nth toss, what is the probability that a head is thrown on the (n + 1)st toss? (b) Consider this process as having two states, heads and tails. By computing the other three transition probabilities analogous to the one in part (a), write down a \transition matrix" for this process. (c) Now assume that the process is in state "heads" on both the (n - 1)st and the nth toss. Find the probability that a head comes up on the (n + 1)st toss. (d) Is this process a Markov chain? Solution is given here: Chapter 11 question 19 :http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Answersodd-10-14-08.pdf Also, what is the best approach in assigning state i a state? "state 1=..."? Also, what is the difference between the Markov process and Markov chain?


title: Is the beta-binomial distribution a conjugate prior for some distribution?

text: The beta distribution is conjugate prior for the binomial distribution. Is the beta-binomial distribution a conjugate prior for some sampling distribution?


title: Effects in panel models "individual", "time" or "twoways"

text: Panel estimators such the one implemented in the R package `plm` allow to estimate `"individual"`, `"time"` or `"twoways"` effects. See http://cran.r-project.org/web/packages/plm/vignettes/plm.pdf (page 5). When do I use which of the three possible specifications?


title: Examples for teaching: Sometimes we CAN infer causality.

text: There have been threads here before which posted links of the media attributing causality to correlational studies, and links to those studies have been posted. It seems as if we are always focusing on showing students how to identify correlations and not to confuse it with causality. We are missing something important in the process though: Some students assume afterwards that causality can never be inferred, and absolutely everything is a correlation! Does anyone have any links to studies in the media which actually DO show causal links? Student understanding of the difference between correlation/causation in studies reported by journalists, can only be inferred when "causal relationship" sometimes is actually the correct answer.


title: Basic standard error issue with short event

text: I'd like to measure the duration of a function call. The function call has one parameter, `n`. If I were to graph an average of the function call's duration, with `n` on the x axis, I would also provide the standard error for each `n`. The function call duration is, however, extremely short. It seems more appropriate to time, say, 100 calls. Is this approach completely orthogonal to the earlier use of the standard error? Does an average of the 100 calls now count as 1 measurement? Does _it_ have a standard error?


title: Percent error or percent difference?

text: 1. Using either percent error or percent difference, I want to compare one of my measured values from a set to the mean of the set. By reading the Wikipedia article on percent difference, it's still not quite clear which I should choose.   2. What if I wanted to compare an instantaneous slope to the slope of a linear regression, do I use percent error or percent difference? If it isn't clear yet, my misunderstanding lies in the comparison of one value from a set to some computed aggregated value from that set. Thanks in advance. -JP


title: How to log transform Z-scores?

text: The data for my variable is in the form of Z-scores only. I'd like to log transform the scores, but I don't know the mean or standard deviation in order to covert to raw scores. Can I assign an arbitrary mean and standard deviation, use that to convert to raw scores, then log transform the raw scores? Or is there some other way to log transform Z-scores? Thanks.


title: Why genes are assumed to follow multivariate normal?

text: I wonder why gene expression data are very frequently modeled by multivariate normal distributions. What is the reason for those strong assumptions that the genes follow multivariate gaussian? Are there any reasons specific for genetics other than the reasons for general gaussian assumptions (mathematical convenience, etc.)?


title: For what broad types of statistics are integral calculus useful?

text: I understand how differential calculus is useful for basic Maximum Likelihood estimation techniques. However, my question is: what broad types of statistics require an understanding of integral calculus?


title: Testing simultaneous and lagged effects in longitudinal mixed models with time-varying covariates

text: I was recently told that it was not possible to incorporate time-varying covariates in longitudinal mixed models without introducing a time lag for these covariates. Can you confirm / deny this? Do you have any references on this situation ? I propose a simple situation to clarify. Suppose that I have repeated measures (say over 30 occasions) of quantitative variables (y, x1, x2, x3) in 40 subjects. Each variable is measured 30 times in each subject by a questionnaire. Here the final data would be 4 800 observations (4 variables X 30 occasions X 40 subjects) nested in 40 subjects. I would like to test separately (not for model comparison) for :   * simultaneous (synchronous) effects : the influence of x1, x2, and x3 at time t on y at time t.   * lagged effects : the influence of x1, x2, and x3 at time t-1 on y at time t. I hope everything is clear (I'm not a native English speaker !). For instance, in R lmer{lme4}, the formula with lagged-effects is :               lmer(y ~ lag1.x1 + lag1.x2 + lag1.x3 + (1|subject))      where `y` is the dependent variable at time t, `lag1.x1` is the lagged independent variable x1 at the individual level, etc. For simultaneous effects, the formula is :               lmer(y ~ x1 + x2 + x3 + (1|subject))      Everything is running well and it gives me interesting results. But is it correct to specify a lmer model with synchronous time-varying covariates or have I missed something ? _**Edit:_** Moreover, is it possible to test both simultaneous and lagged effects at the same time ? , For instance :               lmer(y ~ x1 + x2 + x3 + lag1.x1 + lag1.x2 + lag1.x3 + (1|subject))      Theoretically, it makes sense to test competition between concurrent vs. lagged effects. But is it possible with `lmer{lme4}` in R, for example ?


title: Derivation of distance in two-step clustering

text: I am working with the two step cluster process in SPSS Modeler (Clementine) and trying to get a sense for the distance function used. It is a log- likelihood function (as stated in docs) but I am unsure for even the continuous variables (the function handles continuous and nominal variable) how this is a log likelihood (it is missing most of the elements of a Gaussian). Below is a screen shot of the documentation describing the distance formula. Has anyone seen the derivation of this distance function? ![enter image description here](http://i.stack.imgur.com/LlSix.jpg)


title: How to test difference between two groups of data in one line

text: Usually, we can test difference of two groups of data by adding the dummy variable. But my question is, if these two groups of data in the same line, but they concentrate at the different part of this line. Is there any way we can test difference? For example,in group one, our data set is $(x_{1i},y_{1i})$, in group two, our data set is $(x_{2j},y_{2j})$. Both of the groups can regress to the same line $y=a+bx$,where $a$ and $b$ are constant. But in group one $x_{1i}\in(1,3)$, in group two $x_{2j}\in(6,7)$. Thus if we plot group one and group two together, we can see they centralize to different part of line $y=a+bx$. Apparently, dummy variable does not work well here. So is there test can show the difference? Thanks.


Can I use a chi-squared test to compare two empirical distributions? Signs point to "yes." R's chisq.test allows two vectors x and y. It says: > cases with missing values are removed, the objects are coerced to factors, > and the contingency table is computed from these. Then Pearson's chi-squared > test is performed of the null hypothesis that the joint distribution of the > cell counts in a 2-dimensional contingency table is the product of the row > and column marginals. See also: same question about KS test.


title: Feature selection for classification, controlling for sub-population

text: I have a bunch of points that belong to one of population P1, P2, ... Pn AND to class A or B. Within each population I'll be doing classification between A and B, and I want to select features that discriminates the best between A and B. Now, my features are also correlated with population membership, but I don't care about that, I only want to know how well the measures discriminate within a population. So, is there a "proper" way to measure that? Within a population, I can just check how much the metric is correlated with class membership (biserial correlation), but if I get the correlation on the whole population, Simpson's paradox gets in the way and I may end up measuring the wrong thing. I could just calculate correlations for each population, and average those or something, but that sounds ugly to me (this is for proramming real-world applications, not for publishing in a paper, so I don't mind if my method isn't _that_ rigorous as long as the result is OK). So, what's the proper way to control for population membership?


title: When do uniformly most powerful (UMP) tests fail to exist and how to prove it?

text: Is there a simple example for a pair of hypotheses for which a uniformly most powerful (UMP) test does not exist? How does one go about proving this? Is it true that for tests with double-sided hypotheses, a UMP test does not exist? Is it true that for multivariate parameters, a UMP tests does not exist?


title: Distance Variable with unequal intervals

text: Suppose that $x$ takes the following values: 0 feet from start, 10 feet from start, and 15 feet from start. Let the outcome variable be $y$. How would you interpret the regression coefficient of $x$ ($\beta_1$) in a linear regression? Would it be: a unit increase in $x$ is associated with an expected increase of $\beta_1$ in $y$? It seems that this would be a problem because a unit could be an increase from $x=0$ to $x=10$ or an increase from $x=10$ to $x=15$.


title: Normalizing constant in Bayes theorem

text: Folks, pardon my noobness but I have not touch maths for some time and need a refresher into statistics: I read that in Bayes rule, the denominator of `P(params|data) = L(params)P(params)/P(data)` is called a normalizing constant. Can anyone explain why does integrating f(data|params)f(params) as in the case of the denominator, does not depend on the parameters?


title: What is v-fold cross validation?

text: What is v-fold cross validation in relation to k-fold cross validation? Also is there a more common way in which v-fold cross validation is referenced? I'm struggling to find resources on this method.


title: Time series autocorrelation

text: The definition for the $k$-th lag auto correlation is $Cov(y_t,y_{t-k})/Var(y_t)$.   My question is why should not it be $Cov(y_t,y_{t-k})/[Var(y_t)\cdot Var(y_{t-k})]^{0.5}$.   In another words, what is it different from the correlation coefficient between $y_t$ and $y_{t-k}$?


title: Linear mixed effect model vs. Ordered Probit vs. Ordered Logit with ordinal response

text: I have a set of data with an ordinal response ranging from `1-5` (worst to best) and a categorical predictor with five unordered levels. The experiment is a language experiment whereby subjects are asked to rate different sentence types. In the literature it seems that people fit `lmer()` most of the time by scaling the ordinal response per subject (i.e. taking the mean response and standard deviation of that subject and dividing the every single response by it.). On the other hand, it seems a reasonable approach to use ordered linear mixed effects probit models or ordered linear mixed effects logit models as implemented in e.g. `MCMCglmm()` or `clmm()`. (1) I find myself asking what the best method would be. When is one method to be prefered over the other in an experiment as I have sketched above? (2) Why do people prefer fitting an `lmer()` model with a scaled response? What is the advantage? (3) If it is advisable to use an ordered probit or logit model, how do I decide between an ordered probit or ordered logit model? Thanks for any help!


title: If a statistics doesn't reveal a significance do I have to calculate power for it?

text: Following the design and data described in this question, I did a simple one- way within-subjects repeated-measures (RM) ANOVA and found some significant p-values. I then applied non-orthogonal post-hoc Tukey's HSD tests, and when I got significant results I applied Holm-Bonferroni (1979) correction. Whenever some p-values survived the FWER correction, I calculated 95% CIs and mean for the associated pairwise comparisons. My question is: If I don't observe a significant result at any of the above steps, do I have to carry out a power analysis for the RM ANOVA, apply Tukey's HSD test or Holm-Bonferroni adjustments, or do I simply report results from the RM ANOVA without doing the power analysis? The problem is that I'm starting to immerse in biostatistics only after my experiments, and unfortunately I didn't run a power analysis beforehand.


title: How to compute goodness of fit after applying logarithmic transformation?

text: I want to estimate the parameters of a function of general form $y = a \cdot x^b$. I applied a log-log transformation to obtain a linear function of the form $\log y = \log a + b \times \log x$. I have fitted the linear model in MATLAB. MATLAB computes the goodness of a fit in terms of sum of squares error (SSE) and (adjusted) R bar. I want to report these numbers but they are probably meaningless given that in fact $\log a$ and not $a$ has been estimated. How do I fix this?


title: Is this a valid use case for Euclidean distance?

text: I have a set of points which is a count of links that users have clicked on :                        link1,link2,link3,link4     user1,   0    ,2    ,3    ,2          user2,   1    ,4    ,0    ,5          I am not sure which distance function I should use to calculate the distance between users. I think Euclidean distance is a good fit (reading http://en.wikipedia.org/wiki/Euclidean_distance) ?


title: How does I bias my standard errors estimates when I forget to correct for autocorrelation of errors?

text: Suppose that I have a linear model with autocorrelated errors. Is there any results telling me that if I assume iid errors I overestimate or underestimate my standard errors ?


title: Likert scale coding

text: I have a question on coding a Likert scale in SPSS. Assuming we have a 7-point Likert scale, a "well-being" score may rank from 1 through 7. While doing some reading on multi-variate analysis, I've come across a reference that suggests a "deviation score". That is, it's ranging from [-3, -2, -1, 0, 1, 2, 3]. For later analysis of survey results, which method ("well-being score" vs. "deviation score") will serve me best? Or are these "apples and oranges"? Thanks, EEH


title: Linear regression with linear scale as explanatory variable

text: I have a dataset containing values representing monthly sales in hypothetical e-commerce website. First columns contains order of particular month (eplanatory variable) and the second column (dependent variable) is value of sales in this month.               month_order value         1   2591         2   2262         3   2531         4   2379         5   2106         6   2500         7   2345         8   2200         9   2385         10  2991         11  2219         12  2600         13  3000      Equation of linear regression for this data set is:               f(x) = 27,42x + 2278      **Question:** Because I would like to make a very simple prediction of sales value in next month and I have no another more suitable eplanatory variable, may I use month_order and predict linear trend? In my case it would be:               f(14) = 27,42 * 15 + 2278 = 2689      is this correct in a theory of a linear regression and can I count it in this way?


title: Interpreting exp(B) in multinomial logistic regression

text: This is somewhat of a beginner's question, but how does one interpret an exp(B) result of 6.012 in a multinomial logistic regression model? 1) is it 6.012-1.0 = 5.012 = 5012% increase in risk? or 2) 6.012/(1+6.012) = 0.857 = 85.7% increase in risk? In case both alternatives are incorrect, can someone please mention the correct way? Ive searched many resources on the internet and I get to these two alternatives, and I am not entirely sure which one is correct.


title: Is there a reference that legitimises the use of the unpooled z-test to compare two proportions?

text: The z-test to compare two proportions is $\newcommand{\p}{\hat{p}}\newcommand{\v}{\mathrm{Var}} z=\frac{\p_1-\p_2}{\sqrt{\v(\p_1-\p_2)}}$. Usually it is defined that $$\v(\p_1-\p_2)=\p(1-\hat{p})(1/n_1+1/n_2),$$ where $$\p=\frac{n_1 \p_1+n_2 \p_2}{n_1+n_2}.$$ Is there any written reference that legitimizes me instead to use the unpooled variance, that is $$\v(\p_1-\p_2)=\frac{\p_1(1-\p_1)}{n_1}+\frac{\p_2(1-\p_2)}{n_2}?$$


title: Geodesic distance and mean

text: My data-set consists of points in globe. Suppose a User visits locations $l_1,l_2,\dots l_n$ (each location in $(lat, long)$ in the city with probability $p_1,p_2,\dots,p_n$ and I want to calculate the expected location visited by the user. How should I go about calculating that? If I use geodesic distance between points then do we have a formula or recipe to calculate geodesic expectation? What possibly can go wrong if I convert the individual $(lat,long)$ values to cartesian co-ordinates and calculate their expectation?


title: Post Hoc test for between subject factor in a repeated measures ANOVA in R

text: I have done an experiment of how different animal species affect nutrient fluxes in sediments. I had a number of experimental units (sediment boxes) to which different animal species were added. I then measured nutrient fluxes in these units. Each unit was measured five times. I have used the ezANOVA package for R to do a repeated measures ANOVA with two factors (between subjects factor “Species” and within subject factor “Time”). Here's an example of my input:               > ezANOVA(data=NoP_3_5,dv=.(AcPO4),wid=.(Subject),within=.(Time),between=.(Species),return_aov=T)     $ANOVA         Effect DFn DFd           F            p p<.05        ges     2      Species   2   6 10.60384830 1.072453e-02     * 0.50875722     3         Time   2  12 27.88590570 3.081718e-05     * 0.76667542     4 Species:Time   4  12  0.08093179 9.867291e-01       0.01871588          $`Mauchly's Test for Sphericity`             Effect         W           p p<.05     3         Time 0.1439855 0.007866779     *     4 Species:Time 0.1439855 0.007866779     *          $`Sphericity Corrections`             Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05     3         Time 0.5387889 0.001349263         * 0.5630403 0.001103048         *     4 Species:Time 0.5387889 0.933492515           0.5630403 0.939138882                    $aov          Call:     aov(formula = formula(aov_formula), data = data)          Grand Mean: -0.004276663           Stratum 1: Subject          Terms:                          Species    Residuals     Sum of Squares  0.0004410742 0.0001247870     Deg. of Freedom            2            6          Residual standard error: 0.004560464      Estimated effects may be unbalanced          Stratum 2: Subject:Time          Terms:                             Time Species:Time    Residuals     Sum of Squares  0.0013994206 0.0000081229 0.0003011028     Deg. of Freedom            2            4           12          Residual standard error: 0.005009181      Estimated effects may be unbalanced      I want to use the $aov output to do a Tukey HSD post hoc test of the between factor ("Species"). However, it does not work:               mod <- ezANOVA(data=NoP_3_5,dv=.(AcPO4),wid=.(Subject),within=.(Time),between=.(Species),return_aov=T)         > TukeyHSD(mod$aov)     Error in UseMethod("TukeyHSD") :        no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"      I have also tried to use the aov-command directly:               > mod2 <- aov(AcPO4~(Species)+Error(Subject/Time)+(Species),data=NoP_3_5)     > TukeyHSD(mod2)     Error in UseMethod("TukeyHSD") :        no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"      I have tried to use lme() to specify the anova-model with the aim to produce something that TukeyHSD() will accept without success. Any help would be greatly appreciated!


title: Why must an estimator be independent from the parameter?

text: ![enter image description here](http://i.stack.imgur.com/IBHdy.png) This is an excerpt from "Modern mathematical statistics with applications" by Devore et al. What puzzles me is that the estimator cannot help being dependent on $\theta$, since the sample depends on the parameter.


title: What tasks should I do when performing a model validation?

text: I have a mixed model built using Data A below, and now I want to validate the model using Data B. What should I do to actually "validate" the model?               model <- lme(Y~1+X1+X3, random=~1|School, method="ML", data=A)      I know how to use "predict(model,B) - B$Y" to get residuals, but what else should I also perform in this validation process? (I will not apply the cross-validation technique since my data are hierarchical structured.)


title: Chi square based measure for fairness?

text: I have some project categories $(PC_1, ... ,PC_n)$ and some institutions $(I_1, ..., I_n)$ and would like to determine how fairly projects are assigned to institutions and their project categories by some external sources. I was thinking about using some chi square based measure for this. Basically I have a two dimensional contingency table where one dimension corresponds to the project categories and the other to the institutions. Each cell of the contingency table counts the number of projects that each institution ‘runs’ for a particular project category. I could determine $(O – E)\cdot|O-E|/E$ for each cell. Here:   * $O$ is the observer number of projects   * $E$ is the expected number of projects IMHO this would give me a value for each cell which indicates how much the observed value differs from the expected value and thus how fairly the resources were assigned for an institution and project category. In other words a negative value indicates ‘unfairness’, a positive value indicates ‘favouritism’ and a value close to zero indicates ‘fairness’. What do you think?


title: Can anyone explain why I have obtained an anti-predictive Support Vector Machine?

text: I'm playing with support vector machines (SVM) using the e1071::svm() function in R, and I encountered a scenario where I asked it for a leave-one-out cross- validated classification of a 2-category response and obtained a total accuracy of 38% (35/90), which, given 90 samples, ends up with a 95% confidence interval that is **below** chance. Should I consider this a fluke, and if not, how is it possible for an SVM to become anti-predictive? In case it matters, I used default values for the cost and gamma parameters, and the data predicting the response was a 8192 item vector representing 500 milliseconds of electroencephalogram data collected across 64 electrodes.


title: Compairing 2 ideas - multiple dependent variables

text: I am trying to compare 2 ideas to see which one is better/more effective - I am looking to see if a significant difference. (2 conditions) We used a questionnaire to gather feedback from hundreds of people. (Questions include: how likely would you buy it, how likely would you tell your friends and family, I want to learn more about this product...etc. An swers are multiple choice: strongly agree,somewhat agree, neutral, somewhat disagree, strongly disagree & no advantage, some advantage, slight advantage, a great advantage) **WHAT SHOULD I DO? AND IS THIS THE RIGHT WAY OF DOING IT? (using SPSS)**   1. t-test: by combining the data (would I just give each answer a value like 1-4, add values up for total score for each participant)   2. Run cross tabs and then do something?   3. Other suggestions? THANKS! :)


title: Why only three partitions? (training, validation, test)

text: When you are trying to fit models to a large dataset, the common advice is to partition the data into three parts: the training, validation, and test dataset. This is because the models usually have three "levels" of parameters: the first "parameter" is the model class (e.g. SVM, neural network, random forest), the second set of parameters are the "regularization" parameters or "hyperparameters" (e.g. lasso penalty coefficient, choice of kernel, neural network structure) and the third set are what are usually considered the "parameters" (e.g. coefficients for the covariates.) Given a model class and a choice of hyperparameters, one selects the parameters by choosing the parameters which minimize error on the training set. Given a model class, one tunes the hyperparameters by minimizing error on the validation set. One selects the model class by performance on the test set. But why not more partitions? Often one can split the hyperparameters into two groups, and use a "validation 1" to fit the first and "validation 2" to fit the second. Or one could even treat the size of the training data/validation data split as a hyperparameter to be tuned. Is this already a common practice in some applications? Is there any theoretical work on the optimal partitioning of data?


title: Two envelope problem revisited

text: I was thinking of this problem. http://en.wikipedia.org/wiki/Two_envelopes_problem I believe the solution and I think I understand it, but if I take the following approach I'm completely confused. Problem 1: > I will offer you the following game. You pay me \$10 and I will flip a fair > coin. Heads I give you \$5 and Tails I give you \$20. > > The expectation is \$12.5 so you will always play the game. Problem 2: > I will give you an envelope with \$10, the envelope is open and you can > check. I then show you another envelope, closed this time and tell you: This > envelope either has \$5 or $20 in it with equal probability. Do you want to > swap? > > I feel this is exactly the same as problem 1, you forgo \$10 for a \$5 or a > \$20, so again you will always switch. Problem 3: > I do the same as above but close the envelopes. So you don't know there are > $10 but some amount X. I tell you the other envelope has double or half. Now > if you follow the same logic you want to switch. This is the envelope > paradox. What changed when I closed the envelope?? EDIT: Some have argued that problem 3 is not the envelope problem and I'm going to try and provide below why I think it is by analysing how each views the game. Also, it gives a better set up for the game. Providing some clarification for Problem 3: From the perspective of the person organising the game: > I hold 2 envelopes. In one I put \$10 close it and give it to the player. I > then tell him, I have one more envelope that has either double or half the > amount of the envelope I just gave you. Do you want to switch? I then > proceed to flip a fair coin and Heads I put \$5 in and Tails I put \$20. And > hand him the envelope. I then ask him. The envelope you just gave me has > twice or half the amount of the envelope you are holding. Do you want to > switch? From the perspective of the player: > I am given an envelope and told there is another envelope that has double or > half the amount with equal probability. Do I want to switch. I think sure I > have $X$, hence $\frac{1}{2}(\frac{1}{2}X + 2X) > X$ so I want to switch. I > get the envelope and all of a sudden I am facing the exact same situation. I > want to switch again as the other envelope has either double or half the > amount.


title: What does averaging do for the noise level of my signal

text: We are measuring data with a high sample rate (20 kHz) and calculated a big standard error due to our system setup. Currently we are only interested in slow signals (in the order of Hz's) is it valid to use averaging (once per 20.000 samples) and thus lower our standard error with an order 20.000?


title: How to draw a balanced random sample given prediction confidence and model performance?

text: **INTRO** I have trained a SVM model based on 300 training cases in order to build a filter that should help me extract a bigger and evenly balanced training set, that is to be validated by human judges. **DATA** My training set is:               288 negative cases     12 positive cases      My model performance is:               True positives = 6     False positives = 0     True negatives = 288     False negative = 6      I have used this model on a 440000 cases data set and I am given confidence values for zeroes and ones (confidence.1. = 95% -> confidence.0. = 5% and confidence.1. = 50% -> confidence.0. = 50%) When the model is applied to the 440000 cased it predicts 1499 of the cases to be positive. **PROBLEM** As explained I am about to take a random sample of 3000 cases that will hopefully be balanced in order to train a better model. Thus my formal problem becomes: _"How can one extract a random sample of 3000 and approximate the balance of the dataset given the information I have stated?"_ **OWN THOUGHTS** I am fairly good with R but I simply do not have the knowledge or the skills to figure this out. I dont even know if it is possible. The only thing I do know is that it seems a bit doll to just put a finger in the air and guess on which cases to extract. The easy way would be to take the top 3000 sorted on confidence.1. which is nothing more that a qualified guess. Happy new year Kasper


title: Leave-one-out cross-validation for MLR

text: I am doing cross-validation using the leave-one-out cross-validation method (total 10 runs). I have predicted $\hat{y}$ and observed $y$ from all the runs. I have applied the following equation to calculate $R^2$. Can anyone please confirm whether I am doing right? I am confused because the $R^2$ value appeared negative. \begin{align} SSE&=∑(y-\hat{y})^2 \\ SST&=∑(y–\bar{y})^2 \\ R^2&=1-(SSE/SST) \end{align} (N.b., $\bar{y}=$ is the average of $y$.)


title: What are the most popular domain adaptation methods (for transfer learning)?

text: I understand supervised and unsupervised learning well, and would be able to identify some 'basic' examples of, for example, supervised classifcation as:   * SVMs   * Random Forests   * Logistic Regression These are key works in the field which have lots of code and publications available. I am now starting to look at domain adaptation in supervised learning, where the distribution over the data at learning and testing time are known to be different. Despite reading some of the literature, I haven't spotted any similar 'basic' methods which come up time and time again. In contrast, there seem to be a wild array of completely different methods for achieving the goal, many of which have only been in the literature for a few years. Are there such key, established methods for domain adaptation? What are the most popular methods currently used?


title: How the 'NA' values are treated in glm in R

text: I have a data table T1, that contains nearly a thousand variables (V1) and around 200 million data points. The data is sparse and most of the entries are NA. Each datapoints have a unique id and date pair to distinguish from another. I have another table T2, which contains a separate set of variables (V2). This table also has id and date pair that uniquely identify entries in T2. We suspect that the data in T1 can be used to predict values of variables in T2. To prove this, I thought to apply 'glm' model in R and check if we can really find some variable in T2 that is dependent on variables in T1. For each variable in T2, I started pulling out all data in T1 having same id and date pair which resulted in much smaller ~50K data points for some of test variables. The problems that I am facing now with application of glm is as follows.   1. In some cases, it shows me an error 'fit not found' and warning 'glm.fit: algorithm did not converge '. I am not sure why is it shown?    2. How the NAs are treated in glm? Does it remove all records involving 'NA' first and then perform fitting?   3. Is it a good strategy to remove all NAs first and then call 'glm'. I fear that this may reduce the datapoints significantly as most of them are NAs.    4. Which method is used to calculate the coefficients. I couldnot find any website or paper or book that discuss how the output is calculated. I tested glm with and without 'NAs' and found difft answers which points that NAs are considered while fitting the data: Example 1:               > tmpData       x1 x2 x3        Y     1  1  1  1        3     2  1  0  4        5     3  1  2  3        6     4  0  3  1        4          Call:  glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse = "+"))),          na.action = na.exclude)          Coefficients:                           (Intercept)  as.numeric(unlist(tmpData["x1"]))  as.numeric(unlist(tmpData["x2"]))                               5.551e-16                          1.000e+00                          1.000e+00       as.numeric(unlist(tmpData["x3"]))                               1.000e+00            Degrees of Freedom: 3 Total (i.e. Null);  0 Residual     Null Deviance:      5      Residual Deviance: 9.861e-31    AIC: -260.6       Example 2:               'x1'    'x2'    'x3'    'Y'     100000  1   NA  2     1   1   1   3     1   NA  -1124   2     1   0   4   5     1   2   3   6     0   3   1   4                    Coefficients:                           (Intercept)  as.numeric(unlist(tmpData["x1"]))  as.numeric(unlist(tmpData["x2"]))  as.numeric(unlist(tmpData["x3"]))                              -2.3749044                         -0.0000625                          0.6249899                          1.8749937            Degrees of Freedom: 5 Total (i.e. Null);  2 Residual     Null Deviance:      13.33      Residual Deviance: 1.875    AIC: 20.05


title: Is there a package in R that calculates the Durbin-Watson test assuming a weighted regression?

text: There are several packages that can apply the Durbin-Watson test for serial correlation. However, I do not see a package that supports the calculation in the case that once has a GLS weighted regression. For example, CRAN package lmtest notes in their changelog that they explicitly do not support weighted regressions (yet). Before the recent release, lmtest would not throw an error when passing weighted regressions. My concern is that the other dwtest packages may also not be explicitly dealing with this scenario.


title: I don't know how should I calculate meaningful error bars in this special case

text: I have **6** masses as input data which they are expressed as m+/- delta m. Via a complicated processes I calculate **decay width** as a final result. It takes a long time for the code to compute decay width for example for **m** values. Let's call it **(decay width)_central**. Let's call **m + delta m= m_max** and **m- delta m =m_min**. (For 6 masses: For example m_max means I substitute all 6 masses with their maximum value). Then decay width for m_max is **(decay width)_max** and decay width for m_min is **(decay width)_min**. Can I have an estimate for error bar with only these 3 values of decay widths? (It is very time consuming to find more than these 3 results.) Or may be you suggest a better way? Thanks a lot. My problem is how to calculate error bars for decay width. Errors originate from 6 experimental data which I used as inputs. The question is that how should I find error bars for decay width? I can't use the usual methods such as deviation for errors, because it is so time consuming to calculate decay width for the masses in experimental range and then find mean value. Can I find an estimate for error bars only with 3 values of decay width? The decay width which I found by substituting all the masses with m + delta m, The one by substituting all the masses with m - delta m and the one which I found by substituting the mean values for m?


title: Aligning sample data

text: Given a data matrix like this:                       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]      [1,]  9.520 11.137 16.576 18.225 20.576 25.861     NA      [2,]  9.005  9.491 11.106 16.530 18.184 20.495 25.773      [3,]  9.437 11.050 20.393 25.711     NA     NA     NA      [4,]  9.442 11.058 20.411 25.711     NA     NA     NA      [5,]  9.431 11.045 20.421 25.707     NA     NA     NA      [6,]  9.461 11.052 20.319 25.657     NA     NA     NA      [7,]  9.245 10.819 20.253 25.628     NA     NA     NA      [8,]  9.229 10.801 20.216 25.594     NA     NA     NA      [9,]  9.234 10.805 20.258 25.619     NA     NA     NA     [10,]  9.241 10.814 20.264 25.626     NA     NA     NA     [11,]  9.248 10.819 20.281 25.649     NA     NA     NA     [12,]  9.231 10.800 20.219 25.567     NA     NA     NA      How do I get a data frame like this?:                S     p1     p2     p3     p4     p5     p6     p7      1     NA  9.520 11.137 16.576 18.225 20.576 25.861           2  9.005  9.491 11.106 16.530 18.184 20.495 25.773      3     NA  9.437 11.050     NA     NA 20.393 25.711      4     NA  9.442 11.058     NA     NA 20.411 25.711      5     NA  9.431 11.045     NA     NA 20.421 25.707      6     NA  9.461 11.052     NA     NA 20.319 25.657      7     NA  9.245 10.819     NA     NA 20.253 25.628      8     NA  9.229 10.801     NA     NA 20.216 25.594      9     NA  9.234 10.805     NA     NA 20.258 25.619     10     NA  9.241 10.814     NA     NA 20.264 25.626      11     NA  9.248 10.819     NA     NA 20.281 25.649      12     NA  9.231 10.800     NA     NA 20.219 25.567       In the matrix, each row represent the timepoints from 1 sample. The matrix is balanced with an arbitrary number of NAs as needed. In the dataframe, each column represents all the timepoints that are functionally equivalent. Timepoints in a single sample cannot be functionally equivalent. I think I should be able to solve this with some form of kmeans clustering that restricts each cluster from having more than one member from each sample. Any ideas?


title: Bank loan model (spm?)

text: I am trying to fit a model investigating the amount of the loan (or a transformation of it) as a function of the the variables `income`, `gender`, `customer` and `age`. Fitting a standard linear regression gives poor results (R-squared of only 0.27). Besides, the relationship between income and loan seems quite weird. It has a sort of left-rotated Z-shape Would be appropriate to use a semiparametric model here? And so yes, how would you model it? I tried this:               model <- spm(loan~f(income, basis="trunc.poly")+age+factor(gender)+factor(customer))      But this gave poor results (no significant variables). Anyone has suggestions? Many thanks in advance.


title: Is my Theory on the Law of Truly Large Numbers and Ancestors Valid?

text: I am an anthropology student. I am researching the inter-relatedness of most if not all humans as part of my studies. I personally have an interest in statistics and probability theory and I have used it in regards to this topic. Before I explain my theory, I received the following from a computer scientist that I wish to share to better give a background idea of what I’m trying to put forward:   If there were random intermixing, then we would each have ~1 million ancestors living in 1500 AD, out of a world population of ~500 million. So the fractional overlap between two people would be about 1/500th. But the probability that two people share at least one common ancestor would be essentially 100%. Basically, you are choosing a random number between 1 and 500 a million times and you're asking whether you ever choose number 500. In a million trials, we expect this to happen 2000 times. So that it happens at least once is guaranteed. If we get rid of the random intermixing, the fractional overlap will drop to much less than 1/500th. But I suspect that the probability of at least one overlap will remain very high. If the population in 1500 was 500 million, and it is 6 billion today (12x larger). If the average generation length is 30 years, there are 17 generations in 500 years. So the average number of surviving children per mother is $\exp((\log 12)/17) = 1.157$ Since a child has two parent, the average number of surviving children per person is 2 * 1.157 = 2.315 So this is the average growth rate per generation for the descendants of a person in 1500. $2.315^{17}$ = 1.575 million. So an average person in 1500 has about 1.5 million offspring alive today. Sampling from the whole world, the probability that a random person from 1500 is an ancestor of a random person in 2000 would be 1.5 million / 6 billion = 0.025%. If you were only considering people in a region like Europe, it would probably be something like 1.4 million / 700 million = 0.2%. As I have said above, this is from a computer scientist. Next, is my own theory. From a probabilistic standpoint, due to the Law of Large Numbers, me, you or *almost anyone alive today had an ancestor in like say feudal Japan in 1500 CE. I chose Japan because of its remoteness yet it still had contact with the rest of the world at that time and I chose the year 1500 because of each person had in theory about a million ancestors- lower than that because of inbreeding but still a sizable number nonetheless. I have shared this theory with several of my colleagues and they are intrigued by it yet I have not yet shared it with an actual statistician. I know it seems unlikely, but from a strictly statistical standpoint, what do you think of it? *please note, this excludes populations that have not have had contact.


title: Evaluating unbiased errors on the Test set when cross validation scores are close

text: I trained a few different models, (Perceptron, Stochastic Gradient Descent and Naive Bayes), each with different parameters. I then scored their accuracy on a cross validation set. The scores on the best parameter Perceptron, SGD and NB models were 93%, 91% and 94% respectively. I didn't expect such similar results and I'm at a bit confused because I feel that the possibility of variance makes choosing the NB as the best model questionable. Am I supposed to test all 3 on the test set and use the model with the best unbiased error? Or is that implicitly cherry picking the best model ?


title: Expected value of certain exponential transformation of standard normal variable

text: This is in reference to the Girsanov theorem however question is general. If $X$ is a standard normal variable $N(0,1)$, why is expectation of $e^{-\mu X - \mu^2/2}$ equal to 1? Shouldn't it be $e^{-\mu^2/2}$?


title: Decision tree for numeric dependent variable?

text: I have data on commute times over a specified route over different days during different conditions. Some of the conditions are categorical (e.g., weather, traffic), and some of them are numeric (e.g., time departing from origin). I'd like to find out which of these conditions most strongly correlate with the shortest commute times on any given day of the week (M-F). Is there a decision tree-like algorithm that would discover the conditions most strongly correlated with the shortest commute times?


title: Comparing two linear regression models

text: I would like to compare two linear regression models which represent degradation rates of a mRNA over time under two different conditions. The data for each model collected independently. Here is the dataset.               Time (hours)    log(Treatment A)     log(treatment B)     0   2.02    1.97     0   2.04    2.06     0   1.93    1.96     2   2.02    1.91     2   2.00    1.95     2   2.07    1.82     4   1.96    1.97     4   2.02    1.99     4   2.02    1.99     6   1.94    1.90     6   1.94    1.97     6   1.86    1.88     8   1.93    1.97     8   2.12    1.99     8   2.06    1.93     12  1.71    1.70     12  1.96    1.73     12  1.71    1.76     24  1.70    1.46     24  1.83    1.41     24  1.62    1.42      These are my models:               Exp1.A.lm<-lm(Exp1$Time~Exp1$(Treatment A))     Exp1.B.lm<-lm(Exp1$Time~Exp1$(Treatment B))                    Call:     lm(formula = Exp1$Time ~ Exp1$(Treatment A))          Residuals:         Min      1Q  Median      3Q     Max      -6.8950 -1.2322  0.2862  1.2494  5.2494           Coefficients:                        Estimate Std. Error t value Pr(>|t|)         (Intercept)           74.68       6.27   11.91 2.94e-10 ***     Exp1$(Treatment A)   -36.14       3.38  -10.69 1.77e-09 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 2.97 on 19 degrees of freedom     Multiple R-squared: 0.8575, Adjusted R-squared:  0.85      F-statistic: 114.3 on 1 and 19 DF,  p-value: 1.772e-09               Call:     lm(formula = Exp1$Time ~ Exp1$(Treatment B))          Residuals:        Min     1Q Median     3Q    Max      -7.861 -3.278 -1.444  3.222 11.972           Coefficients:                           Estimate Std. Error t value Pr(>|t|)         (Intercept)             88.281     16.114   5.478 2.76e-05 ***     Exp1$(Treatment B)  -41.668      8.343  -4.994 8.05e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 5.173 on 19 degrees of freedom     Multiple R-squared: 0.5676, Adjusted R-squared: 0.5449      F-statistic: 24.94 on 1 and 19 DF,  p-value: 8.052e-05      To compare these two models, I used this following code.               anova(Exp1.A.lm,Exp1.B.lm)                    Analysis of Variance Table          Model 1: Exp1$Time ~ Exp1$Exp1$(Treatment A)     Model 2: Exp1$Time ~ Exp1$Exp1$(Treatment B)       Res.Df    RSS Df Sum of Sq F Pr(>F)     1     19 167.60                           2     19 508.48  0   -340.88 My question is why the ANOVA analysis doesn't show an F statistics and a p.val. My apologies if this is a naive question. Based on different slopes, the rate of degradation is different in these two models, but I would like to know how statistically significant this difference is. I hope that this makes sense. Thank you-


title: Nonparametric Topic Models

text: I am confused between what types of problems these three models capture, and their applications:   * Latent Dirichlet Allocation (LDA)   * Dirichlet Processes and Pitman-Yor processes   * Hierarchical Dirichlet Process (HDP) & Hierarhical Pitman-Yor Processes More than formal definitions (which I can find on Wikipedia), I am looking for the **intuition** behind each of them ( how they build on each other). Perhaps, more specifically:   * Is the difference between LDA and HDP that LDA is parameteric (i.e. I need to pre-specify the number of topics) whereas HDP is non-parametric? (and therefore I don't need to know how many topics I have)   * What is the difference between a Dirichlet Processes and Pitman-Yor Processes?    * What is the difference between a non-parametric process (e.g. DP) and a hierarhical non-parametric process (e.g. HDP)?


title: Poisson State Space with AR(1) latent process

text: I have been trying to use sspir R package to estimate the following Poisson model: $Y_{t}\sim Po(\exp(\lambda_{t}));$ such that $\lambda_{t}=X_{t}\beta +\gamma_{t}$ and $\gamma_{t}=\theta\gamma_{t-1}+u_{t},u_{t}\sim NID(0,\sigma ^{2}).$ In State-space form, we have $\lambda_{t}=Z_{t}F_{t},$ $Z_{t}=\left[ X_{t}\: 1\right] ,$ $F_{t}=\left[ \beta^{\prime}\: \gamma_{t}\right] $ Then, $F_{t}=\left( \begin{array} [c]{ll}% 1 & 0\\ 0 & \theta \end{array} \right) F_{t-1}+\left( \begin{array} [c]{l}% 0\\ u_{t}% \end{array} \right) $. Standard initial conditions are also set. I have been trying to use `sspir` package in R to estimate this model, but I was not able to. If I let $\gamma_{t}$ to be a random walk, then it is standard. Any suggestion? I am kind of new with this models so any help is very much appreciated.


title: Gap in plot between forecasted and real data

text: I'm trying to plot forecast and real data on the same plot. But there is always gap between them on the image : ![enter image description here](http://i.stack.imgur.com/ni7FR.png) I use for the forecasting following code: y=boardings[,1] ## Simple Exponential smoothing #predict existing values results=HoltWinters(y,beta=FALSE, gamma=FALSE) #print data print(cbind(y,results$fitted)) #draw the plot with fitted values plot(results) #predict future values results2=forecast.HoltWinters(results,h=12) print(results2)               #draw the plot with prediction     plot.forecast(results2,lwd=1,xlab='time',ylab='log.boardings',main='12-month prediction by Exponential Smoothing')      The data is taken from the timeSeries package with name `boardings`. I would really appreciate any help.


title: When are the asymptotic variance of OLS and 2SLS equal?

text: Assume the model $ \ y = X\beta + u \ $ with $\ W \ $ is a $ \ n\times l \ $ so called matrix of instruments. The following assumptions hold. There is a law of large numbers (LLN) for 1.,2.,3. and 4. such that   1. $\text{plim}_{n\to\infty} \ \left(\frac{X^TX}{n}\right) = m_{X^TX}$ holds where $m_{X^TX}$, a finite, non stochastic matrix with full column rank, exsists.   2. $\text{plim}_{n\to\infty} \ \left(\frac{W^TW}{n}\right) = \text{lim}_{n\to\infty} \left(\frac{{\mathbb E}\left(W^TW\right)}{n}\right)$ where we assume, that the RHS exists and is finite as well as positive definite.   3. $\text{plim}_{n\to\infty} \ \left(\frac{W^TX}{n}\right)$ where we assume, that the limit $W^TX$ exists, is finite and $W^TX$ has full column rank i.e. $\text{rk}\left(W^TX\right)=k$.   4. $\text{plim}_{n\to\infty} \ \left(\frac{W^Tu}{n}\right) = \text{lim}_{n\to\infty} \left(\frac{{\mathbb E}\left(W^Tu\right)}{n}\right) = 0$ where we assume that $\text{lim}_{n\to\infty} \left(\frac{{\mathbb E}\left(W^Tu\right)}{n}\right)$ is equal to $0$. The asymptotic variance for OLS will, under the assumption of homoscedastic errors $u$ i.e. that ${\mathbb E}\left(uu^T|X\right) = \sigma^2_0 I_n$ and that the model is actually correctly specified , be $\text{plim}_{n\to\infty} \text{Var}\left[ n^{\frac{1}{2}}\left(\boldsymbol{\widehat{\beta}}_{\text{KQ}} - \boldsymbol{\beta_0}\right)\big| X\right] = \sigma_0^2 \ m_{X^TX}^{-1}$ The corresponding asymptotic variance for the 2SLS case will look like $\text{plim}_{n\to\infty} \text{Var}\left[ n^{\frac{1}{2}}\left(\boldsymbol{\widetilde{\beta}_{\text{2SLS}}} - \boldsymbol{\beta_0}\right)\big| X\right] = \sigma_0^2 \ \text{plim}_{n\to\infty} \left(\frac{X^TP_W X}{n}\right)^{-1}$ If I now consider the precision matrix rather then the variance matrix and take their difference I'll get $\text{plim}_{n\to\infty} \ \left(\frac{X^TM_WX}{n}\right) = \left(\text{plim}_{n\to\infty} \ \frac{X^TX}{n}\right) - \left(\text{plim}_{n\to\infty} \ \frac{X^TP_WX}{n}\right)$ which is a positive semidefinite matrix. **My question: When are both asymptotic variances of OLS and 2SLS equal** My guess is that we need to set $X^TX = X^T P_W X$. The only possibility for such a thing is when $P_WX=X$. So all the columns of $X$ must be in the image of $P_W$. This means that all the columns of $X$ are actually viable instruments. By that I conclude that there are (at least asymptotically) no endogenous regressor within $X$ if their asymptotic variances are equal. But this seems to good to be true. There are downfalls of using a 2SLS-approach if theres no endogeneity involved. Again, this would mean that if the variances of the 2SLS and OLS estimator are equal asymptotically it does not matter if we choose 2SLS or OLS since both estimators are consistent for $\beta$. Any thoughts about this proble would be very much appreciated Druss


title: Does it exist an Bayesian "version" of the Theil-Sen estimator?

text: From what I read the Theil-Sen estimator seems to be a really cool estimator; robust and easy to understand. **Would it be possible to device (or does it already exists) a Bayesian "version" of the Theil-Sen estimator so that you provide a prior for the slope and get a probability density back instead of just a point estimate? How would such an estimator look/work?**


title: How to reduce autocorrelation in Metropolis algorithm?

text: I've been using a Metropolis/Gibbs sampler combination to generate a joint density for some parameters(it is a hierarchical model, with $y_i\sim Poisson(\lambda_i)$, $\lambda_i\sim Gamma(\alpha,\beta)$). What techniques can I use to lower autocorrelation(it is present in $\alpha$ and $\beta$)? I have been using thinning, but even when I use huge lags(4950, which is reaching the memory limit on my computer to use) there is still significant autocorrelation. Is there something I could do with my step size distribution to help with this? I have been drawing the new values of $\alpha$ and $\beta$ from a normal distribution, with mean the current parameter value, and standard deviation 1. Thanks!


title: How do I make this plot using R?

text: I'm familiar with very basic plotting in R, but I'm not sure how best to create the reasonably-complicated plot described below. I have developed a density estimation method that essentially fits an empirical distribution to multiple univariate samples, under distinct conditions. To be clear, let $\mathbf{x}_i^\mathrm{j}$ be a vector of univariate values for the $i^\mathrm{th}$ sample under condition $j$. My method allows me to approximate the distribution from which the values in $\mathbf{x}_i^\mathrm{j}$ were drawn. Although in principle, $j$ could take on many values, in practice I am interested in just two conditions (e.g., $j$ specifies whether the sample is for a pre- or post-intervention condition). In particular, the method I have developed attempts to establish correspondences between common features of the fitted density estimates (e.g., modes). I want to be able to plot curves for all the density estimates in such a way as to allow comparison between the two conditions, and between corresponding features of the distributions. I would like to create a grid of subplots (e.g., _m_ × _n_ if there are a total of 2 × _m_ × _n_ samples; note that _m_ and _n_ are only specified so that the overall plot has a convenient aspect ratio, they do not relate to the design of the experiment etc.). Each subplot would show two curves: one for the distribution underlying $\mathbf{x}_i^\mathrm{pre}$ and the other for $\mathbf{x}_i^\mathrm{post}$. Further, to allow the correspondences to be compared, I would like to color the line used to draw the curves using a “rainbow” (or other color map), so that one could, for example, look at the red regions of two curves and visually determine whether they do in fact correspond. To this end, I can provide vectors ( _x_ , p( _x_ ), _c_ ), where _x_ is an arbitrary value on the horizontal axis of the density estimate plot, p( _x_ ) is the density estimate for _x_ , and _c_ is a value in (0, 1), such that a particular value of _c_ should correspond across all density estimates, and should, therefore, be plotted in the same color (e.g., a value of 0.2 might sit at the peak of the mode of the distributions). Lastly, in order to allow the two curves in each subplot to be distinguished, I'd like to be able to style each differently, for example by using a thin line for the pre condition and a thick line for the post condition. Thanks.


title: How to combine two hazard ratios from the study for a meta-analysis?

text: I've been searching for a way to combine two hazard ratios from the same study for a meta-analysis, but found nothing. Does anyone know how to do this? Any thoughts would be great. Good Day, Simon


title: Updating the probability distribution when removing links from a Bayesian network

text: I have spent a fair amount of time trying to solve this problem but I can't find the solution. More specifically, I have the following matrix:               P(D|E&F) = [ 0.5 0.3 0.5 0.2 ;                   0.5 0.7 0.5 0.8 ]      All the variables are binary (two states) D, E and F are nodes of a Bayesian network. E and F are the parents of D and they are independent. Now, new evidence comes in and we know that D and E are independent as well (the link between D and F remains but the arc from E to D is removed). How do I go about finding P(D|F) from P(D|E&F) ?


title: How to interpret the model parameters of libsvm via MATLAB interface?

text: I used the MATLAB interface of libsvm for doing binary classification of 997-dimensional training data. I am trying to understand how the resulting model is used to compute the predicted output (which we get by calling `svmpredict`) The model contains fields (it has linear kernel):               nSV = [546; 246]; totalSV=792; rho = 0.093     and svCoeff [792x1 double] and SVs [792x997 double]      I thought that we must be simply multiplying svCoeff with SVs to get a [997x1] matrix which we then multiply with the actual feature, before shareholding by rho. But that's not the case. Can someone illustrate with a simple equation how these parameters are used to do classification?


How can I measure separability between different number of instance of one feature vector? For example the main vector is `V=[1 1 2 3 4 5 7 8 10 100 1000 99 999 54]` and three sub-vector with different sample lengths are `t1=[1 1 2 3 99 1000]` or `t2=[1 10 1000]` or `t3=[2 3 4 10 100 99 999 54]`. Which one is more separable and more informative? I mean I am looking for a sub-vector with minimum length and maximum separability power. If I put it in GMM, the sub-vector with less samples has better probability which is not fairly approach.               train=[1 2 1 2 1 2 100 101 102 99 100 101 1000 1001 999 1003];      No_of_Iterations=10;     No_of_Clusters=3;     [mm,vv,ww]=gaussmix(train,[],No_of_Iterations,No_of_Clusters);     test1=[1 1 1 2 2 2 100 100 100 101 1000 1000 1000];     test2=[1 1 2 2 100 99 1000 999];     test3=[1 100 1000];     [lp,rp,kh,kp]=gaussmixp(test1,mm,vv,ww);     sum(lp)     [lp,rp,kh,kp]=gaussmixp(test2,mm,vv,ww);     sum(lp)     [lp,rp,kh,kp]=gaussmixp(test3,mm,vv,ww);     sum(lp)      The results are as follow :               ans =          -8.0912e+05          ans =          -8.1782e+05          ans =          -5.0381e+05      I will really appreciate, if you could help me.


title: Quantitative methods and statistics conferences in psychology?

text: What are the most significant annual conferences focusing on quantitative methods in psychology? This could include but is not limited to psychometrics, mathematical psychology, and statistical methods in psychology. Rules:   * One conference per answer   * Include a link to the conference


title: Calculating statistical significance with unequal sample sizes and unequal variances

text: I have two samples, one with $n_1 = 41,000$ and the other with $n_2 = 881$; the larger sample has a standard deviation of $13.74$, and the smaller has an $SD=10.75$. The means are different, and when I run a Welch's t-test, I get a $p < .001$. I'm not sure if that's the appropriate test. I checked the skew for both samples; it was $29$ for the large sample and $9$ for the small one. Should I use a Mann-Whitney U test, or do I have enough data to assume normally distributed samples? In the end I need to know if the means of the samples are statistically different and be able to say one mean is $X$ times larger than the other.


title: What is the name of this normalization?

text: This may be a stupid question but... is there a specific name for normalizing some data so that it has mean=0 and sd=1? Or do I just say "data was normalized to have mean=0 and sd=1"? thanks nico


title: What is the expected value of the sample variance under a linear regression with omitted variables of an AR(2) process?

text: Lately, I have been interested in phenomenons related to omission of variables. For example, it can be shown that the expected value of the sample variance under the inclusion of one variable $x_1$ but omission of one variable $x_2$ is $\mathbb{E}(s^2|x_1,x_2)= \sigma^2 + \frac{\sigma^2}{n-1} RSS_{x_1, \beta_2 x_2}$ where $\beta_2$ is the true coefficient of $x_2$ and $RSS_{x_1, \beta_2 x_2}$ is the residual sum of squares when running a regression with $x_1$ as predictor and $\beta_2x_2$ as outcome. Now, I am interested in getting a better expression for this when we have an $AR(2)$ process. Let an $AR(2)$-process be given by $x_t = ax_{t-1}+bx_{t-2} + \epsilon_t$ (with the $\epsilon_i$ being independent and normal with standard deviation $\sigma$). If we run a regression with $(x_2,\ldots, x_T)$ as the response vector and $(x_1,\ldots, x_{T-1})$ as the predicting variable, what can be said about the residual sum of squares? That is, if we (falsely) think that the process is an $AR(1)$-process, what can be said about the expected residual sum of squares? Since I am quite unschooled in statistics (my only academic background is in mathematics), I would be interested in references (articles and books) as well as answers, even if they are only tangentially related to this question.


title: The order of variables in ANOVA matter, don't they?

text: Am I correct to understand that the order in which variables are specified in a multifactorial ANOVA makes a difference but that the order does not matter when doing a multiple linear regression? So assuming an outcome such as **measured blood loss, y** and two categorical variables   1. **adenoidectomy method, a**   2. **tonsillectomy method, b** The model _y~a+b_ is different to the model _y~b+a_ (or so my implementation in R seems to indicate). Am I correct to understand that the term here is that ANOVA is a **hierarchical** model since it first attributes as much variance as it can to the first factor before trying to attribute residual variance to the second factor? In the example above the hierarchy makes sense because I always do the adenoidectomy first before doing the tonsillectomy but what would happen if one had two variables with no inherent order?


title: Clustering Road Networks at various Radii

text: I have numerous road segments that belong to a city. Each road segment as a weight that is given to it due to its level of 'connectivity' (‘impact’from a small side street to a motorway). So, small neighbourhood road networks won’t have much weight at the regional level. These are calculated at different radii from each segment going from 400 meters at 400 meter intervals to 20,000 m. So each line segment will have a ‘measure’ at each radius (400, 800, 1200m) that indicate its level of impact on the total network. Now here’s my question: I want to be able to 'cluster' these so I can say that these segments 'cluster' around the R400 (400m) radius and therefore are local BUT also have activity around 800m and 1000m so could be local and a little peripheral activity. So below are 6 records just showing the activity at 400, 800 and 1200 meters - in reality it'll go up to 20,000m at 400m intevals.                 400m       800m        1200m          1 1.2982806 1.3231481 1.3453017 2 1.2048655 1.2961186 1.3369008 3 1.2557751 1.3028713 1.341819 4 1.3692737 1.367762 1.3970656 5 1.3519189 1.3815914 1.3904856 6 1.3940601 1.3942301 1.4031725 So the higher the score the greater the 'connectivity' at that level of radius. So 400m radius of the 2nd record start at 1.20 then at 1200m is 1.33 - so it has a higher connectivity at 1200m. Therefore, I think the 'shape' of the 'line' and how it clusters with similar lines is important. But I want the data to dictate this and not me starting to create 'categories' by forcing the data to do things that I want to see - this is my concern. I was thinking using the K Mean measure (But this asks for number of categories...which I think dictates to the data what I want to see) or comparing areas under the curves Or shapes of the curves once of plotted them out – basically grouping the once that look similar, but I know that there are issues where because they could be same areas and shapes but inverted. thanks for your patience, clearly, I am rather lost... would be grateful for any advice! with kindest regards, Atakan


title: Appropriate probability distribution for censored ratio data

text: I have two bacterial markers, I'll just call them _X_ and _Y_. _X_ codes for virulence, whereas _Y_ just indicates the bacteria is present. Consequently, _X_ will not show up without _Y_ although not all bacteria have _X_. Both also have a detection limit, although I only care about the detection limit for _X_ because I want to know the probability distribution for the ratio of _X_ / _Y_ (pathogenic to normal bacteria if present). Because the data are censored and a ratio, I was considering fitting a censored beta distribution. However I realized that the ratio itself is not censored - for example, _X_ could be below detection limit but still some large % of _Y_ if _Y_ is also fairly low. Does anyone have any suggestions for how to address this problem? I am assuming using a ratio distribution won't work because _X_ and _Y_ are obviously not independent.


title: prequential evaluation - classification

text: I perform prequential evaluation like this: start with a training set, classify a number of examples, then add the correctly classified examples in the training set and continue to classifying the next number of unseen examples. Is this supposed to increase performance as examples are added to the training set or this doesn't apply to every case? By increasing performance I mean if the average F1, not only accuracy, of the second piece of unseen examples must be higher than that of the first or is it possible for a latter piece to have worse average F1 than a former? And if it has worse, what does this could be possibly mean? Could it mean a problem with training data? This paper Sentiment Knowledge Discovery in Twitter Streaming Data describes prequential evaluation. It experiments with Naive Bayes and not SVMs, probably for the reasons mentioned in comments below. Thanks a lot in advance!


title: Test on the properties of two betas?

text: I have a set of observations. When I do an OLS I get a $\beta$. After I remove some values I also get a $\beta$. Now I want to test if the properties of two betas are the same? Please give me a hint how to do that?


title: Methods for imputation of missing values with spatially autocorrelated data

text: I am looking at the spatial patterns of turnover in aquatic assemblages using gradient forest and generalized dissimilarity models in R. I have species and environmental data for more than 400 sites sampled from streams across the country. However, I also have many missing values in my predictor (environmental) variables, up to 30% for some variables, so removing the cases with incomplete data sets or replacing them by the mean will result in bias and loss of information. About the missingness pattern, my data is really an assembly of data collected by different county administrations, and the decisions about which variables to sample are made at the county level. For example, some counties have monitored total and dissolved nutrients but others have only routinely monitored total nutrients. The missingness pattern is thus affected by those decisions. The next matrix is an example of how my data would look like:                    County V1  V2  V3  V4  V5  V6      [1,] 10     52  6   35  294 48  25      [2,] 10     22  7   41  53  42  NA      [3,] 10     118 NA  55  82  59  NA      [4,] 10     150 8   13  91  63  15      [5,] 10     500 NA  NA  NA  102 9       [6,] 9      58  7   NA  22  73  7       [7,] 9      9   6.5 NA  38  152 17      [8,] 9      9   7   NA  14  224 11      [9,] 9      142 5.5 NA  57  64  11     [10,] 9      90  6   NA  102 66  NA     [11,] 6      30  7   9   NA  NA  11     [12,] 6      420 4.5 8   NA  NA  NA     [13,] 6      43  4.5 3.5 NA  NA  NA     [14,] 6      50  6.5 116 NA  NA  14     [15,] 6      10  NA  13  NA  NA  8      >       where "County" is the different county administrations that have provided data, and "V1 to V6" are the environmental variables. In county 10 all variables are sampled but some are missing at random. In county 9 variable V3 is not routinely monitored. In county 6 variables V4, and V5 are not monitored, in addition there are missing values in the routinely measured variables due to, e.g. failure in the sampling device. I cannot forget to mention that my data are also spatially autocorrelated at small spatial scales (<10 to 100km). I would like to estimate those missing values but I don’t know which method is the most appropriate to do so, and which R packages are recommended.


title: Identifying populations and samples from a dataset

text: I'm writing a proposal for a restrospective cohort study using registry data. Am I right to think of the registry data as the whole population? If it is the population, then am I right in thinking there is no need to calculate required sample size (as is required for other cohort studies) or confidence intervals? If it is a sample, what would be the study population, and how would I go about calculating sample size? Additionally, with this in mind, how do I treat missing data? If there are 1000 heart transplant patients in a registry but only 500 of them have complete data (i.e. all the variables have values), do I treat these 500 patients as a sample? Furthermore, if it is a sample, it isn't random. Therefore is it worth calculating the confidence intervals from this actual sample size?


title: Expectation Subscript

text: What is meant by subscripting the expectation with some distribution, e.g. $\mathrm{E_{f}}[h(X)]$? If it's any help, here's the context: In M.C. Simulation, we wanted > $\mathrm{E[h(X)] = \int{h(x)f(x)dx}}$ so we used the law of large numbers and central limit theorem based on samples: > $\frac{1}{n} \sum_{i=1}^{n}h(X_i) \longrightarrow E_{f}[h(X)]$, with $X_{1}, > \ldots, X_{n} \sim f$


title: Can I detrend and difference to make a series stationary?

text: I have a dataset that is clearly increasing as time goes on (exchange rate of a currency, monthly data over 20 years), my question is: Can I detrend the data and then difference it also to make it stationary, if the detrending in itself doesn't achieve this? And if so, would this be considered twice differenced, or just detrended and once differenced?


title: In caret what is the real difference between cv and repeatedcv?

text: This is similar to question Caret re-sampling methods, although that really never answered this part of the question in an agreed upon way. caret's train function offers `cv` and `repeatedcv`. What is the difference in say doing:               MyTrainControl=trainControl(         method = "cv",         number=5,         repeats=5     )      vs               MyTrainControl=trainControl(        method = "repeatedcv",        number=5,        repeats=5     )      I understand `cv` breaks the set into k-folds (parameter `number`), and then starts over and runs it parameters `repeats` number of times. The only thing I could think of is that maybe regular `cv` with `repeats` uses the same exact indexes for the folds each time? essentially running the `cv` on the same exact folds each time, vs perhaps `repeatedcv` selects new folds each times? Can someone clarify?


title: How to analyze these data?

text: The biological data is listed as following:                  V1    V2    V3    V4    V5    V6     0.064 0.014 0.016 0.012 0.013 0.023     0.056 0.000 0.000 0.008 0.010 0.000     0.042 0.014 0.024 0.008 0.017 0.023     0.031 0.014 0.016 0.008 0.013 0.023     0.068 0.000 0.008 0.004 0.020 0.000     0.081 0.000 0.000 0.004 0.010 0.000     0.060 0.014 0.016 0.006 0.010 0.023      or you can download the data from http://www.mediafire.com/?6yp9l9m47jv433a.               A<- dat[,1]      B<- dat[,2:6]      I want to compare the difference between the first column to other columns of the data.Because only dat[,2] and dat[,6] not subject to normal distribute,I used wilcox.test instead of t.test function to caculate in R. But the warning messages rised up,such as "In wilcox.test.default(A, B[, 1]) : cannot compute exact p-value with ties". Could you give me some suggestions? Thank you.               wilcox.test(A,B[,1])          Wilcoxon rank sum test with continuity correction          data: A and B[, 1] W = 49, p-value = 0.00184 alternative hypothesis: true location shift is not equal to 0          Warning message: In wilcox.test.default(A, B[, 1]) : cannot compute exact p-value with ties


title: Calculating probability for bivariate normal distributions based on bootstrapped regression coefficients

text: Dear all, I was encouraged to ask this question here as well as on stackoverflow and would be very appreciative of any answers... Due to hetereoscedasticity I'm doing bootstrapped linear regression (appeals more to me than robust regression). I'd like to create a plot along the lines of what I've done in the script here. However the `fill=int` is not right since `int` should (I believe) be calculated using a bivariate normal distribution.   * Any idea how I could do that in this setting?    * Also is there a way for `bootcov` to return bias-corrected percentiles? sample script:               library(ggplot2)      library(Hmisc)      library(Design) # for ols()          o<-data.frame(value=rnorm(10,20,5),                   bc=rnorm(1000,60,50),                   age=rnorm(1000,50,20),                   ai=as.factor(round(runif(1000,0,4),0)),                   Gs=as.factor(round(runif(1000,0,6),0)))           reg.s<-function(x){               ols(value~as.numeric(bc)+as.numeric(age),data=x,x=T,y=T)->temp          bootcov(temp,B=1000,coef.reps=T)->t2               return(t2)          }           dlply(o,.(ai,Gs),function(x) reg.s(x))->b.list      llply(b.list,function(x) x[["boot.Coef"]])->b2           ks<-llply(names(b2),function(x){          s<-data.frame(b2[[x]])          s$ai<-x          return(s)          })                ks3<-do.call(rbind,ks)      ks3$ai2<-with(ks3,substring(ai,1,1))           ks3$gc2<-sapply(strsplit(as.character(ks3$ai), "\\."), "[[", 2)                k<-ks3      j<-dlply(k,.(ai2,gc2),function(x){          i1<-quantile(x$Intercept,probs=c(0.025,0.975))[1]          i2<-quantile(x$Intercept,probs=c(0.025,0.975))[2]               j1<-quantile(x$bc,probs=c(0.025,0.975))[1]          j2<-quantile(x$bc,probs=c(0.025,0.975))[2]               o<-x$Intercept>i1 & x$Intercept<i2               p<-x$bc>j1 & x$bc<j2               h<-o & p          return(h)          })           m<-melt(j)      ks3$int<-m[,1]             ggplot(ks3,aes(x=bc,y=Intercept,fill=int)) +       geom_point(,alpha=0.3,size=1,shape=21) +       facet_grid(gc2~ai2,scales = "free_y")+theme_bw()->plott      plott<-plott+opts(panel.grid.minor=theme_blank(),panel.grid.major=theme_blank())      plott<-plott+geom_vline(x=0,color="red")      plott+xlab("BC coefficient")+ylab("Intercept")


title: Interpreting interaction effects in a multilevel model

text: I conducted a multilevel analysis using repeated measures across four time points. The model contains intercept, linear slope and quadratic slope. I am interested in examining the extent to which variable A interacts with other identified variables in my model, so I added the interaction terms to intercept, linear slope and quadratic slope to estimate the fixed effects of these interactions. I then eliminated the non-significant interactions and retan the model until the interaction terms in the model were all significant. My questions are:   1. When delete non-significant interactions from the model, if the interaction was significant on linear slope but not quadratic slope, should I still keep the non-significant one on quadratic slope? Or should I only keep the significant one on linear slope?   2. How should I interpret an interaction on linear slope but not quadratic slope (or vice versa)? Thank you for your response in advance! Any advice is appreciated.


title: How to treat binary outcomes with categorical explanatory variables

text: I am trying to apply a GLM in R. I have a binary response (success vs failure), and 3 categorical explanatory variables : Sex (male or female), Food (present or absent) and Wind (none, low, high). I arranged my data to end up with a count of number of success for each possible combination of explanatory variables. ( _*_ Sorry for the pictures instead of codes, but I wasn't able to format it so it could be understandable...) ![enter image description here](http://i.stack.imgur.com/ZZd6m.png) I tried to run a GLM on this, taking the Success column as the response variable and treating it as count data : `model2<-glm(Success~Wind*Sex*Food, data=data, family=poisson)` But the summary of it gives something unexpected : The residual deviance shows 0 degrees of freedom. ![enter image description here](http://i.stack.imgur.com/8H0ko.png) Why is it the case ? (i.e. is it unusual as I think or I'm just freaking out for nothing) Am I using an inapropriate model or distribution ? Is the trouble coming from the way I managed my dataset ? I don't know if I am using an appropriate technique... Basically, I want to know if any of my explanatory variables or their interaction has an effect on the success of the trial. I don't know how else I could get that info... Thanks


title: What is the distribution of $\chi^n_k$?

text: $\chi^n_k=\sum_{i=1}^kx_i^n$ where $x_i$ are Gaussian variables and $n>2$?


title: Hidden Markov Models relationship

text: I have a question regarding a small investigation that I have been conducting into the relationship between the length of observation sequence, T, on which two decoders (BCJR and classic Viterbi) decode in a fixed small (2 state, 3 symbol) HMM. When I plot $T$ against the average decoding rate accuracy, I get unsurprising results, the BCJR algorithm moderately outperforms the Viterbi on symbol-by- symbol decoding and the data shows essentially no relationship between $T$ and accuracy rate. However, for the same sequences, I have plotted T against $\log P(O,Q\vert\lambda)$, for both the sequence from the BCJR decoder and the sequence from the Viterbi decoder, where $O$ is the observation sequence. I have also plotted $T$ against $\log P(O\vert Q,\lambda)$. In both cases, for the same set of results, there is a 'perfect' negative linear relationship between $T$ and these log probabilities. This would accord to a 'perfect' relationship: as the length of observation sequence increases, the log probabilities exponentially decrease and the exponential decrease does not vary over $T$, which would be surprising to say the least. Does anyone have any reasons as to how we know these results are wrong? What would be the expected relationship between $T$ and these log probabilities? Any help, no matter how small or trivial - even if just an idea or a thought, would be really appreciated as I just would like to understand this situation further.


title: Cognitive processing/interpretation of data visualizations techniques

text: Does anyone know of research which investigates the effectiveness (understandability?) of different visualization techniques? For example, how quickly do people understand one form of visualization over another? Does interactivity with the visualization help people recall the data? Anything along those lines. An example of visualizations might be: scatter plots, graphs, timelines, maps, interactive interfaces (like Parallel Coordinates) etc. I'm particularly interested in research within a lay-person population.


title: Univariate- Variance preserving, order reversing transformation

text: This is a soft question: How can the order of a sample univariate data be reversed while preserving the variance?


title: avova one way vs multiway

text: Hello I'm running anova with R and I am wondering what the differce is between multiway anova and single anova. I know single anova gives a different answer than multiway anova but I dont know how to quantify multiway anova. Could someone explain in words the difference between:   * base~col1   * base~col1+col2   * base~col1*col2 http://www.gardenersown.co.uk/Education/Lectures/R/anova.htm#anova_1 How do I interpret the results of col1+col2 and col1*col2 because col2 is giving a P(>F) of 0.473 on single anova but P(>F) 4.21e-07 in multiway anova. col2 does correlate well with base having only a 2% coloration for 711 individuals. According to http://vassarstats.net/rsig.html 2% coloration isn't significant. Thank you edit: i know my col1 and col2 have some interactions is that causing the low multiple anova p value? i am using               aov.ex2 <- aov(base~col1+col2,data=data1)      summary(aov.ex2)     aov.ex2 <- aov(base~col1,data=data1)      summary(aov.ex2)     aov.ex2 <- aov(base~col1,data=data1)      summary(aov.ex2)      and the results are not simular


title: How does one estimate the number of trees and the depth when using random forests?

text: I am using random forests for a class. I am predicting weight training. In scikit learn I always used a rule of thumb of a depth of 5, the depth x features rounded. I had 53 features, rounded to 50 x 5 and got 250 for number of trees. I picked 250 trees and 5 node depth and got great results. I got 98% accuracy for my test set (100% for training). The data set is for coursera machine learning. The data set is at d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv I removed seven features that were unrelated. The predictor is classe. Well, it works again, but I have no formal basis. I would appreciate any guidance,


title: Factorial design for interactions with ordinal DV

text: I want to make a 4 question, multiple choice survey in which each question asks about an analogous range of actions for slightly different scenarios. Each participant would only be administered the survey once. For example: "In situation W how would you behave? (5[very aggressively] -- 1 [not aggressively])." "In situation X how would you behave? (5[very aggressively] -- 1 [not aggressively])." "In situation Y how would you behave? (5[very aggressively] -- 1 [not aggressively])." "In situation Z how would you behave? (5[very aggressively] -- 1 [not aggressively])." I then wanted to sum the responses to the 4 questions, and use that total score as my DV. I plan to have 2 IVs. The first will be relationship status, and consist of 7 levels (e.g. Married, single, divorced, dating, widowed, etc). The second will be, say, partner hair-color preference, and will consist of 3 levels (e.g. blond, brunet, redhead). So, my understanding is that this will be a 3x7 factorial design, with a total number of 21 individual conditions, and an ordinal DV. Ideally, I'd like to do a 3x7 factorial ANOVA to test for main effects and interaction effects, but I don't think my summed, ordinal survey response DV is amenable to this. Is there some other analysis i could use in this situation? Would the Kruskal-Wallis be appropriate? Thanks!


title: Logistic regression model internal validation in SPSS

text: Please help! I have developed a logistic regression model for prediction of a binary outcome at a range of different times following a set point. This gave me different models with different variables lists, although many variables overlapped and were seen to be significant at different times. I've then created a combined model to use the same set of variables at all time points and have run this for each of the 6 times in question. This model is therefore not the best-fit for any individual time point, but is the best-fit across all time points. I've been looking at how to internally validate this model within SPSS, and am stuck. I have been able to bootstrap to obtain variable estimates at each of the time points, but am looking for a measure of overall model performance (i.e. a mean AUC to calculate model optimism or similar). However, my attempts at 10-fold cross-validation have failed as I cannot work out how to apply the model created on 9/10th of the dataset to the other 1/10th. Every time I save an .xml file, SPSS tells me it is not recognised and cannot be used. The additional problem is that the model was tested on different population numbers at different time points, due to missing data issues. Is there a simple way to obtain overall model performance data to internally validate a logistic regression model in SPSS, given the problems I have outlined? Thanks


title: Testing certain contrasts: Is this provably a hard problem, or not?

text: I posted this to mathoverflow and no one's answering: Scheffé's method for identifying statistically significant contrasts is widely known. A **contrast** among the means $\mu_i$, $i=1,\ldots,r$ of $r$ populations is a linear combination $\sum_{i=1}^r c_i \mu_i$ in which $\sum_{i=1}^r c_i=0$, and a scalar multiple of a contrast is essentially the same contrast, so one could say the set of contrasts is a projective space. Scheffé's method tests a null hypothesis that says _all_ contrasts among these $r$ populations is $0$, and given a significance level $\alpha$, rejects the null hypothesis with probability $\alpha$ given that the null hypothesis is true. And if the null hypothesis is rejected, Scheffé points out that his test tells us _which_ contrasts differ significantly from $0$ (I'm not sure the Wikipedia article I linked to points that out). I would like to know if one can do something similar in a different sort of situation. Consider a simple linear regression model $Y_i = \alpha + \beta x_i + \varepsilon_i$, where $\varepsilon_i\sim\operatorname{i.i.d.}N(0,\sigma^2)$, $i=1,\ldots,n$. The null hypothesis I want to consider concerns a different sort of contrast. It says there is no subset $A\subseteq\lbrace 1,\ldots,n\rbrace$ such that $E(Y_i) = \alpha_1 + \beta x_i$ for $i\in A$ and $E(Y_i) = \alpha_2 + \beta x_i$ for $i\not\in A$, where $\alpha_1\ne\alpha_2$. If the subset $A$ is specified in advance, then an ordinary two-sample $t$-test does it, but we want something that considers all subsets and holds down the probability of rejecting a true null hypothesis. One could figure this out if efficiency were not a concern: find a test that goes through all $2^{n-1}-1$ possibilities. Even then it's problematic; two contrasts would not be independent. I asked an expert on outlier detection about this and he just said it's a combinatorial nightmare. Then I asked if one could _prove_ that there's no efficient way to do it, perhaps by reducing an NP-hard problem to it. He just said he stays away from NP-hard problems. So: Can one prove either that this problem is "hard" or that it's not?


title: Is a snapshot of last login times a type of censored data?

text: I have a database of users with a field of the last time they logged in. Would this kind of data be considered "censored", and if so is there anything I can do with it beyond a histogram/density of how many people logged in within the last month, two months, etc? (To be clear: I only know the last time they logged in, so they could have logged in every day for years, up until a week ago, and all I know is they haven't logged in since last week.) It feels to me that this is censored data, but in a reverse sort of way from the definitions I read online. But perhaps this is just left-censored data? Or is this some kind of logit problem? For now, all I'm able to say is the percentage of users who have logged in during the last month, two months, etc, but I know people will want to hear more than that, and I also fear that it's easy to misinterpret what this data means. (It feels very particular, like explaining a confidence interval.) Ultimately, I may have to (manually) dump the data on a weekly basis to be able to take it any farther. Any references appreciated.


title: Trend estimation for participation?

text: I have logs for users and posts in a blog platform for 3 years. I can easily find out how many posts each user have made per day/month/year/etc. What I want to find out is if the frequency of posts is growing or not, in other words, if the users are contributing more over time or not. New users could join at any time and old ones could abandon the system (only the first one is logged). Is a trend estimation (http://en.wikipedia.org/wiki/Trend_estimation) suitable for this ? More over, is there a python package for this kind of time analyse ?


title: Purpose of the link function in generalized linear model

text: What is the purpose of the link function as a component of the generalized linear model? Why do we need it? Wikipedia states: > It can be convenient to match the domain of the link function to the range > of the distribution function's mean What's the advantage of doing this?


title: Finding a leading indicator of a time series

text: I am interested in finding a leading indicator of $Y_t$. Is it sufficient to find a variable $X$ for which its lagged value is correlated with $Y$? Do I have to give consideration to the spurious regression phenomenon? Or is it correct to say that when it comes to forecasting, I can ignore this problem?


title: Relationship between expected score and expected rank in a round robin tournament

text: Suppose you have $n$ players play a round robin tournament. A win adds 1 to the winning score of a player while a loss adds nothing. After the round robin tournament is finished the players are assigned ranks based on their scores, with ties being broken at random, the top player receiving a rank of $n-1$ and the bottom player receiving a rank of 0. What is the relationship between the expected value of the $\textbf{score}$ of player $i$ and the expected value of the $\textbf{rank}$ of player $i$? Sample result of a tournament with 5 players: $scores=(3,3,2,1,1)$ and $ranks=(3,4,2,1,0)$ Addition information: A player $i$ has an attribute $\lambda_i$. The probability that player i beats player $j$ in a game is $\frac{\lambda_i}{\lambda_i+\lambda_j}$. Numerical simulations suggest that taken across players, the two quantities perfectly correlated, but I don't know how to prove this. Any references would be appreciated.


title: Conjugate prior for a Gaussian model with shifted variance

text: Consider a set of observations $ \{ y_i \}$ and assume a Gaussian model for these data: $y_i \sim \mathcal{N}(\mu, \sigma^2)$. Suppose the mean parameter $\mu$ is known, but the variance parameter $\sigma^2$ is not. The conjugate prior for $\sigma^2$ is then the scaled inverse chi-squared distribution; consequently, we have a nice posterior. Now, suppose instead of $\mathcal{N}(\mu, \sigma^2)$ our model is $\mathcal{N}(\mu, \sigma^2 + \sigma_0^2)$ where $\sigma_0^2$ is known. What can we say regarding the prior and posterior distributions of $\sigma^2$ in this case? Thank you. Regards, Ivan


title: Can monthly seasonality adjusted data still exhibit seasonality?

text: I am working on monthly seasonally adjusted data. To the best of my knowledge (maybe I'm wrong), data is seasonally adjusted to take out any form of seasonality. After doing an ACF of the second difference of my data with R, I saw a spike from the lag one point that crosses the horizontal dashed line. This led me to suggest a nonseasonal MA(1) component. I then saw a lag from the 12th point that lies right on the horizontal dashed line. Can this mean there is a seasonal MA(1) component with period 12? And if so, can monthly seasonally adjusted data still exhibit seasonality?


title: How do I incorporate an innovative outlier at observation 48 in my ARIMA model?

text: I am working on a data set. After using some model identification techniques, I came out with an ARIMA(0,2,1) model. I used the `detectIO` function in the package `TSA` in R to detect an _innovative_ outlier (IO) at the 48th observation of my original data set. How do I incorporate this outlier into my model so I can use it for forecasting purposes? I don't want to use the ARIMAX model since I might not be able to make any predictions from that in R. Are there any other ways I could do this? Here are my values in order:               VALUE <- scan()       4.6  4.5  4.4  4.5  4.4  4.6  4.7  4.6  4.7  4.7  4.7  5.0  5.0  4.9  5.1  5.0  5.4       5.6  5.8  6.1  6.1  6.5  6.8  7.3  7.8  8.3  8.7  9.0  9.4  9.5  9.5  9.6  9.8 10.0       9.9  9.9  9.8  9.8  9.9  9.9  9.6  9.4  9.5  9.5  9.5  9.5  9.8  9.3  9.1  9.0  8.9       9.0  9.0  9.1  9.0  9.0  9.0  8.9  8.6  8.5  8.3  8.3  8.2  8.1  8.2  8.2  8.2  8.1       7.8  7.9  7.8  7.8      That is actually my data. They are unemployment rates over a period of 6 years. There are 72 observations then . Each value is to at most one decimal place


title: When to divide data into training & test set in logistic regression?

text: I am using Logistic Regression in a `low event rate` situation.   Overall universe: 46,000   Events: 420 Conventional logistic regression models divide the data into training and test sets and compute the error rates. The final coefficients and threshold levels are chosen and a model is created. OTH, I'm just trying to prove that so and so coefficient is significant and has positive association with the event in study. I'm not developing a model as of now. I don't focus on error rates (too many true negatives!) and chose my threshold level ~ hit rate. Should I consider dividing my universe into 2 samples, the conventional way? With such a low event rate, I'm worried that doing this to bias my coeff. estimates.


title: One tailed Fisher's exact test and the hypergeometric distribution

text: I'm working on a software and need to do some statistics. In my case I'm comparing two different networks that have a certain number of nodes with attributes and want to test with a Fisher's exact test if some of these attributes are enriched in the smaller one.   I also have a good math library that is able to calculate hypergeometric distributions and want to use it. So after reading a bit it seemed to me that you need to calculate the probabilities for your obtained data and more extreme cases and then add these. If this is correct and my table looks e.g. like               -------------------------------------     |                 |Network1|Network2|     -------------------------------------     |has Attribute    |a=7     |b=10    |     -------------------------------------     |has not Attribute|c=17794 |d=1107  |     -------------------------------------      How more extreme tables do I have to calculate for a one-tailed test? I guess until a=0, b=17, c=17801 and d=1100? Thanks a lot!


title: Why is variability measured relative to a point?

text: Why are measures of dispersion calculated relative to some central point? Why wouldn't, for instance, all possible non-repeated, pairwise differences in the dataset be a valid measure of spread?