title: Order of Matrix Operations in Mahalanobis Calculations

text: I'm teaching myself to translate equations to code after many years of letting my math skills atrophy, and am trying to do it on my own as much as possible. I've run into a couple of difficult equations though that I just can't solve without a few pointers though, particularly the one for Mahalanobis distances mentioned at http://en.wikipedia.org/wiki/Mahalanobis_distance (where the equation depicted below was taken from). ![enter image description here](http://i.stack.imgur.com/2vEXF.png) ![enter image description here](http://i.stack.imgur.com/YqRyr.png) ![enter image description here](http://i.stack.imgur.com/O6tha.png) The main thing that's slowing me down is the order of operations, in large part because the are matrix operations mixed in. $S^{-1}$ is an inverse covariance matrix, which I've already managed to write code for. The square root it trivial, of course. I'm not sure if I'm reading the part under the square root correctly, but my first guess is that I should compare every single data point in the X vector to its corresponding mean, transpose that resulting matrix and then multiply it by the inverse covariance matrix, then multiply that result by another matrix consisting of each X vector value compared to the corresponding mean again. If I understand this correctly, it will return one Mahalanobis distance for each X value in the vector. Is there anything in this order of operations I might be doing incorrectly? Since this is a more advanced topic that some of the stats I've coded to date, it's difficult to find examples on the Web with dummy data I can check my own code against, let alone examples of how to calculate it out step by step. If anyone knows off-hand of any such examples, that also might be helpful. I have the impression that anyone familiar with such equations would find it trivial to solve them, but I don't have any experience with 'em yet.


title: Interpreting p-value significance

text: If you have performed a statistical test using the standard 95% confidence threshold then a p-value of < 0.05 will indicate significance. I was reading a blog post about typesetting some correlation data using the `R` library `xtable`. The author uses a series of `ifelse` statements to replace p-values with stars to indicate whether the result was less than `0.001`, `0.01`, or `0.05`. I was wondering if anyone could comment on the validity of this? If you perform a statistical test at a 95% confidence level and you get a p-value of `< 0.01` does that mean that you could say that the result is significant to the 99% level. Or can you only assume a stronger confidence level if you explicitly performed the test with that criteria?


title: Contrast for hypothesis test in R (lmer)

text: I'm running a lmer mixed effects model with a four-level factor (levels "0","10","100","1000") as the fixed effect.               lmer(free ~ reward.f + (1|S), longdata)      I know that by default, R uses treatment contrasts and the levels 10, 100, and 1000 are compared to level "0". I would instead like each level to be compared to the previous one, to test a monotonic decrease in "free" across levels of "reward.f" I would do this:               contr.mat <- matrix(c(c(-1,1,0,0),c(0,-1,1,0),c(0,0,-1,1)),4)     colnames(contr.mat) <- c(10,100,1000)     contrasts(longdata$reward.f) <- contr.mat         contrasts(longdata$reward.f)          10 100 1000     0    -1   0    0     10    1  -1    0     100   0   1   -1     1000  0   0    1      Is this the correct contrast matrix for these comparisons?


title: Random multinomial sample, with some restrictions

text: Given a list of $k$ non-negative integers $m_i$, $i=1,...,k$, and a positive integer $n$, with $n\ge m +k$, where $m=\sum m_i$. I need to generate $k$ random positive integers $n_i$, such that $\sum n_i = n$, $n_i > m_i$, and $n_i=1$ whenever $m_i=0$. What's the simplest way of accomplishing this using an ordinary multinomial generator?


title: Is there a multivariate version of logistic regression?

text: Based on readings with logistic regression, it appears that you could use this analysis to make predictions about categorical variables. Does logistic regression allow you to predict multiple dependent variables with one independent variable, or just one? For example, if I were looking at predicting the gender of someone based on their emotional intelligence scores, could I also predict gender and race at the same time? How does this work? ### Reference Stevens, J. P. (2009) _Applied Multivariate Statistics for the Social Sciences_ (5th Edition) New York: Routledge Academic. Chapter3.


title: "NAs introduced by coercion" in seqecreate

text: I am using `seqecreate` to specify an event sequence object using `TraMineR`. The dataset has 1271 rows and 4 columns, an id, timestamp (timed), and 2 event variables (one numeric and one string). Most of the data have only 1 event per id. I've done other analyses with these data in SAS without issues. When I run `seqecreate` the procedure stops at id 425 with an error message:               > etv.ev <- seqecreate(id=etvev$id, timestamp=etvev$timed, event=etvev$event2, use.labels=FALSE)          Warning message:          In seqecreate.internal(data = data, id = id, timestamp = timestamp,  :       NAs introduced by coercion      The offending data in this part of the dataframe look ok (row 535 is the last before the error):               row     id       event2 timed event1     534 2144155702   xETV     0      1     535 2149998801   xETV     0      1     536 2150098301   xETV     0      1     537 2150098301   ETV2    27      2     538 2152463101   xETV     0      1      This error occurs with both event variables. I'd like to understand the problem and how to fix it.


title: Choice of referent twin in twin difference model

text: Carlin (2005) points out that mixed effects models specifically for twin data can be simplified by calculating differences between paired clusters. This allows for modeling specifically the within cluster effects and cancelling out any between cluster effects. http://ije.oxfordjournals.org/content/34/5/1089.short The mixed model: \begin{equation} E(Y_{ij} | X_{ij}) = \beta_0 + \beta_W (X_{ij} - \bar{X}_{i\cdot}) + \beta_B\bar{X}_{i\cdot} \end{equation} Becomes \begin{equation} E(D_i^{Y}) = \beta_W D_i^{X} \end{equation} with $D_{i}^X = X_{i1} - X_{i2}$ The results of this analysis depends on how you code twins as 1 and 2. The author simply states, "The oldest twin is chosen as the referent [i=1] twin" without mentioning whether it is a global rule or a deliberate choice for this analysis. I think there must be some kind of causal framework aspect to this. It seems like, for some analyses (especially post-infancy), the choice of order could be different and this could give vastly different results. How do you chose who is twin #1?


title: Comparing two distributions

text: I have more than five hundred thousands samples of a continuous variable measured in two groups: a treatment and a control one. I would like to decide whether the measurements follow the same (unknown) distribution in the two groups. Graphically this seems to be the case: the density functions (obtained using R) for the two groups overlap entirely. I tried to asses this in a more formal way using the Kolmogorov-Smirnov test. However I obtain a p.value of 2e-13 which suggests that the the two groups do not follow the same distribution. I am not convinced though and think that this is due to the fact that I have so many measurements that even the slightest difference leads to reject the null hypothesis. I tried to check if the mean is the same (using Wilcoxon-Manney test) and again the p.value is very low p < 1e-6. However the difference between the mean values in the two groups is quite low (0.006 and the values can range from 0 to 1) which is, for each practical purpose, identical. Am I using the wrong statistical tests? How can I assess in a formal way whether the two distributions are the same or not?


title: Why is my t-test intuition so far from reality?

text: Interpretation of t-test value (2 tailed): ~20% chance that these results could have been obtained randomly. My intuition: MEAN_2 - MEAN_1 = ~3% of MIN(STD_1, STD2). Calculate area under normal normal PDF between -3% and +3% of std from mean which is ~2.5%. Why is the t-test result so far off from my intuition??? * * * COUNT_1 = 3572 MEAN_1 = 702 STD_1 = 2926 * * * COUNT_2 = 4596 MEAN_2 = 621 STD_2 = 2730 * * * t = 1.26 df = 7403


title: Incorporating non-collinearity into RANSAC affine transform estimation

text: The MATLAB code at this link can be used to estimate the affine transform between two sets of 2D points. The code makes the implicit assumption that neither of the point sets contain collinear points. However, I have data in which either the source or destination point sets (p0 and p1 in the code) _can_ contain collinear points. Instead of trying to somehow remove collinearity in the sets, I wish to modify the RANSAC procedure so that collinear points become outliers. In the initial model guess phase(code snippet from link below), 3 point pairs are selected at random. If I were to check here for collinearity of point pairs, would it **guarantee** that collinear points will be treated as outliers at the end of the procedure ?               for k = 1:num_guesses       % Select 3 random constraints       c = [ceil(rand * n)];       for j = 2:3          k = ceil(rand * n);          while any(c == k)             k = ceil(rand * n);          end          c = [c; k];       end         ....(more code)      end % loop on k          Some additional details: Analysis of my data reveals that even if the input points(p0) are not collinear, using the above code ~~will~~ can select inliers in such a way that output(p1) points can be collinear. For example, in the toy data set at this link, no two points of the same color should be selected. Put another way, for each color, 2 points must become outliers. Is it possible at all to enforce non-collinearity into this RANSAC based procedure?


title: Which is more appropriate in this case -- PCA or Factor Analysis?

text: Suppose I have an experiment that has 8 factors. These 8 factors are probably related and can hopefully be reduced. For each combination of these 8 factors that I test, I get a single output. My goal is to reduce the 8 factors to a smaller set of variables. In this case is PCA or Factor Analysis more appropriate? I am not too clear on the differences between the two methods. My limited understanding is that PCA is best to reduce multiple outputs into a smaller set of the most important outputs and that Factor Analysis is best for finding relations between many factors. If this is true I would expect Factor Analysis to be the best method for this situation.


title: Should coin flips be modeled as Bernoulli or binomial draws in RJags?

text: What is the best way to model coin flips as a hierarchical model? Do you say coin draws are a series of draws from Bernoulli trials or as one draw from a binomial distribution? That is something like this:               model {       p ~ dunif( 0, 1 )       for( i in 1 : n) {         h[i] ~ dbern( p )       }     }      or, this:               model {       p ~ dunif( 0, 1 )       numberOfHead ~ dbinom (totalTrials, p)     }


title: How to determine appropriate number of features and also which features to select?

text: So I have a dataset which I am using K fold cross validation on to select the number of features and which features should be selected. As I understand it, I would set the number of features to be selected as a constant, then let my feature selection method pick whatever N features it wants to for each fold of cross validation. E.g. Set N = 1, evaluate CV performance for all classifiers I tested with a single feature Set N = 2, Repeat with 2 features... Select the N with the best CV performance. So I now have K sets of N features each. To train a final model on all of the data, I need to settle on the "best" N features. Do I select (1) the pair of features which appeared together the most out of K folds, or do I select (2) the top N features which were selected the most out of all K folds? You can think of (2) as a histogram of all the features which were ever selected in all K folds and each fold provides a vote if it picked that feature. (1) would be like (2) except each bin is a pair of features, rather than a single feature.


title: Effect of features that are highly correlated with each other on a decision tree

text: I have a dataset of roughly 500 features and am training a binary classifier using GBM - gradient boosted machines, an ensemble of decision trees. Of these 500 variables, I am sure some are highly correlated with each other, though probably not to the extent where they are linearly dependent. For example, one variable might be average age of people in city X which was collected by survey 1, and another variable is the average age of people in city X collected by survey 2. How does such a massive set of features affect the decision trees? In the regression setting, this should increase prediction variance, but can also be mitigated by regularization.


title: Testing for population similarity on non-parametric data

text: I'm looking at similarity of samples for non-ordinal/non-parametric data. I've included a dummy dataset below for reference in the question. Neither A, B, or C is ordinal data: it's entirely categorical. Where C is filled in, it shows a Pareto distribution. The actual dataset is large and a very small subset of the population has C's data filled in (< 10%). Essentially, the columns are not completed in the dataset. However, I'm trying to determine if there the "shape" of A (and, independently, B) is different for where C is filled in and where C is not filled in. Phrased another way: "Do id's who have a marked C behave the same way in relation to (either/and A, B) as id's who don't have a marked C". And, "do id's who have a marked C behave the same way (as indicated by A,B) as the entire set of id's?" My readings[1] indicate that what I am looking for is Pearson's Chi-squared test, as I will be testing a non-parameteric dataset against another dataset (in this case, the same dataset in itself); there is no average/stddev meaningful for A,B, or C. [1] I'm boning up on analytics stats on the fly; this isn't my area of expertise.               id  A   B   C     1   0   1   a     2   0   0   a     3   0   1   a     4   0   0   a     5   0   1   a     6   1   0   a     7   1   1   a     8   1   0   a     9   1   1   a     10  1   0   b     11  0   1   b     12  0   0   b     13  0   1   b     14  0   0   c     15  1   1   c     16  1   0   c     17  1   1   d     18  1   0   e     19  0   1   f     20  0   0   g     21  0  1  null     22  0  1  null       ... null C's continue on for a long time.


title: OLS with a heavily skewed independent variable

text: I am regressing a log-normally distributed dependent variable (`wage`) on a heavily skewed independent variable and I want to make sure I handle it in the best possible way. The independent variable `award` is the percentage point (measured 0-100) likelihood that a given person will win a specific award in their career. The measure is heavily skewed to zero (~70.00% of the data takes a 0 value), however, it does not look log normal and not quite bimodal. There is no other value with more than 1.00% of the population density, but many values are around 0.25%, and a good number are between 0.50% and 1.00%. The current model is an OLS regression, `log_wages = b0 + b1*award + b2*control_1 + b3*control_2 + b4*control_3`. Is this the best model? Do I need to transform `award` in any way? Are there any tests I should run to test the distribution?


title: How to denoise a "Poissonous" time series

text: I have $N$ time series each of which can be modeled as $$y_{kt}=Ax_{kt}+b+\varepsilon_{kt}\quad(1\le k\le N,1\le t\le T),$$ where $x_{kt}\sim\text{Pois}(\lambda\Delta t)$ and $\varepsilon_{kt}\sim N(0,\sigma^2)$. Parameters $A$, $b$ and $\sigma^2$ are unknown and I want to extract the uncontamined version of all $y_{kt}$. Several techniques come into my mind, like ICA (two sources, uncontamined data and noise) or wavelet (perform DWT on data and keep the largest coefficients, for example, or use wavelet shrinkage, etc.). But I am not sure about that:   1. Can ICA handle poisson noise inherent in poisson processes (which is described by $\lambda$ but not $\varepsilon_{ki}$)?   2. Will wavelet denoising ruin the properties of the model? Or are there better choices? Any hints are appreciated.


title: Output of dredge and model.avg functions of MuMIn

text: I am running a GLM of about 15 predictor variables and using `dredge` and `model.avg` functions of `MuMIn`. The output, however, is confusing (attached below). I would appreciate your help understand it and/or organize it. Cheers!               > dd<-dredge(model1, beta = FALSE, eval = TRUE, rank = "AICc",     +                  fixed = NULL, m.max = 15)     There were 50 or more warnings (use warnings() to see the first 50)          > dd<-dredge(model1, beta = FALSE, eval = TRUE, rank = "AICc",     +                  fixed = NULL, m.max = 15)     There were 50 or more warnings (use warnings() to see the first 50)          > model.avg(get.models(dd, subset = delta < 4))          >Call:     model.avg.default(object = get.models(dd, subset = delta < 4))          >Component models:      ‘1/2/3/4/5/8/9/10/11/12’     ‘1/2/3/4/5/7/8/9/10/11/12’   ‘1/2/3/4/5/9/10/11/12’      ‘1/2/3/4/5/6/8/9/10/11/12’   ‘1/2/3/4/5/6/8/9/11/12’           ‘1/2/3/4/5/6/8/11/12’        ‘1/2/3/4/5/8/9/10/11’        ‘1/2/3/4/5/6/7/8/9/11/12’    ‘1/2/3/4/5/6/7/8/11/12’      ‘1/2/3/4/5/9/10/11’               ‘1/2/3/4/5/7/9/10/11/12’     ‘1/2/3/4/5/8/10/11/12’       ‘1/2/3/4/5/6/7/8/9/10/11/12’ ‘1/2/3/4/5/7/8/9/10/11’      ‘1/2/3/4/5/6/9/10/11/12’          ‘1/2/3/4/5/10/11/12’         ‘1/2/3/4/5/8/10/11’          ‘1/2/3/4/5/6/11/12’          ‘1/2/3/4/5/10/11’            ‘1/2/3/4/5/6/8/11’                    >Coefficients:              (Intercept)         lensweightn               sexf2              sitef2              sitef3              sitef4             weightn              yearf2              2.486550299         0.136111214        -0.462421682         0.375999064         0.291747922         1.821189287        -0.097745693        -1.699335700                   yearf3  lensweightn:yearf2  lensweightn:yearf3        sexf2:sitef2        sexf2:sitef3        sexf2:sitef4       sexf2:weightn       sitef2:yearf2             -0.363010762        -0.080879247        -0.146453405        -1.305819037        -0.182946730        -0.348561920         0.041144419         0.430176195            sitef3:yearf2       sitef4:yearf2       sitef2:yearf3       sitef3:yearf3       sitef4:yearf3      weightn:yearf2      weightn:yearf3 lensweightn:weightn              0.203700885        -1.434464341        -1.251719523         0.416009083        -0.793133787         0.094364720         0.100236695         0.001200883        lensweightn:sexf2              0.028792820


title: What are the distributions of these random variables?

text: What would be the distribution of $y$, when:   1. $y = x^2$ and $x\sim\mathcal{N}(\mu, \sigma^2)$.   2. $y = x^2$ and $x\sim$ Log-$\mathcal{N}$.


title: Similarity matrix and multiple-regression

text: Let, $S_{n*n}$ represent a similarity matrix, among $n$ observation, my case n = 215. and $Y=\{y_1, y_2, ...,y_n\}$ contains a response value for each $x_n$ observation. For each observation we have $t$ features,my case t = 50. Since, the number of observation is not enough vs features to do a nice multiply regression, I am wondering, if there is any clever way to use that Similarity matrix ! I had an idea, which I would like to add, I was thinking that, for each test data $x^*$, I can look at my training similarity matrix,$S$ and take, say top 10 high similar observation $X_{top10}$, and try to guess $y^*$ for my test data ${x^*}$, from those 10 observation only. Theoretically, it makes a lot of scenes, due to similar observations are ending up at similar response values, but the problem is how to make that regression and also identify most predictive features.


title: Time series trend

text: I have a time series which has a very strong upward trend for the first half, then very strong downward for the second half and finishes pretty much back where it started. Should I split the data in two for analysis - or can I still account for a trend which is net neutral?


title: Logistic regression for time series

text: I would like to use a binary logistic regression model in the context of streaming data (multidimensional time series) in order to predict the value of the dependent variable of the data (i.e. row) that just arrived, given the past observations. As far as I know, logistic regression is traditionally used for postmortem analysis, where each dependent variable has already been set (either by inspection, or by the nature of the study). What happens in the case of time series though, where we want to make prediction (on the fly) about the dependent variable in terms of historical data (for example in a time window of the last $t$ seconds) and, of course, the previous estimates of the dependent variable? And if you see the above system over time, how it should be constructed in order for the regression to work? Do we have to train it first by labeling, let's say, the first 50 rows of our data (i.e. setting the dependent variable to 0 or 1) and then use the current estimate of vector ${\beta}$ to estimate the new probability of the dependent variable being 0 or 1 for the data that just arrived (i.e. the new row that was just added to the system)? To make my problem more clear, I am trying to build a system that parses a dataset row by row and tries to make prediction of a binary outcome (dependent variable) , given the knowledge (observation or estimation) of all the previous dependent or explanatory variables that have arrived in a fixed time window. My system is in Rerl and uses R for the inference.


title: Looking up and inserting values from a table

text: I have a huge dataset of reaction times from a psycholinguistic experiment, and in separate tables, I have information that I'd like to add in about each word's frequency and things like that. The main dataset looks something like this:               maindataset <-structure(list(item = c(101, 103, 102, 104, 104, 102, 103,     102, 101, 103, 104, 101), react = c(512, 510, 506, 499, 515, 516,     517, 518, 509, 599, 520, 523), prime = c(1, 2, 3, 1, 3, 1, 3,      2, 2, 1, 2, 3)), .Names = c("item", "react", "prime"), row.names = c(NA,      -12L), class = "data.frame")      And the frequency information looks like this:               frequencyinfo <-structure(list(item = c(101, 102, 103, 104),      frequency = c(10,30, 40, 50)), .Names = c("item", "frequency"),      row.names = c(NA,-4L), class = "data.frame")      So basically, for each value of item in maindataset, I want to look up the frequency from the frequency table, and store it in a new column in maindataset. I've been able to do it in the past using for loops, but it's horribly inefficient (it takes minutes to go through the whole 10,000 row dataset). Anyone know the best way to go about this?


title: Fixed effects have df=0 and chi-sq=0 when tested using anova()

text: I have a made a linear mixed model using `lmer()`. I began with a maximal model containing 6 terms and all two-way interactions (15 total), I also have a random effect of subject ID. Of these terms, 11 interactions were significant and so were kept in the final model. I found this by dropping interactions and comparing to previous model using `anova()`. As all 6 fixed terms are part of a significant interaction they remained in the model. My problem is that when I tried to test the fixed terms to find a chi-squared value for them by removing the term and comparing to the full model, the `anova` output said `chisq = 0, df=0` but that `p< 2.2e-16`. I don't understand how p can be so small if `chisq = 0`. I also don't understand why `df=0` which would suggest I am comparing the same model when actually I have removed a term. What does this mean? Have I fit too many interactions to my data (I only have 160 data points). Is there any point in finding the chi- squared value for these fixed effects?


title: Chart for seeing trends, correlations, triggers and patterns

text: I've scraped a set of data, and I was wandering what the best way to start analyzing it was. My initial thoughts were to use a multi y axis graph, like this: ![http://imgur.com/bFbDIxz](http://i.stack.imgur.com/wrexa.png) The reason I was thinking about a multi axis graph was because if they were using the same axis some of the data would be way overstretched for the axis. That is, some of the data sets are 10-100 others are 10,000-100,000, but I wasn't sure if this is best suited. I'm comparing tweet sentiment of a particular stock, to the performance of that stock, I've got up to 7 data sets to compare -   * tweet volume   * sentiment-pos   * sentiment-nue   * sentiment-neg   * stock trade volume   * stock buy price   * stock sell price Any idea of the best way to visualize this, is to allow me to view patterns and correlations? Also is this usually done on a chart or can it be done by running `VBA` or other scripting language to find correlations in the data?


title: Approximation of partition function (normalizer)

text: Say we have a nasty probability distribution like, $$ P(x) = \frac{P^*(x)}{Z} $$ where we can easily compute $P^*(x)$ for a given $x$ but not $P(x)$ because partition function $Z$ is expensive to compute. There are approximation techniques like Monte Carlo which can be used to estimate the expectation of functions under this distribution. Instead of some expectations, if we want a probability of some $x$, can we get it using some sampling techniques? Can we approximate the partition function $Z$ itself? I can't see resources which do these approximations, probably because it is not a standard thing to do.


title: Popular implementations of one vs all logistic regression for mutiple classes in R

text: Is there any popular implementation of logistic regression in R that defaults to using a one vs all approach for a categorical dependent variable with more than two classes? I have specific reasons to favor a one vs all approach vs simply performing a multinomial logistic regression. Also, while I know I can write my own routine to iteratively perform the classification myself, for the purposes of my project it is better if I don't do that. I'd be grateful for the knowledge!


title: Simple definition of elliptical distribution

text: I'm looking for a simple explanation of what an elliptical distribution is and how it differs from a non-elliptical distribution. My knowledge of stats is very basic; Wikipedia was not much help on this topic (too technical for my level of understanding). Thanks.


title: Method for generating correlated non-normal data

text: I'm interested in finding out a method for generating correlated, non-normal data. So ideally some sort of distribution that takes in a covariance (or correlation) matrix as a parameter and generates data that approximates it. But here's the catch: the method I'm trying to find should have the flexibility to also control its multivariate skewness and / or kurtosis. I'm familiar Fleishman's method and the use of the power method of normal variates, but I believe most of those extensions only allow the user for certain combinations of _marginal_ skewness and kurtosis, leaving the multivariate skewness / kurtosis just out there. What I was wondering is if there is a method that helps specify the multivariate skewness and / or kurtosis, alongside with some correlation / covariance structure. About a year ago I took a seminar on copula distributions and I remember the professor casually mentioning that through the use of vine copulas, one could generate data which is, say, symmetric in each one of its 1-D marginals but jointly skewed and vice-versa. Or, even further, that any lower-dimensional margins could have some skewness or kurtosis while keeping the highest dimensions symmetric (or not). I was marveled by the idea that such flexibility could exist I've been trying to find some sort of article or conference paper that describes said method but I have been unsuccessful :(. It doesn't have to be through the use of copulas, I'm open to anything that works. Edit: I have added some R code to try to show what I mean. So far I am only well-acquainted with Mardia's definition of multivariate skewness and kurtosis. When I first approached my problem I naively thought that if I used a symmetric copula (Gaussian in this case) with skewed marginals (beta, in this example), univariate tests on the marginals would yield significance but Mardia's test for multivarite skewness/kurtosis would be non-significant. I tried that and it didn't come out as I had expected:               library(copula)     library(psych)     set.seed(101)          cop1 <- {mvdc(normalCopula(c(0.5), dim=2, dispstr="un"),                  c("beta", "beta"),list(list(shape1=0.5, shape2=5),                  list(shape1=0.5, shape2=5)))}                      Q1 <- rmvdc(cop1, 1000)                 x1 <- Q1[,1]                 y1 <- Q1[,2]               cop2 <- {mvdc(normalCopula(c(0.5), dim=2, dispstr="un"),                  c("norm", "norm"),list(list(mean=0, sd=1),                  list(mean = 0, sd=1)))}                      Q2 <- rmvdc(cop2, 1000)                 x2 <- Q2[,1]                 y2 <- Q2[,2]          mardia(Q1)            Call: mardia(x = Q1)          Mardia tests of multivariate skew and kurtosis     Use describe(x) the to get univariate tests     n.obs = 1000   num.vars =  2      b1p =  10.33   skew =  1720.98  with probability =  0     small sample skew =  1729.6  with probability =  0     b2p =  22.59   kurtosis =  57.68  with probability =  0          mardia(Q2)     Call: mardia(x = Q2)          Mardia tests of multivariate skew and kurtosis     Use describe(x) the to get univariate tests     n.obs = 1000   num.vars =  2      b1p =  0.01   skew =  0.92  with probability =  0.92      small sample skew =  0.92  with probability =  0.92     b2p =  7.8   kurtosis =  -0.79  with probability =  0.43      Upon inspecting the contours for 'cop1' VS 'cop2' as well as the empirical bivariate density plots, I can also see that none of them look symmetric at all. That's when I realized this is probably a little more complicated than I thought. I know that Mardia's is not the only definition of multivariate skewness/kurtosis, so I'm not limiting myself to finding a method that only satisfies Mardia's definitions. thank you!


title: Generating Bayesian network graph with dsc file

text: I'm using a free version of Bayesian network software called Netica. It allows only 15 nodes for the free version. Do you know any other software or R package that generate a kind of graph below using a dsc file? ![network graph](http://i.stack.imgur.com/uetJn.png) I know R package called bnlearn has a function read.dsc for dsc files, but am not sure how to make a graph. If you know how you can get R to make this kind of or similar graph (with probabilities of influence in each node box), please walk me through R scripts.


title: Canonical correlation analysis with continuous and binary data

text: I came across interesting article on application of canonical correlation analysis (CCA). Authors apply classical CCA on a mixed variables dataset (both independent and dependent sets include continuous and binary variables). To sum up their main contribution: "usual" correlation matrix (i.e., based on Pearson's correlation) which is particularly designed for continuous measurements is not appropriate, so they proposed "new" measure to capture the correlations between different data-type variables. This part of the analysis seems clear to me. I'm struggling to grasp the procedure to calculate canonical variate scores and canonical loading in the context of mixed data. In usual settings (i.e., when all variables are continuous) canonical variate scores are found by multiplying raw data with the canonical weights. Canonical loadings are subsequently found by correlating the raw variable scores with the variate scores. I wonder to know, how to derive variate scores in these particular settings. It's proper to simply multiply canonical weights with the corresponding raw matrix to obtain variate scores? Thanks in advance for any suggestions.


title: Univariate time series analysis: sine wave with a constant period?

text: > **Possible Duplicate:**   >  Sine wave with a constant period: what else is constant? Which statistic that you can calculate on a window remains constant on a sine wave with a constant period? This is for a window function for a proprietary filter.


title: Fitting a pdf against Weibull pdf

text: I have a pdf function as follows: $$\dfrac{1}{s+a-b} [bs e^{-bt} + (a-b)(s+a)e^{-(s+a)t}]$$ I want to fit this against a weibull pdf with shape=1.12 and scale=461386. I want to calculate the values of s,a and b. Is there any standard tool for doing this? Prasenjit


title: Ordinal logistic model in R

text: I have the neuropsychiatric questionnaire scores of 300 individuals, of which 200 are normal, and 100 have the disease. The questionnaire is divided into 12 categories (delusion, agitation, ... etc). And each category begins with a question along the lines of: does the patient experience delusions? (yes=1, no=0, a binary variable). If yes, the person proceeds to specify the severity (1: Mild; 2: Moderate; 3: Marked) and frequency (1: Occasionally; 2: Often; 3: Frequently; 4: Very frequently), and these are ordinal Likert scale variables. To start off I have some trouble with organizing the data.frame. I want the binary variable, status, to be the response variable (0 = normal, 1 = disease) of my ordinal logit model. Below is how I envision part of my data.frame to look like (not sure if this is the best way to organize the data)                   status  delusion   severity   frequency  agitation  severity  frequency    ...         0       0          NA         NA         0          NA        NA           ...         0       1          2          1          0          NA        NA           ...         1       0          NA         NA         1          1         4            ...         0       0          NA         NA         0          NA        NA           ...         1       0          NA         NA         0          NA        NA           ...         .       .          .          .          .          .         .             .         .       .          .          .          .          .         .             .         .       .          .          .          .          .         .             .      What should I do with the NAs? Should I just code them as 0? Also, most (90%) of the participants do not experience any of the 12 symptoms (so most of my data.frame will consist of NAs). So would using a logit model be ideal here? I'm interested in seeing how each of the 12 symptoms (at what frequency and severity) contribute to the disease status of the participant.


title: How to use residuals to train a model?

text: I've read in a couple places about people using residuals from predictions of other models to improve the model's accuracy. Anyone have a link about how to carry this out in R, ideally with some example code? UPDATE: This is in the context of reading the papers about the Netflix prize where k-nearest neighbors (KNN) models were trained using residuals from other KNN models.


title: Sequential Sum of Squares, degree of freedom when number of variables greater than number of samples

text: I would like to use the Sequential Sum of Squares test. But the degree of freedom for the denominator is (n - p - 1), where n = number of samples, and p = number of variables in the full model. What do I do when p > n?


title: Measuring significance of differences in coefficients of regression models - same independents used to predict different independents

text: I am currently constructing a regression model meant to predict a given dependent variable - call it IN. I will also be constructing a model with the same independent variables to predict a second dependent variable - call it OUT. I would like to examine the difference between coefficients for the same independent variable across models. E.g., if the coefficient for VarX is 0.2 in the IN model, and 0.3 in the OUT model, I want to be able to test whether this difference is significant. What is the best way to do this? Additionally, I will want to run the IN model using slightly varied independent variables. My data varies over time, so I have fields for population size every year from 2003 to 2012. I want to run different versions of my model, using IN as my dependent variable, with each version using the population size of a different year. This will give me a series of coefficients - I want to be able to test whether the difference between coefficients for the "population size" variables is statistically significant. As above, what is the best way to do this? If answers to these questions already exist on the board, please direct me to them - I'm new to CrossValidated!


title: Bootstrapping residuals: Am I doing it right?

text: **First of all:** From what I understood, bootstrapping residuals works as follows:   1. Fit model to data    2. Calculate the residuals   3. Resample the residuals and add them to 1.   4. Fit model to new dataset from 3.   5. Repeat `n` times, but always add the resampled residuals to the fit from 1. Is that correct so far? * * * **What I want to do** is something slightly different: I want to estimate parameter and prediction uncertainty for an algorithm that estimates some environmental variable. What I have is a error-free time-series (from a simulation) of that variable, `x_true`, to which I add some noise, `x_noise`, in order to generate a synthetic dataset `x`. I then try to find optimal parameters by fitting my algorithm with the sum of squares `sum((x_estimate - x_true)^2)` (! not `x_estimate - x` !) as an objective function. In order to see how my algorithm performs and to create samples of my parameters' distributions, I want to resample `x_noise`, add it to `x_true`, fit my model again, rinse and repeat. Is that a valid approach to assess parameter uncertainty? Can I interpret the fits to the bootstrapped datasets as prediction uncertainty, or do I have to follow the procedure I posted above? Sorry if I'm not very clear with terminology, English isn't my first language and I'm pretty new to this. /edit: I think I haven't really made clear what my model does. Think of it as essentially something like a de-noising method. It's not a predictive model, it's an algorithm that tries to extract the underlying signal of a noisy time- series of environmental data. /edit^2: **For the MATLAB-Users** out there, I wrote down some quick & dirty linear regression example of what I mean. This is what I believe "ordinary" bootstrapping of residuals is (please correct me if I'm wrong): http://pastebin.com/C0CJp3d1 This is what I want to do: http://pastebin.com/mbapsz4c


title: Reproducing ARIMA error terms

text: When forecasting a moving average (MA) model using R's `forecast`, why does using `residuals(fit)` produce different results than `forecast` does? The question is very similar to (MA terms in arima), but the thread does not answer my question. It seems that, even if MA is invertible, using `residuals(fit)` won't produce the same result as `forecast`, except when the coefficient of MA is small. I understand that `forecast` is a wrapper of `predict`, which uses `KalmanForecast`. Does `KalmanForecast` use residuals at all to forecast under MA? How can I extract the error terms of a moving average (MA) model that is used by `KalmanForecast`? I would like to use an Excel sheet to forecast the future values using the R parameters. Can I just replace error terms with residuals?               library(forecast)         T <- c(-1.6749530, -1.8200980, -0.9128676, -1.0099365,  1.2365373, -1.9131794,                -1.4843611, -0.3176719, -0.4673834, -0.3048056) # Sample Data     fit <- arima(order=c(0,0,1), T) # Fit MA(1)     fit$coef[2] + fit$coef[1]*residuals(fit)[10] # Mannual Calculation     forecast(fit, 1) # forecast is a wrapper of predict     predict(fit, 1) # predict uses KalmanFilter     KalmanForecast(1, fit$model)$pred + fit$coef[2] #So, all of them produce the same results


title: How to do centering if I have a quadratic term?

text: I have been trying to run a multilevel model with both a linear and a quadratic term for income as my main variables of interest. It looks something like: \begin{eqnarray} &&y_{ij}=\beta_{0j}+\beta_{1j}\text{Income}+\beta_{2j}\text{Income}^2+\epsilon_{ij}\nonumber\\ &&\beta_{0j}=\gamma_{01}\text{Level-2 variables}+u_{0j}\nonumber\\ &&\beta_{1j}=\gamma_{11}\text{Level-2 variables}+u_{1j}\nonumber\\ &&\beta_{2j}=\gamma_{21}\text{Level-2 variables}+u_{2j}\nonumber \end{eqnarray} It has suggested to me by someone that it is useful to center the income variables for interpretation, yet I am not entirely sure how to do that when a quadratic term if present. I tried to first demean the linear term and then square it separately for the positive and negative part, but this would fix the tipping point at 0, which is not something I want. Is there a better way of doing this? Or should I give up on centering at all? Thanks JY


title: Why doesn't ML point estimate equal MAP point estimate even though I'm using uniform prior?

text: I asked a previous question about why the ML and MAP estimates are the same when using a uniform prior (How does a uniform prior lead to the same estimates from maximum likelihood and mode of posterior?) However, I am playing around with it some more and I tried a simple example which doesn't make sense to me. Let's say I flip a coin and it comes up heads. I now want to estimate the probability `p` of the coin coming up heads. Using maximum likelihood: we get p = 1/1 = 1 Using MAP estimate (with Beta(1, 1) prior, which is uniform): p = (1 + 1) / (1 + 1 + 1) = 2/3 So howcome the estimates aren't the same even though I'm using a uniform prior?


title: Orthogonal contrasts for four levels - are Helmert contrasts orthogonal?

text: I am planning to do a 2-way mixed ANOVA with the within-subjects factor having four levels. Since the exact data do not matter, I hope to suffice by only providing a general question. On my endeavours through the internet and different textbooks on statistics I have come to the point that I would like to compute this with ezANOVA from the package ez (as described by Andy Field in his "Discovering Statistics Using R") and since I will be interested especially in interactions between the 2-level between subjects factor and the 4-level within subjects factor I have decided to use type 3 sums of squares rather than the ezANOVA-standard type 2 sums of squares. In order to calculate type 3 SSs, however, according to Field, I need to construct orthogonal contrasts, since in case I do not, R will not really return type 3 SSs to be (but something else instead, I suppose, on which the textbook does not elaborate. Any explanations of what R does in case will be highly appreciated.) I have a couple of questions:   1. Can I generally compute a mixed ANOVA using ezANOVA and type 3 SSs without specifying any contrasts at all and still get interpretable results? I am left confused since the intriguingly easy usage of ezANOVA appears to allow this. However, as stated above, Field writes that one needs to use orthogonal contrasts for type 3 SSs (Andy Field, Discovering Statistics Using R, 2013 reprint, p. 476). IGNORE: 2. In case I truly need to construct orthogonal contrasts, how can I achieve this with a 4-level factor? I have tried to do this and I have come up with a Helmert contrast (as I have noted post hoc, when searching for a solution online). However, I am confused again, because as far as I know, the sum of products of the contrasts should add up to zero, which for Helmert contrasts does not hold, since the sum of products is not equal to 0:                   condition   C1      C2      C3      Cproduct             1        3       0       0       0             2       -1       2       0      -2             3       -1      -1       1       1             4       -1      -1      -1      -1                                       sum = -2      The only explanation that I could come up with for myself is that possibly the "sum of products = 0" criterion is only one way of being sure that contrasts are orthogonal, but that they might also be orthogonal even if this criterion does not hold. Or is there something special about Helmert contrasts that I am unaware of, some circumstance that qualifies them for orthogonality? STOP IGNORING Thanks in advance for any help. EDIT: Greg and Aniko have made me see the obvious mistake in basic arithmetics that I have made regarding the second question in my original post, so it might as well be ignored. Thanks for that one! However, maybe there is someone who could provide an answer to my first question regarding the behaviour of ezANOVA when using type 3 SSs without orthogonal contrasts. This would be great.


title: Tonga: low rank approximation of the natural gradient, question regarding Le Roux et al. 2007

text: This question concerns an implementation of the topmoumoute natural gradient (tonga) algorithm as described in page 5 in the paper Le Roux et al 2007 http://research.microsoft.com/pubs/64644/tonga.pdf. I understand that the basic idea is to augment stochastic gradient ascent with the covariance of the stochastic gradient estimates. Basically, the natural gradient approach multiplies a stochastic gradient with the inverse of the covariance of the gradient estimates in order to weight each component of the gradient by the variance of this component. We prefer moving into directions that show less variance during the stochastic gradient estimates: $ng \propto C^{-1} g$ Since updating and inverting the covariance in an online optimisation setting is costly, the authors describe a ``cheap'' approximate update algorithm as described on page 5 as: $C_t = \gamma \hat{C}_{t-1} + g_tg_t^T$ where $C_{t-1}$ is the low rank approximation at time step t-1. Writing $\hat{C}_{t} = X_tX_t^T$ with $X_t =\sqrt{\gamma} X_{t-1}\ \ g_t]$ they use an iterative update rule for the gram matrix $X_t^T X_T = G_t = $ ($\gamma G_{t-1}$      $\sqrt{\gamma} X^T_{t-1}g_t$ $\sqrt{\gamma} g^T_t X_{t-1}$      $g_t^Tg_t$) They then state ``To keep a low-rank estimate of $\hat{C}_{t} = X_tX_t^T$, we can compute its eigendecomposition and keep only the first k eigenvectors. This can be made low cost using its relation to that of the Gram matrix $G_t= X_t^T X_T$: $G_t = VDV^T$ $C_t = (X_tVD^{-\frac12})D(X_tVD^{-\frac12})^T$'' Because it's cheaper than updating and decomposing G at every step, they then suggest that you should update X for several steps using $C_{t+b} = X_{t+b}X_{t+b}^T$ with $X_{t+b} = \left[\gamma U_t, \ \gamma^{\frac{b-1}{2}g_{t+1}},...\ \ \gamma^{-\frac12}g_{t+b-1}, \ \gamma^{\frac{t+b}{2}g_{t+b}}\right]$ I can see why you can get $C_t$ from $G_t$ using the eigendecomposition. But I'm unsure about their update rule for X. The authors don't explain where U is coming from. I assume (by notation) this is the first k eigenvectors of $C_t$, correct? But if so, why would the formula for $X_t$ be a good approximation for $X_t$? When I implement this update rule, $X_t$ does not seem to be a good approximation of the ''real'' $X_t$ (that you would get from $X_t = [\sqrt{\gamma} X_{t-1}\ \ g_t]$) at all. So why should I then be able to get a good approximation of $C^{-1}$ from this (I don't)? The authors are also not quite clear about how they keep $G_t$ from growing (The size oft he matrix $G_t$ increases at each iterative update). I assume they replace $G_t$ by $\hat V\hat D\hat V^T$ with $\hat V$ and $\hat D$ being the first k components of the eigendecomposition? So in summary:   * I tried implementing this update rule, but I'm not getting good results and am unsure my implementation is correct   * why should the update rule for $X_t$ be reasonable? Is $U_t$ really the first k eigenvectors of $C_t$? (Clearly I cannot let $X_t$ and $G_t$ grow for each observed gradient $g_t$)   * This is far fetched, but has anyone implemented this low rank approximate update of the covariance before and has some code to share so I can compare it to my implementation? As a simple example if I simulate having 15 gradients in matlab like:               X = rand(15, 5);     c = cov(X);     e = eig(c);     %if the eigenvectos of C would give a good approximation of the original , this should give an appoximation of c, no?:     c_r = e'*e; %no resemblence to c      So I'm quite certainly doing it wrong, I guess U might actually not be the eigenvectors of C, but then what is U? Any suggestions or references would be most welcome! (sorry about the terrible layout, looks like only a subset of latex is supported, no arrays for matriced and embedded latex doesn't look particularly good, all the formulas are much more readable on page 5 of the referenced paper :) Also, is this considered off-topic here? It's really more related to optimisation and machine-learning...)


title: A Bayesian perspective on omitted-variable bias (and other covariate-selection bias problems)

text: As I know OVB, from a frequentist education, when you leave a variable $(z)$ out of your control set $(X)$ that is correlated with both your independent variable of interest (treatment $T$) and your dependent variable of interest ($Y$), your coefficient estimate will be biased because the explanatory power of the missing variable $z$ is distributed to the coefficients of included variables $(\hat\beta_i X_i$). **How does a Bayesian perspective view OVB?** For instance, if we use a data- based, rather than theory-based, variable-selection algorithm -- LR, AIC, BIC -- is it hard to conceive of OVB? Furthermore, how would one formally integrate the awareness of $z$ into our conditional probability statement? I mean, in Bayesian inference we want to estimate $P( \text{model} \mid \text{data}) = P(\theta \mid X)$. If we acknowledge some important but unobserved $z$, would we write $P(\theta \mid X, z)$? Furthermore, **how would a Bayesian perspective interpret other classes of covariate-selection bias problems?** I thinking about covariate-selection issues as elaborated by Pearl and others including,   * Pearl: Controling for an instrument can augment bias in linear models, and introduce and augment bias in non-linear models: http://arxiv.org/pdf/1203.3503.pdf   * Greenland: Conditioning on pre-treatment variables may induce bias (M-bias, aka Collider-Stratification bias): http://journals.lww.com/epidem/Abstract/2003/05000/Quantifying_Biases_in_Causal_Models__Classical.9.aspx (Sorry, paywalled)


title: Simulating correlated lifetimes at pre-specified correlation level

text: I am trying to simulate remission times for 100 patients from an exponential distribution with mean 1 year. I also want to simulate after-remission times for these 100 patients. But it is very much likely in the real life scenario that those who remit quickly are likely to have comparatively longer after- remission times and those who remit slowly are likely to have comparatively shorter after-remission times. So the after-remission time is likely to be correlated negatively to the remission time. I want to test the performance of an estimator at certain levels of correlation (say, `r=-0.5`, `r=-0.8` etc.). Could you please suggest how should I simulate the remission times and the after-remission times imposing the condition of certain levels of correlations? Should I consider the relationship to be linear or non-linear? I know from the real life data that the average after-remission time is about 1.5 years. How can I reflect this information in my simulated data, imposing the level of correlation?


title: Density function for a multivariate Bernoulli-like distribution

text: I'm looking for a distribution to model a vector of $k$ binary random variables, $X_1, \ldots, X_k$. Suppose I have observed that $\sum_i X_i = n$. In this case I do not want to treat them as independent Bernoulli random variables. Instead, I would like something like the multinomial: $P(X_1=x_1, \ldots, X_k=x_k) = f(x_1, \ldots, x_k; n, p_1, \ldots, p_k) = \frac{n!}{x_1! \cdots x_k!} \prod_{i=1}^k p_i^{x_i}$ but instead of the $x_i$ being nonnegative integers, I want them restricted to be either 0 or 1. I have been trying to see if the multivariate hypergeometric is appropriate, but I'm not sure. Thanks in advance for any advice.


title: Power analysis for moderator effect and mediation analysis

text: I want to do an a priori analysis of minimum sample size with G*Power and I have two questions about this: First, I will test an interaction effect for a moderator variable. I am not sure which statistical test to choose in G*Power as I am less interested in the overall $R^2$ and more interested in whether the interaction is significant. Second, I will do a mediation analysis following the Baron & Kenny approach. Which statistical test do I have to choose here to calculate the minimum sample size? I know about the G*Power tutorials, but I couldn't find anything about moderation- and mediation-analysis there. So I would be grateful if anybody could help me!


title: Any good book for learning probability programming

text: Are there any good books for me to learn probability programming? For example, I am new to Latent Dirichlet allocation (LDA) and Gibbs sampling. I have read some books about the techniques, but it still feels hard to develop programs to apply them. Are there any books that can help? I have taken some courses in machine learning, statistics and data mining. I am familiar with discriminant classifiers such as Support Vector Machines and the basic knowledge in generative classification. So I understand the concept in Naive Bayes Classification and Language models in Information Retrieval. I am new to the field of Latent variable analysis and graph modeling. Specifically, are there any books about how to develop programs related sampling, bayes inference? For instance, if I want to implement LDA using Gibbs sampling, any book can help? It is preferable if it has some easy understandable codes using toy examples.


title: How do betting sites update odds during a sporting match in real-time?

text: I recently played on a betting site (soccer games) and wondered how they change odds so quickly on sporting events. I get the basics of betting, but wonder how they change live betting odds so quickly? Do they use variables as morale, ball possession, ... ? And if thats the case how do they calculate it in real-time?


title: How to handle count data (categorical data), when it has been converted to a rate?

text: I am working on disease infection data, and I am puzzled on whether to handle the data as "categorical" or "continuous".   * "Infection Count"      * the number of infection cases found in a specific period of time, the count is generated from categorical data (i.e. no. of patient tagged as "infected")   * "Patient Bed Days"      * sum of total number of day stay in the ward by all patients in that ward, again, the count is generated from categorical data (i.e. no. of patient tagged as "staying in that particular ward")   * "infection per patient bed days"     * "infection count" / "patient bed days" both were originally count data, but now becomes a rate ### Question:   * Can I use Chi-Square here to assess whether the difference in "infections per patient bed days" is statistically significant or not? * * * ### Updates I have found that I can compare the incidence rate (or call it infection rate), but doing something like "incidence rate difference" (IRD) or "incidence rate ratio" (IRR). (I found it from here)   * What is the difference between IRD and t-test?    * Is there any statistical test complementary for IRR?


title: Likelihood of a set of ranks

text: I'm sorry if I cannot formulate the problem precisely. Let us consider an ordered set $Q$ with $n_q$ elements, and a number of its (possibly overlapping) subsets $B_1, \dots, B_k$. For a given $i$, one is interested in the question whether the ranks of set $B_i$ are uniformly distributed in $Q$, or whether the ranks of the elements of $B_i$ are lower than of the elements of $Q$. Note that the number of elements $n_i$ of the set $B_i$ can vary wildly. So, say below we have Q, sorted from rank 1 to 10. Elements of Q which do not belong to some $B_i$ are shown as "o", elements which _do_ belong to that $B_i$ are shown as "O": $B_1$: ooOooOooOO   $B_2$: OOOOOooooo   $B_3$: OOoOoooooo   In the real-world application, $Q$ has $10^4$ to $10^5$ elements, and the sizes of $B_i$ vary between 10 and $10^4$. Some of the $B_i$ are clearly biased compared to a more or less uniform distribution -- like in example $B_2$ above. Some of them less so, like $B_3$ above. My ultimate goal, however, is to order $B_1 \dots B_k$ by their deviation from a uniform distribution of ranks in $Q$. One possibility is just to use a Mann-Whitney U test on each of $B_i$, and sort $B_i$ by the respective p-values. The problem is that the power of this test will be very different for different $n_i$. Thus a subset with many elements can slightly deviate from the uniform distribution and nonetheless have a p-value lower than a much smaller set with a larger deviation. Consider this example: $Q = \{ 1 \dots 1000 \}$ The elements of $B_1$ are the first 10 elements of $Q$: $B_1 = \{ 1 \dots 10 \}$, $n_1=10$, $p_1 < 10^{-7}$ The elements of $B_2$ are sparsely distributed within the first 500 elements of $Q$: $B_2 = \{ 50, 55, 60, \dots, 450\}$, $n_2=81$, $p_2 < 10^{-15}$ So if I order them by $p$-value, I will often get results which have a modest effect size (i.e., modest difference in ranks), but a large number of elements in that particular $B_i$, **ranked before** other, more interesting cases where size of $B$ is smaller. Another option is to sort the $B_i$ according to an estimate of the effect size, for example, the U statistics. This, however, gives random flukes at the beginning of the list, which are not interesting. What is interesting? It is a combination of two things: statistical significance and effect size. As a biologist, I can only describe my practice: I am _not_ interested in something which has a large effect size, but is non- significant. Likewise, I'm not interested in something with a small effect size, but a small p-value. Look, imagine you have been doing a number of t-tests and you ask yourself the question: which of them are giving me the most relevant (i.e., interesting for the experimentalist) results? The ones with a low p-value may nonetheless correspond to a very small difference in means, while the ones where the means show a huge difference might be not significant at all. If I want to sort these results, what sorting should I adopt? I thought that maybe instead of applying a test-like or test-based procedure, I could calculate AIC values or evidence ratios for the different subsets. For this, I would have to obtain a likelihood function for the subsets, and I don't know how to tackle it. **Which is my actual question.**


title: How to estimate the change in prevalence of a trait in the population from non-random samples?

text: Suppose I am interested in the prevalence of drug use in a population over time, but I am unable to conduct audits (sample some transactions at random and check to see if drugs were exchanged). I do have data on seizures, needles or vials found in the park, as well as crime hotline reports, but these are not perfect because the policing technology, drug delivery methods, and reporting norms change, so these cannot reliably tell me about changes over time. If all these were immutable, I could implement a dynamic factor model where all these variables are a function of one or more unobserved latent factors. How is this problem usually handled? Are there any statistical procedures that can identify how the latent factor changes over time?


title: LMS cost function vs cross entropy cost function in neural networks

text: What is difference between using various cost functions: LMS,Cross entropy in neural networks? All of them have same derivative w.r.t final activation and hence all the gradients are still gonna remain same?


title: Should I standardize my variables for this particular case of cluster analysis?

text: I'm trying to cluster a list of records based on a (percentage) frequency distribution of variables which add up to 100%. Like   1. Record1 - VarA(25%) VarB(25%) varC(50%) varD(0%)   2. Record2- VarA(50%) VarB(15%) varc(0%) VarD(35%)  and so on. I have standardized variables before while dealing with different dimensions(lengths and weights) etc. In this case I do not think standardization is appropriate. Also is k means clustering appropriate in this context, I wanted to use k means and use the distribution observed at the centroid of the cluster for the whole cluster. Thanks a ton.


title: Time series analysis with neural networks

text: I'm new to neural networks and machine learning and I was wondering how you use time series data to set the weights of a regular FNN, and how you use the ending weights to forecast the time series. In essence how can you transform the time series data into weights and back again for the output.


title: How to compare or normalize two sets of probabilities?

text: I'm a stats newbie so I hope this question isn't too easy for you! :) I'm not sure how to ask it so I've written two versions: the short version is a pure stats question and the longer version explains my application. ## _The short version:_ I have two sets of probabilities. The first set has a range from 5.8x10-3 to 2.52x10-9 and the second ranges from 3.51x10-4 to 1.59x10-6. Both sets have a quasi-logarithmic distribution. The sets are not the same length. The sets are comprised of similar data but have different ranges due to the sample sizes. I need to normalize the probabilities somehow so that I can compare my input data against both. How do I do this? (remember I am a stats newbie!) If you need more info, please read the longer version: ## _The longer version:_ I am writing a software application that tries to predict the word a user is typing before he/she is finished. I have two source "dictionaries":   1. The first is a list of individual words ("singlets") + a count (the number of times the word appears in a corpus). It contains around 65,000 unique English words.   2. The second is a list of word triplets (e.g. "one of the") and the count of each triplet in a similar corpus. This list contains the most common 100,000 English triplets. Although it has more entries than the first list, it has far fewer unique words. The probability of a dictionary entry is: count / Σcount Both lists have quasi-logarithmic probability distributions, as you can imagine. I must predict words as the user is typing. For example, if the user just typed "one of t", my goal is to predict what the third word might be. The reason I use both lists is because the triplet list provides more accuracy but the singlet list has more unique words. He might be typing a common triplet like "one of the" or he might be typing an uncommon triplet like "one of tomorrow's". As a result I need to combine results from both lists. My result set is the most common singlets that start with "t" and the most common triplets that start with "one of t", sorted by probability. The problem is the probability ranges are very different (see above) due to the different sizes of the corpi and the nature of singlets vs. triplets, so my results are usually skewed toward one list or the other. I don't fully understand the mathematics behind this, but the bottom line is that the predictions are screwy. Thanks in advance...


title: How to plot the contribution of each regression coefficient in a model, with R?

text: I've fit a model with lm(), and I now want to analyze it to see where it's overfit, etc... I'm imagining a plot that has the index of each observation on the x-axis, and the corresponding responses plotted as points on the y-axis. In (vertical) line with each point is a set of colored lines stacked on top of each other. Each vertical, colored line corresponds to a predictor in the model, and it's height is the coefficient on that predictor applied the data point it's in line with. Is there already a function that does this? I don't want to reinvent the wheel (particularly because parsing the output of lm looks unpleasant).


title: error in getting predictions from a lme object

text: I am trying to get predictions for observations from a lme object. This is supposed to be quite straightforward. Yet, since I am get different types of errors for different trials, it seems to me I am missing something. My model is the following:               model <- lme(log(child_mortality) ~ as.factor(cluster)*time +              my.new.time.one.transition.low.and.middle + ttd +              maternal_educ+ log(IHME_id_gdppc) + hiv_prev-1,              merged0,na.action=na.omit,method="ML",weights=varPower(form=~time),              random= ~ time| country.x,              correlation=corAR1(form = ~ time),              control=lmeControl(msMaxIter = 200, msVerbose = TRUE))      It runs fine, fits the data well and the results make sense. Now to get predictions I have tried the following:               test.pred <- data.frame(time=c(10,10,10,10),country.x=c("Poland","Brazil",                 "Argentina","France"),                      my.new.time.one.transition.low.and.middle=c(1,1,1,0),                  ttd=c(0,0,0,0),maternal_educ=c(10,10,10,10),                  IHME_id_gdppc=c(log(5000),log(8000),log(8000),log(15000)),                     hiv_prev=c(.005,.005,.005,.005),                   cluster=c("One Transition, Middle Income","One Transition,                     Middle Income","One Transition, Middle Income","Democracy,                   High Income"))     >     > predict(model,test.pred,level=0)               Error in X %*% fixef(object) : non-conformable arguments      If I exclude, say, France, and only include countries in which cluster="OneTransition, Middle Income" then I get a different error               # create a toy data set     test.pred0 <-         expand.grid(time=20:29,country.x=c("Poland","Brazil","Argentina"))     z0 <-as.data.frame(cbind(my.new.time.one.transition.low.and.middle =                               c(0,0,0,0,0,0,1,2,3,4), ttd=c(0,0,0,0,0,0,1,0,0,0),                              maternal_educ=seq(from=10.0, to=12.0, length.out=10),                              IHME_id_gdppc=log(seq(from=5000, to=8000, length.out=10)),                              hiv_prev=rep(.005,10),                              cluster=rep("One Transition, Middle Income",10)))          z <- rbind(z0,z0,z0)     test.pred <- cbind(test.pred0,z)     # check     head(test.pred)     >  time country.x my.new.time.one.transition.low.and.middle ttd     > maternal_educ    IHME_id_gdppc hiv_prev     > 1   20    Poland                                         0   0     >   10 8.51719319141624    0.005     > 2   21    Poland                                         0   0     > 10.2222222222222 8.58173171255381    0.005     > 3   22    Poland                                         0   0     > 10.4444444444444 8.64235633437024    0.005     > 4   23    Poland                                         0   0     > 10.6666666666667 8.69951474821019    0.005     > 5   24    Poland                                         0   0     > 10.8888888888889 8.75358196948047    0.005     > 6   25    Poland                                         0   0     > 11.1111111111111 8.80487526386802    0.005     >                         cluster     > 1 One Transition, Middle Income     > 2 One Transition, Middle Income     > 3 One Transition, Middle Income     > 4 One Transition, Middle Income     > 5 One Transition, Middle Income     > 6 One Transition, Middle Income          # run the predictions     predict(model,test.pred,level=0)     > Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :     >   contrasts can be applied only to factors with 2 or more levels      In this example, the problem is due to cluster="One Transition, Middle Income" all the time. I don't understand why this is a problem. If I want to get predict() to work I have to include all variables from the model, right? Obviously, the input data in the model's call will not include factor set to the same values for all cases. Yet, if I want to get predictions just for subset of the data, or for new observations, I may be interested in only in cases where some factor is always set to be the same. Does it make sense? How can I get predictions in that case?


title: Proper way to combine conditional probability distributions of the same random variable conditioned on a discrete variable ? (based on assumptions)

text: My question is the following. Let's say I have two probability distributions: $f(x|b), g(x|c)$ $b$ and $c$ are discrete events while $x$ is a continuous variable, i.e., when the button b is pressed there is some distribution for the amount of rain fall the next day, $x$. When the button $c$ is pressed there is a different distribution of rain fall the next day, $x$. Are there any strategies for estimating the distribution of rain fall if both buttons are pressed, i.e., $h(x|b,c)$ ? And, what assumptions do those strategies rest on?


title: How should mixed effects models be compared and or validated?

text: How are (linear) mixed effects models normally compared against each other? I know likelihood ratio tests can be used, but this doesn't work if one model is not a 'subset' of the other correct? Is the estimation of the models df always straightforward? Number of fixed effects + number of variance components estimated? Do we ignore the random effects estimates? What about validation? My first thought is cross validation, but random folds might not work given the structure of the data. Is a methodology of 'leave one subject/cluster out' appropriate? What about leave one observation out? Mallows Cp can be interpreted as an estimate of the models prediction error. Model selection via AIC attempts to minimize the prediction error (So Cp and AIC should pick the same model if the errors are Gaussian I believe). Does this mean AIC or Cp can be used to pick an 'optimal' linear mixed effects model from a collection of some un-nested models in terms of prediction error? (provided they are fit on the same data) Is BIC still more likely to pick the 'true' model amongst candidates? I am also under the impression that when comparing mixed effects models via AIC or BIC we only count the fixed effects as 'parameters' in the calculation, not the actual models df. Is there any good literature on these topics? Is it worth investigating cAIC or mAIC? Do they have specific application outside of AIC? Cheers


title: Most interesting statistical paradoxes

text: Because I find them fascinating, I'd like to hear what folks in this community find as the most interesting statistical paradox and why.


title: Visualizing two scalar variables over time

text: What would be the best way to display changes in two scalar variables (x,y) over time (z), in one visualization? One idea that I had was to plot x and y both on the vertical axis, with z as the horizontal. Note: I'll be using R and likely ggplot2


title: Parametric techniques for n-related samples

text: The title is quite self-explanatory - I'd like to know if there's any other parametric technique apart from repeated-measures ANOVA, that can be utilized in order to compare several (more than 2) repeated measures?


title: Post hoc selection of important features in random forest?

text: I want to guarantee a parsimonious random forest (few features used). What are methods to do this? It was suggested to me to get the feature importance after the model was created, and then create a new model using only the important features. Is this a good idea/what other options are there?


title: Checking If two experiments match

text: I have a noob question for you. I have run two experiments: A and B. Between the two experiments I changed an element which should _not_ change the results. I want to check that the results of experiment B do not contradict the results of experiment A. The results for each experiment are a series of integer counts. I would expect that the results would be approximately normally distributed. I was thinking of using the following method:   1. Calculate mean and standard deviation of experiment A's results   2. The null hypothesis is that the "independent" variable _did not_ affect experiment B   3. Use p-value and the normal distribution confidence interval [erf(n/sqrt(2))] to calculate the probability that a sample as extreme as the sample in B could appear in experiment A   4. Multiply the calculated probability of all samples together to get the p-value Is there a better approach?


title: Getting R square value from ctree

text: Is there a way to retrieve the R square value from a `ctree` model? Or is this impossible with how the function runs? Thanks!


title: How to use ppoints to generate points within 0 and 0.05 for qq plotting in R?

text: I ran Tassel3 and I filtered results with p-value not more than 0.05. This way, it is ok to draw a Manhattan plot. However, for a QQ-plot there is a problem. Say I have 40 thousands SNPs, after filtering, about 1,800 SNPs were kept in the output file. Now how do I draw QQ-plot because of partial p-values? `ppoints` only generate points within (0,1). I think the resulting points should be in the filtered range (0,0.05) with the same number (1,800) too. `runif` can do the trick, but it is random, each time it generates different results. I think this would affect the QQ-plot line. I figured out three ways to solve it but I am not sure if it is right. Hope you guys can give some advices. the following are three methods:   1. `0.05*ppoint(1800)`   2. `ppoints(1800*100/0.05, a=0.05)`, then I only get the former 1/20 results.   3. `ppoints(40k)`, and use those p-values below 0.05 in plotting.


title: Interpreting effect size

text: One of the most widely used interpretation of effect size is Cohen's $d$, as follows:   * $<.10$: _trivial_   * $[.10,.30]$: _small to medium_   * $[.30,.50]$: _medium to large_   * $>.50$: _large to very large_ This appears to be the standard in psychology. **I am wondering if there is such a standard in sociology or social science in general?** I know the size of effect is best interpreted in the context of the research!


title: Goodness of fit of linear model for Segmentation of GPS positions time series

text: I have some **GPS coordinates series** taken in regular time steps and I need to verify whether some chunks of the trajectories fit well as a **straight line** or not. The aim is to perform **segmentation** on each trajectory, I am already able to find hard curves quite well, but the other segments are straight lines and smooth curves, I'd like to separate them based on a criterion of goodness of fit to a straight line. The whole trajectories are not so big (usually not much bigger than 50 miles or kilometres) so spherical geometry corrections and haversine distances are not really required for this, but they would be interesting if simple. I thought about fitting an $Y=aX+b$ model (longitude $X$, latitude $Y$) and looking at the $r^2$, but if a trajectory is nearly south to north I could have different results than the same trajectory rotated by 90 degrees since the values of $a$ and $b$ would be quite unstable. Linearly fitting each coordinate versus time, hoping that if both fits well then I'd have a good fit may solve this issue, but I don't know if this is a valid approach nor if $r^2$ would be the best measure neither how to combine those $r_x^2$ and $r_y^2$. Does anyone know a more appropriate way of verifying the goodness of fit of a straight line in such GPS series?


title: Why do we typically visually assess our assumptions?

text: For example, in linear regression, we assume the errors are normally distributed, have constant variance and have an expected value of 0. These assumptions all refer to the population. But we look at our actual data points to see if these assumptions are satisfied, Why do we do this? How do we know that results hold for the population? If we only have 10 data points, why don't we say that we don't have enough data to test whether these assumptions are satisfied?


title: How to design neural networks for pattern recognition in biometry?

text: Having read numerous texts regarding neural networks and their characteristics, I am getting increasingly confused, paradoxically – I am looking for a brief explanation or references to the right sources. I am trying to implement neural networks using PyBrain to **recognise patterns in biometric data and classify them**. It is just a small project to actually learn the procedures and the goal is to recognise which subject provided which data. I am using a feedforward multilayer perceptron trained using back-propagation and the results are successful with roughly 0.7 probability but honestly I am just guessing the settings like types of layers and such. I realise this question is very broad and very difficult, but could you perhaps **suggest the best approach to such task** and possibly **identify the variables that determine which algorithms to use**? The data are in the time domain so for example, the typical input for training and following classification would be a vector of 256 values representing 2 seconds of ECG input. ~~I should also point out that the same data are classified with a ~0.97 probability using Matlab's Neural Network Toolbox which is why I believe that the configuration could use some work and the data contain the necessary unique patterns.~~ Cheers. * * * **Edit to address the strikethrough:** The Matlab program in question was not mine and it turned out that the input data were terribly flawed during preprocessing so none of it was true. * * * **Edit to answer @lambruscoAcido questions:** I have learned a lot new information since asking this question and some parts of the procedure as well as results have changed but the general problem persists – choosing a correct model. Despite reading through loads of literature, I still lack basic understanding of certain concepts.   * The raw data coming from each subject are processed before feeding them to the NN (MLP) by extracting the same single channel, filtering out everything except beta range, splitting them into chunks of 256 samples with 50% overlap (fs = 128 Hz) and obtaining a PSD using Welch's method so they are skewed and look like this: ![Data](http://i.stack.imgur.com/3mnPr.png)   * They do correlate with each other but I am not sure about the principal components – how can I find out? I am interested in the general methodology of empirically establishing the optimal number of hidden layers, the number of perceptron units and their activation functions (I am currently using tanh for hidden layers and softmax for output). The studies I read did not include any information regarding this procedure.   * The dimensions of my current network are 129-10-6. The input is the power spectral density estimate (as pictured above) computed from each chunk of 256 samples times 30 per each of 6 subjects minus 1 due to overlapping (59 inputs per subject in total). I am currently achieving a mean error of ~15% after 100 epochs so there has been an improvement. The studies I have read were at ~2%.


title: NMDS: why is the r-squared for a factor variable so low

text: I am doing an NMDS ordination. The data come from a number of sites scattered around two lakes. In the plots, I coloured the samples from the two lakes blue and green. There seems to be some pretty decent separation between the lakes in the NMDS ordination plot. However, when I use lakes as a factor variable with envfit, I only get an R squared value of 0.19 (although the pvalue was less than 0.001). The R squared value is way lower than I thought. Am I missing something or am I reading too much into this? ![enter image description here](http://i.stack.imgur.com/C8Jie.jpg)


title: Differences and relation between retrospective power analysis and a posteriori power analysis?

text: From a note > **A Priori Power Analysis.** This is an important part of planning research. > You determine how many cases you will need to have a good chance of > detecting an effect of a specified size with the desired amount of power. > See my document Estimating the Sample Size Necessary to Have Enough Power > for required number of cases to have 80% for common designs. > > **A Posteriori Power Analysis.** Also know as “post hoc” power analysis. > Here you find how much power you would have if you had a specified number of > cases. Is it “a posteriori” only in the sense that you provide the number of > number of cases, as if you had already conducted the research. Like “a > priori” power analysis, it is best used in the planning of research – for > example, I am planning on obtaining data on 100 cases, and I want to know > whether or not would give me adequate power. > > **Retrospective Power Analysis.** Also known as “observed power.” There are > several types, but basically this involves answering the following question: > “If I were to repeat this research, using the same methods and the same > number of cases, and if the size of the effect in the population was exactly > the same as it was in the present sample, what would be the probability that > I would obtain significant results?” Many have demonstrated that this > question is foolish, that the answer tells us nothing of value, and that it > has led to much mischief. See this discussion from Edstat-L. I also > recommend that you read Hoenig and Heisey (The American Statistician, 2001, > 55, 19-24). A few key points: > >   * Some stat packs (SPSS) give you “observed power” even though it is > useless. >   * “Observed power” is perfectly correlated with the value of p – that is, > it provides absolutely no new information that you did not already have. >   * It is useless to conduct a power analysis AFTER the research has been > completed. What you should be doing is calculating confidence intervals for > effect sizes. > I was confused what differences and relation are between Retrospective Power Analysis and A Posteriori power analysis? I think they both have a given sample size and try to estimate the power? What does "the size of the effect in the population was exactly the same as it was in the present sample" mean? Thanks and regards!


title: Are the sequential sum of squares appropriate when treatments must be applied in sequence?

text: I'm working with some modeled future stream flow data that was created in two steps. First, future precipitation predictions (at certain points in a watershed) were created by running historical precipitation data through four precipitation models. The output from those models (four runs per precipitation model) were then used as inputs in four different hydrologic models, which were used to generate future stream flow predictions (at those same points). This was done in such a way that each combination of precipitation model and hyrdologic model were considered (i.e. crossed). I'm trying to determine the relative importance of both precipitation model and hydrologic model choice on modeled stream flow, and have analyzed the output using a two factor linear model. Unfortunately, due to some issues with two of the models, a number of data points had to be thrown out. As a result, the model is unbalanced. So when looking at, say, the sum of squares, the type I and type II SS differ. To note, both variables stay in the model during model selection. I don't mean to resurrect a sum of squares debate. But I'm curious, if in a situation like this, because the data was created by applying one treatment and then another, if it's more appropriate to apply sequential sum of squares than the type II sum of squares. Many thanks! Jon


title: Accuracy and Statistical Significance

text: According to experimental result, One approach has 10% better accuracy over the other. My supervisor said that this should be coupled with an estimative of the statistical significance of this increase (p-value) according to the size of the dataset. Do you guys have any idea how to calculate statistical significance in this context?


title: How is research based on the U.S. Census organized?

text: As I understand, most of the research is centralized through the U.S. Census Bureau. However, I couldn't find any overview of fields of research on their website. Is there any other institution or website that summarizes the research that has been done? I'd especially like to know which data is of most interest to researchers but everything seems to be so restricted and it is hard to find a collection of publications regarding research outcomes based on the census. I appreciate any advise or direction.


title: What does linear stand for in linear regression?

text: In R, if I write               lm(a ~ b + c + b*c)       would this still be a linear regression? How to do other kinds of regression in R? I would appreciate any recommendation for textbooks or tutorials?


title: Back-transformation of regression coefficients

text: I'm doing a linear regression with a transformed dependent variable. The following transformation was done so that the assumption of normality of residuals would hold. The untransformed dependant variable was negatively skewed, and the following transform made it close to normal: $$Y=\sqrt{50-Y_{orig}}$$ where $Y_{orig}$ is the dependent variable on the original scale. I think it makes sense to use some transformation on the $\beta$ coefficients to work our way back to the original scale. Using the following regression equation, $$Y=\sqrt{50-Y_{orig}}=\alpha+\beta \cdot X$$ and by fixing $X=0$, we have $$\alpha=\sqrt{50-Y_{orig}}=\sqrt{50-\alpha_{orig}}$$ And finally, $$\alpha_{orig}=50-\alpha^2$$ Using the same logic, I found $$\beta_{orig}=\alpha\space(\alpha-2\beta)+\beta^2+\alpha_{orig}-50$$ Now things work very well for a model with 1 or 2 predictors; the back- transformed coefficients resemble the original ones, only now I can trust the standard errors. The problem comes when including an interaction term, such as $$Y=\alpha+X_1\beta_{X_1}+X_2\beta_{X_2}+X_1X_2\beta_{X_1X_2}$$ Then the back-transformation for the $\beta$s are not so close to the ones from the original scale, and I'm not sure why that happens. I'm also unsure if the formula found for back-transforming a beta coefficient is usable as is for the 3rd $\beta$ (for the interaction term). Before going into crazy algebra, I thought I'd ask for advice...


title: Fitting logistic function with pymc

text: I've asked this question on stackoverflow too, but no answer yet. This seems a more appropariate place to ask this question: I'm messing around with pymc to understand it a bit better. Now I am trying to fit a logistic function to my generated data. Though the parameter it returns is not even close to the real parameter. Any idea why ?               import pymc as pm     from matplotlib import pyplot as plt     import numpy as np     from pymc.Matplot import plot as mcplot               #data generation     def logistic(x, beta):         return 1.0 / (1.0 + np.exp(beta * x))          #hidden     beta = -5          x = np.linspace(-4, 4, 600)     data = logistic(x, -5)          plt.plot(x,data , label=r"$\beta = -5$")     plt.legend();          plt.show()                         #Priors     beta = pm.Uniform("beta", -10, 10, value=1.)     sig = pm.Uniform("sig", 0.0, 100.0, value=1.)               @pm.deterministic     def logistic(x=data, beta=beta):         return 1.0 / (1. + np.exp(beta * x))               #likelihood     y = pm.Normal("obs", mu=logistic, tau=1.0/sig**2, value=data, observed=True)               model = pm.Model([y, data,beta,sig])     mcmc = pm.MCMC(model)      mcmc.sample(20000, 5000, 1)               beta_samples = mcmc.trace('beta')[:]     sig_samples = mcmc.trace('sig')[:]          print beta_samples.mean()     print sig_samples.mean()          mcplot(mcmc.trace('beta'))     mcplot(mcmc.trace('sig'))     plt.show()


title: Likelihood function of a sample selection model

text: I'm having some problems with the following assignment: I have a sample selection model (tobit II model) which however does not have the standard tobit II representation but rather the sample selection model is described by y=1(xδ+u>0)e^(xb+e), where 1(.) is an indicator function taking the value of 1 if the argument in parenthesis is true, δ,β are coefficients and u, e are the residuals. The indicator function represents the sample selection issue. The standard tobit II model(truncation) is of the form y=1(xδ+u>0)(xb+e), thus without the exponential function. I am essentially asked to find the log likelihood of this model where truncation takes place and i am stuck. Essentially i am trying to find the pdf's from both the truncated part and the non truncated part of the model but so far it hasn't worked out. If anyone could point me towards the right direction i'd really appreciate it!. Thanks in advance for any replies!


title: How to interpret a QQ-plot of p-values

text: I am doing GWAS SNP association studies on diseases by using a software called plink (http://pngu.mgh.harvard.edu/~purcell/plink/download.shtml). With association results I get p-values for all the SNPs that was analyzed. Now, I use a QQ-plot of those p-values to show if a very low p-value differs from the expected distribution of p-values (a uniform distribution). If a p-value deviates from the expected distribution one "may" call that p-value for statistic significant. As you can see in the QQ-plot, at the top tail end, the last 4 points are somewhat hard to interpret. Two of the last points in the grey suggests that those p-values are in the expected distribution of p-values, whilst the other two are not. Now, how to interpret this, the last two points have _lower_ p-values but are not "significant" according to the QQ-plot, whilst the other two points with _higher_ p-values are "significant"? How can this be true? ![enter image description here](http://i.stack.imgur.com/zwSwt.png)


title: Any known example of research project being led by statistician in public health

text: (I am fully aware that my topic may not be relevant here, do let me know if this is true.)   I am a research staff with training in epidemiology and biostatistics, working in government organization on public health matters like disease surveillance and infection control. Being a research staff here keeps giving me a feeling that I cannot lead a research project in this position, because I lack the context needed to ask sensible research questions. Doctors/PhDs in our service, although they know little about statistics, can become principal investigators of research projects, and then assign research staff to handle the technical parts.   I wonder if this kind of "statistician can only take supportive role" curse is universal to other fields also? Any comment or counter argument on this, and any suggestion of breaking this kind of curse?


title: Can not classify Iris dataset by continuous deep belief network

text: We have a case to map real input values into low dimensional binary vectors. So we use CDBN here, but whatever I change the parameter combinations of k, lr or epochs, it always failed to classify our dataset. So we change our implementation from our own code to deeplearning4j, but it fails also. Under such a situation, I wrote a simple iris dataset test to verify my understanding https://github.com/guokr/neuseg/blob/master/test/com/guokr/dbn/IrisTest.java It fails again and again under different params. for one example: under k = 50, lr = 0.001, and epochs = 1000, the output is [0.460480, 0.503034, 0.036486; 0.460480, 0.503034, 0.036486; 0.460480, 0.503034, 0.036486] but it should be [1, 0, 0; 0, 1, 0; 0, 0, 1] or some approximative values. So I am not sure whether my understanding is wrong or we choice a wrong model? Thanks for your help!


title: Given two sets of data, what could explain similar means but different standard deviations?

text: Given two sets of data of user activity, both of which appear to be in an exponential distribution, I have calculated the mean and standard deviations using both a mean/deviation and a sample mean/deviation `(sample size = 30, number of samples = 10k)`: * * * ### A `(size: 627,000)`:   * Raw -- ` μ = 45.947, σ = 114.2, σ/√n = 0.14422 ` ### B `(size:3570)`:   * Raw -- ` μ = 46.43, σ = 116.1 ` Using the above data, it seems that the two means differ by a statistically significant amount, and thus allow us to say with confidence that the average for B is greater than the average of A. * * * ### A   * Sampling -- ` μ = 46.174, σ = 21.256 ` ### B   * Sampling -- ` μ = 46.786, σ = 21.366 ` Using the sampling data standard deviations, we see that the difference in means (`0.612`) is much less than the deviation, it seems that the means do _not_ differ by a statistically significant amount. * * * So given the above, which is right? Can we say that these data sets differ? If the underlying distribution is exponential, are the above tests even accurate?


title: Forecasting time series with missing data and irregular intervals

text: I have a data set of medical drug stock levels at health centres and I want to forecast monthly consumption over the following 3-6 months. However about 30%-40% of the data is missing and some of the reporting intervals are irregular. Several other variables can potentially be used to further predict consumption. I'm looking for suggestions for how to best approach this problem. Some additional information about the data set:   * Total of ~4000 facilities   * Data spans 1 year for most facilities. A small subset has ~3 years of data.   * Facilities submit monthly reports for current stock on hand for up to 20 products.    * Stock receipts occur roughly every 3 months but only ~30% are explicitly reported - the others can be detected by an increase in the stock level   * Less than 1% of all product/facility combinations have 1 year of continuously reported data that can be used for testing    * Occasionally the data contains data points mid-month. These data points are particularly important when a facility reports 0 stock and consequently its consumption is 0.    * A stock-out frequently goes unreported until end of the month. This is hard to detect if combined with an unreported receipt of stock that same month.   * The consumption rate of different drugs is likely proportionally co-related for one facility: if a facility consumes 2 times the amount an average facility consumes for drug X, that facility will likely consume 2 times the average consumed of drug Y.   * The consumption rate is likely correlated to the size of the community it's serving (catchment population). This data is available for a subset of facilities.   * The consumption is likely correlated to the consumption of geographically proximate facilities. GPS coordinates for ~60% of facilities is available. All facilities are grouped into districts of 10-50 facilities of varying geographic size.  I'm assuming it's best to start with a simple model and slowly build up its complexity to include additional factors. What models and methods are most appropriate?


title: Sandwich estimator intuition

text: Wikipedia and the R sandwich package vignette give good information about the assumptions supporting OLS coefficient standard errors and the mathematical background of the sandwich estimators. I'm still not clear how the problem of residuals heteroscedasticity is addressed though, probably because I don't fully understand the standard OLS coefficients variance estimation in the first place. What is the intuition behind the sandwich estimator?


title: Calculating the n value

text: I have some practice homework questions. I did the first one I will go over the steps please tell me If I am doing it right. a) As mentioned earlier, it is claimed that 70% of households in Ontario now own large-screen TVs. You would like to verify this statement for your class in mass communications. If you want your estimate to be within 5 percentage points, with 99 per cent level of confidence, how large of a sample must be acquired? So for this i knew the targe parameter was p so the standard error is $\sqrt{\frac{pq}{n}}$ and that is equal to 0.05 and since we are finding the 99% C.I. our $Z_\alpha /2$ = 2.575 so my formula became $0.05 = 2.575 \cdot \sqrt{\frac{(0.7)(0.3)}{n}} \Rightarrow n = 556.07 = 557$ b) You are to conduct a sample survey to determine the mean annual family income in a rural area. The question is, “How many families should be sampled?” In a pilot study of just 10 families, the standard deviation of the sample was $\$500$. The sponsor of the survey wants to use the 95 per cent confidence level. The estimate is to be within $\$100$. How many families’ should be interviewed? For this I kinda got confused because it said how many family should be interviewed but we are given 10 families? but I was going to follow my procedure as in a but I didnt know where the n would come from because we are give the S.D. ? or would it be this $100 = 1.96 \cdot \frac{500}{\sqrt{n}}$ and solve for n afterwards? c) A student conducted a study and reported that an 80% confidence interval ranged from 48 to 52. He was sure the sample standard deviation was 16 and that the sample was at least 30, but could not remember the exact number (i.e., the size of the sample). Can you help him out? for this since we are finding the 80% CI our area is 0.9 and our critical value now would be 1.285. ans since we know it ranged from 48 to 52 our mean is $\frac{48+52}{2} = 50$ and then that is equal to $50 = 1.285 \cdot \frac{16}{\sqrt{n}}$ and then solve for n only confusion was is that number 50 right? I know this is wrong because my n value is less than 30, please help out Thanks alot for your help appreciate it