title: Generalization of skewness detection in non-uniform case

text: I am not sure if my title reflects what I am attempting to ask. Let me elaborate my question. I want to classify an observation/feature into one of two categories ($A$ and $B$) based on frequency of occurrence of feature in each category. For example, a particular observation occurs $x$ times in category A and $y$ times in category B, and I have same number of sample for each category. If $x/(x+y> c$ for some $c > 0.5$, e.g. $c = 2/3$, then I associate that feature with category A else, if $y/(x+y) > c$ then I associate it with category $B$, else discard the feature as having no discriminating information. Here I was assuming that I have same number of sample from each category. Now, if I had for example, $a$ sample from $A$ and $b$ sample from $B$ to check for the occurrence of features, how do I detect if my observation is skewed towards a category.


title: Difference between Sobol indices and total Sobol indices?

text: Given a mathematical model $Y\widetilde Y(X_i)$, where $X_i=x_i^*$ represents a particular point estimate for input variable $X_i$. In sensitivity analysis, Sobol indices explain the importance of an input factor $X_i$ on the variance of the output $Y$ such that: \begin{equation} S_i=\displaystyle \dfrac{V_{X_{x \sim i}}(E_{X_i}(Y|X_i=x_i^*))}{V(Y)},\quad 0 \leq S_i \leq 1 \end{equation} This represents the variance of the expectation of $X$ in the red box over the total variance of $Y$. ![enter image description here](http://i.stack.imgur.com/wGUNW.png) Total Sobol indices are _seemingly_ given by this Wikipedia entry: $$S_{Ti}=\displaystyle 1-\dfrac{V_{X_{x \sim i}}(E_{X_i}(Y|X_i=x_i^*))}{V(Y)}$$ But what do they really mean? I'm not convinced.


title: Spatial autocorrelation (SAC) while analysing survey data

text: I am confused about some aspects of spatial autocorrelation usind survey data (survey which is repeated every year). I have data from 1991 to 2012 with sampling region pretty consistent every year. I am studying the species distribtion of two species in the area, in this goal I used binomial GLM. After running the GLM, I noted that they were some spatial autocorrelation in the model error residuals (Moran's I index from 0.2 to 0.05 significant on 6 first lags). Therefore I decided to take in account the spatial autocorrelation using Moran's Eingenvector Mapping. However this method is generally described for dataset over 1 year. As the eigenvector are calculated using the coordinates (longitude and latitude) of my points in function of the neighbourhood response I assumed it should take in account the different distribution over the years. Thus I applied a 10 factor on longitude and lagitude to artificially separate points at the same location but sampled at different years (point A year 1 will have (1,1) coordinate; point B year 2 will have (10,10) coordinate. My first question is to know if this is valid? If not what would you advise? The outcomes of this method is that I have different eigenvector value for a point at the same position but at different years. Is this not a problem for the model including MEM variables power of prediction? Indeed following a question not answered yet I posted here (Spatial autocorrelation -- GLM, autocovariate, MEM (Moran's eigenvector mapping)), it appears that my GLM including MEM suffer a loss of 0.25 or more power of prediction (ROC Curve - AUC value). I am digging since one week to understand what is happening and why spatial autocorrelation included using MEM method are so unefficient while they seems to be the best method following litteratur. Any clue or help will be really appreciated, Thank you, Xochitl C.


title: Election forensics using statistical methods in practice?

text: Are you aware of any examples of election forensics in practice? Or at least any applied research on real large-scale datasets (e.g. govermental elections)? Thanks.


title: Predicting raffle winners

text: I don't have a strong background in statistics, so forgive me if I'm making very elementary mistakes here or if this question is too simplistic. I am in a club with ~80 members, and every week, 10 members are selected to win a prize. No member can be selected to win a prize twice in a row, but other than that, anything goes. Some members are selected more often than others, though. For example, a few people are consistently selected as a winner pretty much every second or third week, while several others get picked only a few times a year. There's also a bunch of new members, so there's not enough data to know how often they will get picked (though they automatically get picked the week that they join). I would like to make a program that would give a rough prediction of the chance that each member gets picked in the next week, and in the next month, based on previous data (I have records from every week for the past year). Right now, my idea is to look at the gaps in between weeks that each member gets picked (e.g. "Annie"'s gaps are 3 3 3 3 3 2 4 3 4 3 5 7 4, so she usually gets picked at least once a month), and assume that this is a normal distribution. Then I look at how long has it been since she was last picked, to see what would be the next gap would be if she gets picked next week, and then look at the standard bell curve to determine the % chance that the next gap would be that value or lower (one-sided z-score). The problem that I see is that this doesn't take into account the "no member will get picked twice in a row" rule (which may be significant for those who get picked very very often), and I am not sure how well the standard normal curve works for discrete values like number of weeks (as opposed to something like height, which can be). Is my idea good enough to get at least something reasonably useful, or are there just too many issues with it? If so, is there an alternate model that I can use? Again, I'm not looking for anything rigorous or overly complicated (like Hidden Markov models and such) that takes everything into account perfectly, just something remotely useful to get me a rough idea of the upcoming picks. EDIT:   Yes, there is some sort of mechanism where certain members consistently get a better chance of winning, but we don't know the exact details. It can be approximated by the analogy of "some people having more raffle tickets than others, but nobody knows how many tickets each member has in the raffle". The only information we have on the "weight" given to each member's chances is the past record of winnings. In actuality, there's a lot of different factors that influence who will win the raffle for that week, but almost none of those details are known for sure. I am assuming that these factors remain constant over time, though in actuality this is not always the case. For example, a member could get banned from winning anything for several weeks if they broke certain rules and such... But it'd be impossible to come up with anything usable if we need to take that into account.


title: Standardized residuals vs. regular residuals

text: I've got an easy question concerning residual analysis. So when I compute a QQ-Plot with standardized residuals $\widehat{d}$ on the y-axis and I observe normal distributed standardized residuals, why can I assume that the error term $u$ is normal distributed? I'd think that if $\widehat{d}$ looks normal distributed I just could assume that the standardized error term $d$ should be normal distributed. So why can we assume that $u\sim N$ when we just observe that $\widehat{d}\approx N$. By the way can we assume it?


title: Why is the CI for an odds ratio not always centered on the sample value?

text: Also, what will tend to produce a CI in which the sample value is far above the center, or far below the center?


title: Analysis of covariance in R

text: 1. How do you carry out analysis of covariance using R?    2. How do you interpret the results?    3. A practical example will be highly appreciated.


title: Formal statistical test for whether a process is a white noise

text: Is there a formal statistical test to test if process is a white noise?


title: Differencing weekend fluctuations with R?

text: Suppose a time-serie like this on the left-top corner with weekend and daily fluctuations. This time-series need differencing due to the rising ACF (bottom-left) and portmanteau tests' p -values too small, the weekend must be differenced away. Now I want to differentiate away the weekend-fluctuation in R, how to do that? This puzzle is from the lab here (sorry not in English), data here. ![aika/vrk = time/day, kulutus=consumption, ](http://i.stack.imgur.com/ARH6x.png) _(translation: aika/vrk = time/day, kulutus=consumption)_ > **My ideas (silent discussionhere):** > > A) Friday-Monday -average for the weekend fluctuating values > > B) create some sort of PID -controller and differentiate with that info > > C) ? > > **Possible tools in R** > >   1. `nls()`-`lm()`-`diff()` -problem? > >   2. `diffinv()`? > > > > **Some References** > >   1. Load Forecasting -paper by Eugene A. Feinberg: > www.almozg.narod.ru/bible/lf.pdf >


title: What is a white noise process?

text: What is the best way of defining white noise process so it is intuitive and easy to understand?


title: Groups in linear regression with different intercepts. How do I find the differing variable?

text: This is more of a conceptual question. I have a coefficient estimate of .80 in a linear regression model with one IV and one dependent variable. However, plotting the data I see distinct groups, differing by the intercept. I have about 300 additional variables in the data set that could have (or could not have; maybe it's something we didn't ask) an impact on the different intercept. I want to find the variables that explain the difference between these groups. What I thought about so far:   * I won't find it with a moderation model because the difference is in the intercept, not the slope   * I already tried separating the groups and conducted `t.test()` on all of the 300 additional variables, but isn't there a more formal way to find the differing variable(s)?  I hope I could adequately describe the problem... I include a plot of the two variables. Can someone give me a hint? > ![enter image description here](http://i.stack.imgur.com/YlJXO.png)


title: Random assignment: why bother?

text: Random assignment is valuable because it ensures independence of treatment from potential outcomes. That is how it leads to unbiased estimates of the average treatment effect. But other assignment schemes can also systematically ensure independence of treatment from potential outcomes. So why do we need random assignment? Put another way, what is the advantage of random assignment over nonrandom assignment schemes that also lead to unbiased inference? Let $\mathbf{Z}$ be a vector of treatment assignments in which each element is 0 (unit not assigned to treatment) or 1 (unit assigned to treatment). In a JASA article, Angrist, Imbens, and Rubin (1996, 446-47) say that treatment assignment $Z_i$ is random if $\Pr(\mathbf{Z} = \mathbf{c}) = \Pr(\mathbf{Z} = \mathbf{c'})$ for all $\mathbf{c}$ and $\mathbf{c'}$ such that $\iota^T\mathbf{c} = \iota^T\mathbf{c'}$, where $\iota$ is a column vector with all elements equal to 1. In words, the claim is that assignment $Z_i$ is random if any vector of assignments that includes $m$ assignments to treatment is as likely as any other vector that includes $m$ assignments to treatment. But, to ensure independence of potential outcomes from treatment assignment, it suffices to ensure that each unit in the study has equal probability of assignment to treatment. And that can easily occur even if most treatment assignment vectors have _zero_ probability of being selected. That is, it can occur even under nonrandom assignment. Here is an example. We want to run an experiment with four units in which exactly two are treated. There are six possible assignment vectors:   1. 1100   2. 1010   3. 1001   4. 0110   5. 0101   6. 0011 where the first digit in each number indicates whether the first unit was treated, the second digit indicates whether the second unit was treated, and so on. Suppose that we run an experiment in which we exclude the possibility of assignment vectors 3 and 4, but in which each of the other vectors has equal (25%) chance of being chosen. This scheme is not random assignment in the AIR sense. But in expectation, it leads to an unbiased estimate of the average treatment effect. And that is no accident. Any assignment scheme that gives subjects equal probability of assignment to treatment will permit unbiased estimation of the ATE. So: why do we need random assignment in the AIR sense? My argument is rooted in randomization inference; if one thinks instead in terms of model-based inference, does the AIR definition seem more defensible?


title: Error distribution for linear and logistic regression

text: With continuous data, a linear regression $Y=\beta_1+\beta_2X_2+u$ assumes that the error term is distributed N(0,$\sigma^2$) 1) Do we assume that Var(Y|x) is likewise ~N(0,$\sigma^2$)? 2) What is this error distribution in logistic regression? When the data is in the form of 1 record per case, where the "Y" is 1 or 0, is the error term distributed Bernoulli (i.e. variance is p(1-p) )) and when the data is in the form #successes out of #of trials, is it assumed binimial (i.e. variance is np(1-p)), where p is the probability that Y is 1?


title: Iterative proportional fitting and independent variables

text: I'm trying to understand the "classic" iterative proportional fitting (IPF) algorithm. Does it always assume that the variables being analyzed are independent? If the variables are independent, then can't we just compute each joint count from the product of two marginal counts? For example, in the wikipedia example referenced above, are the two variables (gender and handedness) independent? The final joint counts are               45.24 = 100 * 52/100 * 87/100     41.76 = 100 * 48/100 * 87/100     6.76 = 100 * 52/100 * 13/100     6.24 = 100 * 48/100 * 13/100      so no iterative algorithm is necessary.


title: Is it appropriate to plot the mean in a histogram?

text: Is it "okay" to add a vertical line to a histogram to visualize the mean value? It seems okay to me, but I've never seen this in textbooks and the likes, so I'm wondering if there's some sort of convention not to do that? The graph is for a term paper, I just want to make sure I don't accidentally break some super important unspoken stats rule. :)


title: Strange behavior of predictive model gain curve

text: When building predictive models (binary target), one of the principal methods I use for determining how useful the model is, is to plot the true proportion of y=1 values for each decile of the probability score. [These models are marketing models involving rare events and accuracy or the like are not suitable methods for validation]. For example, I build a model on a training set and use the model to score the hold out test set. Then   * Rank the scored data descending by the prediction (probability the target variable is 1)   * Decile (10 equal sized groups) the scored data by this descending prediction (1=highest , 10= lowest)   * Plot the actual proportion of 1's for each decile. I see from time to time a pattern like the following, which has good ranking discriminatory power but is not strictly monotonically decreasing. Instead, there will be either   1. A hump in the middle OR   2. Misaligned top and bottom of the curve. Here are fictitious examples of case #1 and case #2, followed by the "ideal". **Question:** What causes this type of result - in other words is this a diagnostic for some deficiency in the model that suggest a coarse of corrective action to achieve the ideal? It seems the first shows the model can separate really good and really poor but not the "mushy middle". The second case is less clear to me. I am not sure a correction for either. ![enter image description here](http://i.stack.imgur.com/oTX6b.jpg)


title: How many cases are needed for an extrapolation?

text: In social science it is comon to use sample weight to post adjust the sample so that it fits to a given basci population. There are algorithms to calculate how many cases you need to have a representive sample of the population. But often we would like to analyse a bit deeper than just the sample as a whole e.g. male vs. female or age groups. Theoretically - if I'm not wrong - we need for each field a representive number of cases. So if I need 360 cases in the whole sample, I need for the analysis of male vs. female 360 * 2 (one for each field) If we use weights, the question is on who many cases we the extrapolation is based. But how many cases I need at least to extrapolate my sample? I have the following table with observed cases per field:                                      Male                 ||            Female                     Education                   ||      Education        age   |  low  | middle| high  |  sum     ||  low  | middle|  high |  sum  |  SUM     ---------+-------+-------+-------+----------++-------+-------+-------+-------+------     18 - 24  |   26  |   60  |   23  |  109     ||   33  |   45  |   20  |   98  |  207     25 - 39  |    6  |   78  |   94  |  178     ||    8  |  114  |  147  |  269  |  447     40 - 49  |    1  |  105  |  134  |  240     ||    6  |  124  |  172  |  302  |  542     50 - 59  |    6  |   63  |  117  |  186     ||    8  |  117  |  146  |  271  |  457     60 - 69  |    3  |   48  |  110  |  161     ||   17  |   99  |  122  |  238  |  399     70+      |    7  |   30  |   75  |  112     ||   41  |   89  |   58  |  188  |  300     ---------+-------+-------+-------+----------++-------+-------+-------+-------+-------     SUM      |   49  |  384  |  553  |  986     ||  113  |  588  |  665  | 1366  | 2352      If we take a look at the males with a low education in the different agegroups, I have doubt, that we can extrapolate them. But I can't argue that right now. Even if we take just the low educated males at all, I don't think that the amount of cases are enough to extrapolate them. What I have to do is to analyse injuries in sex, education leavel and age groups. But because of this distribution table I don't know, if the reslaults are representive for the basic population wich is over all about 2 million people. So here my questions:   * How I can figure out if a number of cases is big enough to get extrapolated   * Do I really need, as I mentioned above, the amount of samples per field e.g. 360 males in the age of 18-24 with a low education level to analyse this group?   * On which categories (e.g. age groups, sex or combined) I can run representive calculations? Maybe there is also something with which I can do the extrapolation with R and something like a test how good my extrapolation is... Thanks for y'all help. Dominik


title: What is the correlation of [1,2,3] and [1,5,7] to 8 decimal digits?

text: This stackoverflow post describes computing a Pearson correlation of [1,2,3] and [1,5,7] in several different ways in Python. The most straightforward implementation from the definition in Wikipedia comes up with               0.973328526785      while Excel, R, NumPy, an online calculator, and a different Python implementation (involving what looks like a more numerically unstable calculation to me) come up with               0.981980506      I am just curious to know what you think.


title: How to draw a random sample from distribution of prediction?

text: For my microsimulation, I want to use R to predict values and draw a random sample based on this prediction. To clarify my point: I want to simulate the number of chronic conditions people suffer from ($y_t$) at a certain point in time. I have a few waves of panel data available to estimate a relation between age, sex, number of chronic conditions in the previous observation period (plus some others that I might include in later stages). Suppose my regression model is $y_t = β_0 + β_1 age + β_2 sex + β_3y_{t-1} +u $. Since R provides me with coefficients for the betas, it is easy to predict $y_t$ given the independent variables. However, this is not what I want to do. Instead, I want my population of about 1300 individuals to resemble the variance in the possible outcomes of $y_t$ (otherwise after a few steps my simulated population won’t include those unlucky ones with much more chronic conditions than the average). I believe what I have to do is to draw a random sample from the distribution of the predicted value $y_t$, conditional on the independent variables. I further believe this can be done by drawing random numbers with mean $β$ and variance $var(β)$, multiplied by the actual values of the independent variables. So my question is: Is this the correct approach? Will this produce reliable values? Or do I need to take possible covariation of the independent variables etc. into consideration? **Edit** : Another point came to my mind. Does it make a difference whether $y_t$ or $y_t - y_{t-1}$ is my left hand side variable? Thank you for your ideas.


title: Given P(A) and P(B), what would be the minimum probability of the intersection?

text: Like if I was given a P(A) of .5 and a P(B) value of .4 how would I get the minimum of the P(A∩B)?


title: Lognormal distribution using binned or grouped data

text: I understand the Max likelihood estimators for mu and sigma for the lognormal distribution when data are actual values. However I need to understand how these formulas are modified when data are already grouped or binned (and actual values are not available). Specifically, for mu, the mle estimator is the sum of the logs of each X (divided by n which is the number of points). For sigma squared, the mle estimator is the sum of (each log X minus the mu, squared); all divided by n. (Order of operations is taking each log X minus the mu; square that; sum that over all X's; then divide by n). Now suppose data in bins b1, b2, b3, and so on where b1 to b2 is the first bin; b2 to b3 second bin and so on. What are the modified mu and sigma squared? thank you.


title: Testing two I(1) vectors for a relationship

text: Suppose I have two I(1) time series X and Y, and I want to know whether X and Y are "related" (for some definition of "related"). The standard cointegration approach defines relationship as cointegration, and says that X and Y are cointegrated if some linear combination of X and Y is stationary. To test whether X and Y are cointegrated, you perform a regression on X and Y, and test for stationarity of the residual errors. It seems to me like another approach might be to difference the I(1) time series X and Y, to get new I(0) time series X' and Y', and to use a standard linear regression relationship test on X' and Y' (i.e., perform a regression to get Y' = aX' + b, and use a t-test to see whether a is significantly non- zero). You could then define X and Y to be related if X' and Y' pass this test. Is this second approach valid, or do you get spurious relationships? What's the difference between this approach and the cointegration approach, or what are the advantages of the cointegration definition?


title: How to create one score from a mixed set of positive and negative variables?

text: I have 3,000 observations (administrative communities) characterized by five variables. Four of them work in the direction 'the more, the worse' and one goes in the opposite. I'd like to create one score or an ordered list of these observations that will best take into account all of those five variables. I have tried clustering using MCLUST package in R, and it gives some meaningful results but it's hard to decide about the ordering of observations on the basis of cluster membership. My second attempt was to run PCA and extract the first component, which is closer to what I'd like to get. What other solutions (R- or Stata-based preferably) could I use to deal with this problem?


title: How to identify invalid online survey responses?

text: I have a set of around 500 responses to a online survey that offered an incentive to complete. While most of the data appears to be valid, it's clear that some people were able to get around the (inadequate) browser cookie-based duplicate survey protection. Some respondents clearly randomly clicked through the survey to recieve the incentive and then repeated the process via a couple methods. My question is what is the best way to try and filter out the invalid responses? The information that I have is limited to:   * The amount of time it took to complete the survey (time started and ended)   * The IP address of each respondent   * The User Agent (browser identifier) of each respondent   * The survey answers of each respondent (over 100 questions in the survey) The most obvious sign of an invalid response is when (sorted by time started) there will be a group all from the same IP address or similar IP (sharing the same first three octets for example, 255.255.255.*) which were all completed in a much shorter amount of time than the total average in quick succession. With this information there must be a thoughtful way to weed out the people who were exploiting the survey for the incentive from the rest of the survey population. I know that someone from the community here would have an interesting idea about how to approach this. I'm willing to accept false positives as long as I can be confident that I've gotten rid of most of the invalid responses. Thanks for your advice!


title: What is the variance of the sum of components of a multivariate normal distribution?

text: I'm talking with my advisor about how to compute standard deviations for, say, combined standardized test scores for admissions purposes. For example, we'd be interested to compute the sum of the verbal and quantitative scores from the GRE, which are correlated, and normed to be approximately normal. More formally, say you have a multivariate normal with vector mean and matrix covariance $X \sim N(\mu, \Sigma)$, with $X = (X_1, X_2, ...)$ (and covariances are non-zero). What is the variance of $\sum{X_i}$? If it's hard to compute in general, I'm happy with bivariate for now, a recursive approach or similar.


title: Estimating the parameters of a sum of a Gaussian and an $\alpha$-stable random variable

text: Let's assume I have a set of samples of a random variable $$ X = Y + Z \>, $$ where $Y$ is Gaussian (with a mean of zero and variance $\sigma^2$) and $Z$ has a symmetric $\alpha$-stable distribution with index of stability $\alpha$ and scale parameter $c$. I know $\alpha$, but I don't know $\sigma$ and $c$. How do I go about estimating them?


title: Function timings in the Matrix package

text: In order to test the Matrix package in R, I ran the following piece of code:               library(Matrix)          p <- 50; n <- 500     a0 <- sample(1:(n*p), n*p*0.1)     d1 <- (a0-1) %% n     d2 <- floor((a0-d1)/n)     d0 <- cbind(d1+1, d2+1)     x0 <- Matrix(0, n, p, forceCheck=TRUE)     x0[d0] <- rnorm(n*p*0.1, 10, 1)     x1 <- matrix(rnorm(n*p), n, p)          fx01 <- function(ll, x0) tcrossprod(x0, rnorm(ncol(x0)))     fx02 <- function(ll, x1) crossprod(t(x1), rnorm(ncol(x1)))          system.time(lapply(1:1000, fx01, x0=x0))     system.time(lapply(1:1000, fx02, x1=x1))      Is it normal that the first function takes 4-5 times more time to run than the second when the second does not use a sparse matrix ?


title: Outlier detection for generic time series

text: In this case, "generic" being the entire gauntlet of macroeconomic time-series that private and government statistical offices put out. Some background - I recently started working at a data provider - we collect data releases and repackage them in a presumably more convenient and accessible fashion for our clients, and we have tens of thousands of data series (wouldn't be surprised if we topped a million, actually). As part of our QA process, we run the following outlier detection: $X_t-X_{t-1} = E_t$   $\sigma^2$ is estimated from the resulting sample of $E_t$, and a z-score is calculated based off $E_t\sim N(0,\sigma^2)$ I think we can do better - the math clearly falls apart for everything that isn't a random walk. I initially thought of fitting an ARMA(m,n) based on the peak of the autocorrelation/autocovarience functions of the series and checking the residuals. I'm wary of the robustness of this, and a previous question seems to indicate that autocorrelation is not particularly robust.


title: Bayes classifier

text: Given the prior probability of 2 distributions, $N(x,y)$ and $N(a,b)$, where $N(\mu,\sigma^2)$: How do you make a decision rule to minimize the probability of error, if the prior probabilities are equal? Can you give an example? What if the prior probabilities are different, such as   $P(\text{Distribution 1}) = 0.70$   $P(\text{Distribution 2}) = 0.30$?


title: Statistical project directory structure with multiple languages (e.g., R and Splus)?

text: Building on the post How to efficiently manage a statistical analysis project and the `ProjectTemplate` package in R... Q: How do you build your statistical project directory structure when multiple languages feature heavily (e.g, R AND Splus)? Most of the discussions on this topic have been limited to projects which primarily use one language. I'm concerned with how to minimize sloppiness, confusion, and breakage, when using multiple languages. I've included below my current project structure and methods for doing things. An alternative might be to separate code so that I have `./R` and `./Splus` directories---each containing their own `/lib`, `/src`, `/util`, `/tests`, and `/munge` directories. Q: Which approach would be closest to "best practices" (if any exist)?   * **/data** - data shared across projects    * **/libraries** - scripts shared across projects    * **/projects/myproject** - my working directory. Currently, if I use multiple languages they share this location as their working directory.   * **./data/** - data specific to `/myproject` and symlinks to data in `/data`   * **./cache/** - cached workspaces (e.g., `.RData` files saved using `save.image()` in R or `.sdd` files saved using `data.dump()` in Splus)    * **./lib/** - main project files. Same across all projects. An R project will be run via `source("./lib/main.R")` which in turn runs `load.R`, `clean.R`, `test.R`, `analyze.R`, `.report.R`. Currently, if multiple languages are being used, say, Splus in addition to R, I'll throw `main.ssc`, `clean.ssc`, etc. into this directory as well. Not sure I like this though.   * **./src/** - project-specific functions. Collected one function per file.    * **./util/** - general functions eventually to be packaged. Collected one function per file.   * **./tests/** - files for running test cases. Used by `./lib/test.R`   * **./munge/** - files for cleaning data. Used by `./lib/clean.R`   * **./figures/** - tables and figure output from `./lib/report.R` to be used in the final report    * **./report/** - `.tex` files and symlinks to files in `./figures`   * **./presentation/** - `.tex` files for presentations (usually the `Beamer` class)    * **./temp/** - location for temporary scripts    * **./README**   * **./TODO**   * **./.RData** - for storing R project workspaces    * **./.Data/** - for storing S project workspaces


title: Is E(u|x)=0 is a required condition for estimator consistency?

text: For OLS parameter estimates to be consistent it must be the case that E(u|x)=0. Is it true? E(u|x)=0 is a required condition for unbiasedness. But as far as I understand, unbiasedness does not necessarily mean consistency. Therefore I am really confused.


title: Automatic text quality grading

text: I came up with an idea to use machine learning for automatic grading of topic- specific texts. More specifically, I will first use normal text classification techniques to sort all candidate texts into topics. Then, I want to be able to **judge the quality of texts in specific topics** : For instance, news articles in different topics (Technology, Sports, International, Business, etc). Firstly, every article will be sorted into one topic. And what I want to archive then is to automatically grade this 'Technology' article based on the range of 1-10 or good/poor. The criteria for such scoring scheme may be:   * Elaboration (rich in details)   * Coverage of all key words   * Length? (It is actually a little bit hard for me to quantify my criteria. The grading criteria here can be very similar to those used when human beings instinctively and quickly judge quality of two articles on a same story.) However, as I reckon, text classification based on 'bag of words' may not be able to perform this task well. (Or will it?) I believe such task is well-researched, and even it is a whole field rather than text classification. But I haven't figured out how to do this by searching around. Could you please refer me to related techniques/discussions/names or give me some hints if text classification can handle that. Many thanks in advance.


title: Probabilities in a Pick 3 lottery

text: In a Pick 3 lottery three numbers are drawn from three separate sets of numbers, between 0 and 9. Matching 2 winning numbers, one of them duplicated, can be done in 3 ways, according to all lotteries websites. But I found 18 ways to do it. As an example: If the winning numbers are 1 2 3, then there are 18 possible ways to win:               1   1   2     1   1   3     1   2   1     1   2   2     1   3   1     1   3   3     2   1   1     2   1   2     2   2   1     2   2   3     2   3   2     2   3   3     3   1   1     3   1   3     3   2   2     3   2   3     3   3   1     3   3   2      Since the probability is equal to the number of favourable outcomes divided by the number of possible outcomes, it will be 18 divided by 1,000. The odds would be 1 in 55.55 Please let me know if I am making a mistake in my calculations.


I wish to build Feed-forward GMDH-type Neural Network model in R. The most popular base function used in GMDH is the gradually complicated Kolmogorov- Gabor polynomial as explained in wikipedia. Any guidance with package/code etc. would be great. Regards, purnendumaity


title: Increase Recall Rate for SIFT

text: Is there a better way to increase **Recall Rate** when using SIFT features?   I am thinking a way to replace the NN1/NN2 ratio to account for slightly distorted objects.   Moving towards clustering and using BOW(Bag Of Words) seems a way but I need to do one-on-one match of objects in images rather than training and learning. This refrains me from thinking towards BOW. Anybody got any idea?


title: Sample variance

text: Given the following data: $102, 40, 27, 108, 124, 113, 143, 100, 115, 128$ If $\sum_{i=1}^{10} X^2 = 112600$, what is the sample variance? The correct answer is $1400$. However, when I tried calculating in Excel, I got an answer of $34,047,397$. I took each $X$ value, squared it, subtracted the mean of $112,600$, took the difference and squared that, summed it, and then divided by $10$ terms. I also tried using a chi-squared variance of $2v=20$, but that didn't match the answer either. Can someone please point me in the right direction?


title: Improving Python Exponential smoothing

text: I am going to improve my code to the Exponential smoothing I submitted to Statsmodel which can be found here. The code handles 15 different variation Standard Exponential Smoothing models including Holt-Winters, SES, Brown, Holt, and Damped models. I have several questions what are State Space Models and Error correction forms? I do not really understand them. Secondly how do I solve for the hyperparameters and the damped parameters? is there a good way to fit them? Third there are `for` loops I am trying to eliminate to reduce the speed. How do I do this? Also I want suggestions for summary and other features that can be easily added.


title: Small Sample Size in Kernel Density Estimation

text: I am working on a problem where I have to solve an optimization problem over a dataset of 6 variables (~300 data points per variable). The data set is a historical data set, unfortunately it is small and has fat tails. I could probably get some more points but not as many to make to change the size by an order of magnitude. The data is economic data (asset class data) and captures real life scenarios that have happened. Most parametric models I've tried (for e.g. t-distribution) fail at describing the data especially the joint movement. The data set is small and historical, it's not fully representative of the actual distribution but it's all I have, I was leaning towards a bayesian approach that would allow me to impose the uncertainty associated with the data. But someone suggested using a non-parametric approach to estimate the joint density and re-sample a large number of points from it and then doing the optimization. My concern is that I would be overfitting to this particular dataset by this technique because I'm not accounting for the fact that my sample size is small and that my future realization could be very different. May be the bayesian approach is not going to yield very different results. I was wondering what suggestions people who have dealt with a similar situation have. Is my concern even a valid concern ? Any questions (It's possible I've not described the problem well enough)/thoughts/comments/references would be appreciated.


title: MCMC methods - burning samples?

text: In MCMC methods, I keep reading about `burn-in` time or the number of samples to `"burn"`. What is this exactly, and why is it needed? ### Update: Once MCMC stabilizes, does it remain stable? How is the notion of `burn-in` time related to that of mixing time?


title: Alternative to Chow test for equality of estimates, when the data for the 1st group is not available

text: My goal is to compare published regression weights estimated from a sample in another country, with those of the same model established on our data.There is ample evidence on the web that a Chow test is the best method to test equality of regression parameters, by combining two groups. But a problem emerges when the data for the 1st study is not at hand. Alternatively, z-test for comparison of regression slope has been suggested. but as addressed in a similar question, this test is only robust when the ratio of two SEs for the coefficients does not exceed 3. Unfortunately this ratio of SEs is about 10 in case of my study. Again in the question attached, Welch-Satterthwaite t-test or Cochran test (Cochran and Cox 1950) is recommended for Behrens-Fisher problem. Since I don't have the data for the first group I only can perform the tests by hand, but I could not find the formula for Welch t test or CBF for comparison of slopes. The only Welch or Cochran t test I found is for comparison of the means, not the slopes. So four questions arise here:   1. Are these two tests the same?    2. Are there any alternatives in such a situation?    3. Or can I compare the R squared in the same model across the groups instead?    4. If yes, what test can be applied?


title: Single ARMA model for multiple time series

text: I have 365 days of hourly data (24points each) of a prediction error (realised -pred_day_before). I want to model the evolution of the prediction error as an ARMA process. Matlab System Identification Toolbox offers an ARMA estimation PER time series, whilst I have 365 time series (with supposedly the same characteristics). How to estimate just one ARMA model for all of the 365 rows of 24 values?


title: Admissibility and domination for estimators

text: Watching a video by the "mathematicalmonk" on the web, I was wondering how to answer this kind of questions: Given $X_1,\ldots,X_n\sim \mathcal{N}\left(\mu,\sigma^2\right)$. Assume that $\mu$ is unknown. Using the square loss $\mathcal{L}\left(a,b\right)=\left(a-b\right)^2$: 1) Does $\hat{\sigma}^2=\frac{1}{n}\sum_{i=1}^n(X_i-\bar{X})^2$ dominate* $s^2=\frac{n}{n-1}\hat{\sigma}^2$? Or viceversa? 2) Are $\hat{\sigma}^2$ and $s^2$ admissible**? 3) Is there a "best" $s_c^2=c\hat{\sigma}^2$ for some $c\geq 0$? *given decision rules $a,b$, $a$ dominates $b$ if $\mathcal{R}(\theta,a)\leq\mathcal{R}(\theta,b),\,\, \forall \theta \in\Theta$ and $\mathcal{R}(\theta,a)<\mathcal{R}(\theta,b),\,\,$ for some $\theta \in\Theta$, where $\mathcal{R}(\theta,f)=\mathbb{E}\left[\mathcal{L}\left(\theta,f\left(D\right)\right)|\theta\right]$, and $D$ is the sample. **A decision rule $a$ is inadmissible if there is another decision rule which dominates $a$. Otherwise, $a$ is admissible.


title: Meta-analysis of means and medians in R?

text: I've been asked to review the literature for all studies that estimate the mean and median Blood Lead Level (BLL) in a certain country, and then perform a meta-analysis to come up with one overall value for the mean and median level. Is such a thing typically done? I'm not interested in any "effect size," just a univariate measure. I'd like to create a forest plot and do some of the standard tests of heterogeneity across studies. Is there a package and function in R that will allow me to do a meta-analysis of mean and median BLLs? Just to clarify, this would be a separate meta-analysis for means, and a separate meta-analysis for medians.


title: Recommendations on methods to predict sales data

text: I'm very inexperienced in statistical analysis. I think my question is fairly basic so I hope I can get some good pointers so I will be able to learn more efficiently. I have extensive hourly sales data of a website going back years and I would like to make predictions on what the most important metrics (most likely sessions and order total) should be in the next few hours. The available data covers almost everything about transactions, how much money were spent, how many sessions, how many visits, etc, with an hourly granularity. Obviously these data follow yearly patterns, seasonal patterns, and daily patterns. i.e. growth is trending up year to year, holiday seasons are more busy, and it's more quiet during the night and busy during the day. What are some established statistical model that can be used for the hourly prediction? The predictions only need to be made for the next few hours, but should be fairly accurate; they should be almost always correct within some margins unless major abnormalities happen. I have looked into ARMA/ARIMA methods but I'm unclear if what I want to do is the best fit, since I have more information about my time series. For example, there is clearly a 24 hour cycle in the data. Please let me know what things I should be looking into. Thank you very much!


title: how to deal with autocorrelation

text: My study is related to determinants of corporate liquidity and I am working on eviews. My model gives durbin watson stat value 0.89. after applying fixed effect estimation the value changes to 1.47. should I still worry about autocorrelation and go for another remedy? if I change my model into log- linear form, the D-W stat value changes to 1.71. Intercept and one of the independent variables have also become insignificant. which remedy will be preferable? I need help to get rid of the autocorrelation disease


title: When a random variable has a distribution whose parameter is another random variable

text: Is there a standard name for a situation where a random variable follows a distribution whose parameter is another random variable ? For example a binomial(15,p) variable where the the p is distributed as beta(1,2), or a Poisson(Y) where Y is distributed as exponential(2) Is this called a compound distribution, or ? Then my real question is, given Y is distributed according to some given pdf with parameter X (say `pdf1`), but X is distributed according to another distribution (say `pdf2`), how do I use Bayes rule: $$ f_{X|Y}(x|y)=\frac{f_{Y|X}(y|x) \, f_X(x)}{f_Y(y)} $$ ? $f_X(x)$ must just be `pdf2`, right ? Is $f_{Y|X}(y|x)$ just the pdf of Y (that is, `pdf1`) with the pdf of X substituted in place of X ? How do I work out $f_Y(y)$ ? I hope it isn't asking too much for someone to tell me the general approach and also give an example of this, not necessarily one of those I mentioned above. I have looked in several statistics books but I didn't find the answer.


Suppose that A and B are mutually exclusive events for which P(A) = .3 and P(B) = .5. What is the probability that               (a)  either A or B occurs?      (b)  A occurs but B does not?      (c)  both A and B occur?       Trying to understand how to solve this problem, I solved (a) and got .15 (by .3*.5) for the probability for either of them to occur. However, the back of my text stated the answer for the problem was               3 · 4 · 4 · 3     -------------        (14 choose 4)                       = .1439      Could someone help explain to me what I'm doing wrong? This is my first stats class.


title: Using Pandas and statsmodels for ordinary least squares

text: Apologies in advance for the tedious beginner question. I'm trying to translate a least-squares problem from a manual process (using Excel for matrix transposition and multiplication) to using the Python statsmodels package. In this case, I'm performing an affine transform from a set of observed coordinates to a set of ground coordinates in eastings ( _E_ ) and northings ( _N_ ). I've used the following formula to form the **A** (design) matrix: $$ \begin{equation} f_i(a_{0}, a_{1}, a_{2}) = a_{0} + a_{1}x_{i} + a_{2}y_{i} \\ f_i(b_{0}, b_{1}, b_{2}) = b_{0} + b_{1}x_{i} + b_{2}y_{i} \end{equation} $$ Which gives me a matrix that looks like: $$ \begin{bmatrix} 1 & E_1 & N_1 & 0 & 0 & 0 \\ 1 & E_n & N_n & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & E_1 & N_1 \\ 0 & 0 & 0 & 1 & E_n & N_n \\ \end{bmatrix} $$ and a **b** vector, which is all _x_ co-ordinates, followed by all _y_ co- ordinates: $$ \begin{bmatrix} x_1 \\ x_n \\ y_1 \\ y_n \end{bmatrix} $$ I also construct a square covariance matrix, **W** , using the square of the standard errors of the eastings and northings, arranged in a diagonal. This is used as the weight matrix during the least-squares process (the standard errors are assumed to be independent) I then calculate: $A^TWA$, $A^TWb$, and $(A^TWA)^{-1}$, then multiply $(A^TWA)^{-1}$ by $A^TWb$ to determine a vector **x** , which contains values for $a_0, a_1, a_2$ and $b_0, b_1, b_2$ 5. Multiply **A** by **x** , and subtract **b** from the result, to determine a residuals vector, **v**. I can calculate the unit variance ($\sigma0$) by obtaining the square root of $\frac{v^TWv}{observations - unknowns}$, and multiplying it by the _a priori_ standard error of each co-ordinate, in order to assess the quality ( _a posteriori_ standard error) of the transform. I can also calculate the standard error of my **x** vector (the diagonal values of **Cx** ) by multiplying $(A^TWA)^{-1}$ by $\sigma0^2$ Now that the tedious step-by-step manual explanation is out of the way, let's say I have Pandas DataFrames for **A** , **b** and **W** :               In [124]: A_matrix     Out[124]:      <class 'pandas.core.frame.DataFrame'>     Int64Index: 108 entries, 0 to 107     Data columns (total 6 columns):     As             108  non-null values     Eastings_a     108  non-null values     Northings_a    108  non-null values     Bs             108  non-null values     Eastings_b     108  non-null values     Northings_b    108  non-null values     dtypes: float64(4), int64(2)          In [125]: b_vector     Out[125]:     <class 'pandas.core.frame.DataFrame'>     Int64Index: 108 entries, 0 to 107     Data columns (total 1 columns):     coordinates    108  non-null values     dtypes: float64(1)          In [162]: Weight_matrix     Out[162]:      <class 'pandas.core.frame.DataFrame'>     Int64Index: 108 entries, 0 to 107     Columns: 108 entries, 0 to 107     dtypes: float64(108)      How do I use `statsmodels.ols_regression` to easily calculate my residuals and $\sigma 0$?


title: prediction from incomplete observations

text: Suppose I have a linear model predicting class-membership from a set of predictors. Now, I am going to classify a new observation which has, however, some predictor values missing. How can I deal with such situation? I know there are methods for imputing the missing values but I would like to avoid this and to use only measurements that were really made.


title: Would a Random Forest with multiple outputs be possible/practical?

text: 1. Random Forests (RFs) is a competitive data modeling/mining method.   2. An RF model has one output -- the output/prediction variable.   3. The naive approach to modeling multiple outputs with RFs would be to construct an RF for each output variable. So we have N independent models, and where there is correlation between output variables we will have redundant/duplicate model structure. This could be very wasteful, indeed. Also as a general rule more model variables implies a more overfit model (less generalisation). Not sure if this applies here but it probably does. In principle we could have an RF with multiple outputs. The prediction variable is now a vector (n-tuple). The decision nodes in each decision tree are now splitting the set of target/prediction vectors based on a threshold vector, I figure this threshold is taken to be a plane in the n-dimensional space and that therefore we can determine which side of the threshold vector each of the target vectors is on. The optimal prediction value for each side of the decision split is the mean (centroid) calculated for the vectors on each side. Finding the optimal split point when working with single variables is trivial and computationally fast/efficient. For an n-tuple we cannot find the optimal split (or at least it becomes computationally infeasible as N increases), but we may be able to find a near optimal split using a Monte Carlo type method (or some hybrid of Monte Carlo and local gradient traversal). Would this actually work? That is, would it just map the training pairs without generalising? Does this technique already exist under a different name? You might also want to consider how this relates to neural nets such as Restricted Boltzmann Machines (RBMs) and Deep Belief Networks.


title: use of validation set on lasso cross validation

text: When training a model a train, a validation and test set are used. I was wondering if there is any paper or example that proves that the use of an independent validation set increase the performance of the lasso estimator. I am particularly interested in situations where the penalty value is chosen through cross validation


title: Regression model for predicting sales?

text: We sell machinery. The following is graph is an _approximation_ of units sold over time for a particular piece of equipment1 : ![http://i.stack.imgur.com/3m3rM.png](http://i.stack.imgur.com/uO6E3.png) It starts out slow and slowly grows over time. I tried using linear regression, but if all data points are considered, the RMSE is in the order of 120 days, which is completely unacceptable. If I only include the sales over the last 365 days (last 40 or so data points), the RMSE drops to a more acceptable 10 days. Is there a better way to do this without arbitrarily discarding data? It seems like there's a small seasonal factor, where sales appear slower in some months. How can I incorporate that into the regression? 1 The data points are ${(Date, ID)}$, where $ID$ is a number assigned incrementally for each sale, and $Date$ is the actual date the sale was completed. Because some sales complete faster, lower ID's can end up with slightly later dates.


title: When is the maximum value of chi square achieved for a non-symmetric table?

text: Let us say we have an RxC table, when R and C are not necessarily equal, when is the maximum value of the Chi Square statistic achieved? and how can this be proven?


title: Handling missing attributes when using weighted average

text: I was trying this optimization where I have some predictors and I want to predict the real value using weighted average of the predictors. What I am trying to do is given the predictors' prediction and the ground truth values, I will calculate the weighted average and the weighted average has weights defined k1,k2,...kN for the N predictors such that when I take the weighted average the variance of the error of the weighted average is minimized. Now at times, the predictors' predictions are missing. How can I handle these missing values?


title: How to visualize iterative parameter constraint?

text: I have conducted an analysis in which I start with a set of informed prior parameter distributions, and then conduct sequential analyses that constrain the distributions with data. I am currently using density plots to visualize these results, but I am wondering how I can most effectively convey the process of parameter constraint to a non-bayesian audience. **Questions** Is there a more effective way to visualize the sequential constraint of parameter uncertainty? Specifically:   * is there a way to I indicate the progression within the plot?    * would an alternative approach, such as a boxplot, be more effective in communicating the change in range and central tendency?   * would it help to offset the density=0 line so that the lines are not overlapping?   * would it be effective to use sparklines, e.g. with the 95%CI and medians indicated by a point and a value?   * is it appropriate to exclude the y-axis density scale since the area = 1? I have included two parameters here, although I will actually be plotting between six and fifteen. **Examples** Here are some example data:               n = 10000; set.seed(0)     prior       <- data.frame(theta1 = rnorm(n, 10, 3),                               theta2 = rnorm(n, 20, 1.5))     posterior1  <- data.frame(theta1 = rnorm(n, 11, 0.5),                               theta2 = rnorm(n, 22, 1))     posterior2  <- data.frame(theta1 = rnorm(n, 10.5, 0.3),                               theta2 = rnorm(n, 23, 0.8))      My current approach is along the lines of plotting the densities increasingly dark grey to black:               for(i in 1:2){      plot(density(posterior2[,i]), main = colnames(prior)[i],            xlim = c(0,20*i), xlab = '', ylab = '', col = 'red)      lines(density(posterior1[,i]), col = 'darkgrey')      lines(density(prior[,i]), col = 'grey')     }      These figures are from the above example code: ![enter image description here](http://i.stack.imgur.com/K7ZDB.png) Here is some of the actual data, with the x-axis scaled to focus on posteriors rather than priors, but at this scale the priors, despite being informed, look fairly flat. (These figures also illustrate that the posteriors are from mcmc chains rather than standard densities as in the example above). ![enter image description here](http://i.stack.imgur.com/dOTBs.png)


title: How to obtain confidence interval for lambda and the max using R Box Cox transformation?

text: I did a box cox plot on the ozone data in R. I need to determine the best transformation. Is there a way to get the exact confidence interval for lambda and the max. lambda or by just looking at the graph to estimate.(I dont know how to paste the graph).


title: Likert scale questionnaire and logistic regression

text: I have three independent variables , and one moderating variable. I'm studying the affect of these variables on the decision to adopt electronic commerce. I have a questionnaire, based on this questionnaire which is on a 7-point likert scale (from totally disagree which is coded 1 to totally agree which is coded 7). I have gathered data on on my three independent variables. Now I want to analyze the model. Question is should I calculate mean for every IVs and then do the logistic regression? Regards


title: Aggregating and comparing data across days using paired t-test

text: I have two groups of data (with and without stimulus) from the same subject. For each experimental session, I measure a score (let say reaction times, n_with_stim=50, n_without_stim=50). For each day, the subject participates in four sessions (hence, n_with_stim_total=200, n_without_stim_total=200). My ultimate goal is to test whether the stimulus has any effect. The problems are (i) his baseline score changes across session; The later in the day, the slower his reaction time, and (ii) his baseline score may change across day. How can I aggregate and compare the results across sessions? Can I just combine all the data across sessions within the day to see whether there is an effect on that particular day? Is it ok to combine data across sessions and days to compare? Is paired t-test an appropriate test for me to use? Thank you in advance.


title: What does this distribution look like?

text: I have data with the following empirical distribution![empirical distribution](http://i.stack.imgur.com/5okOh.png). What would be a good model for it? Edit: here is a closer view: ![enter image description here](http://i.stack.imgur.com/ZNnZY.png)


title: Compare regressions among more than 2 groups

text: I would like to know if there is a way to compare the relationship between a factor (length) and a response variable (lipid content) among multiple groups? I have given up on ANCOVA, as length is continuous. I am unsure if a multiple regression allows the use of a categorical variable (group membership) as an IV.


title: Quantile regression

text: I have a question regarding quantile regression. Supposing that I have 10000 observations with one response variable and several predictor variables in a dataset collected each year over several years. I run a multivariate quantile regression at 90% percentile. I want to compare the observed 90th percentile response value for each year (a single value) to the 90th percentile predicted response value for that year (a single value) based on the multivariate quantile regression model. How can I do that?


title: Which regression tree to use for large data?

text: I have a dataframe with 2 million rows and approximately 200 columns / features. Approximately 30-40% of the entries are blank. I am trying to find important features for a binary response variable. The predictors may be categorical or continuous. I started with applying logistic regression, but having so much missing entries I feel that this is not a good approach as glm discard all records which have any item blank. So I am now looking to apply tree based algorithms (`rpart` or `gbm`) which are capable to handle missing data in a better way. Since my data is too big for `rpart` or `gbm`, I decided to randomly fetch 10,000 records from original data, apply `rpart` on that, and keep building a pool of important variables. However, even this 10,000 records seem to be too much for the `rpart` algorithm. What can I do in this situation? Is there any switch that I can use to make it fast? Or it is impossible to apply `rpart` on my data. I am using the following rpart command:               varimp = rpart(fmla,  dat=tmpData, method = "class")$variable.importance


title: Probabilistic modelling MCMC question with pyMC

text: This is my first post and I am a newby in pymc. I am trying to model a non- linear system (see below for a further explanation). I create my synthetic data with:               data = np.random.normal(8., data_std, 1000)     obs = pm.Normal("obs", coupling, value=data, observed=True)      My stochastic variable is:               phi = pm.Uniform("phi", -180, 180., value=-150)      and it is related to the data though the deterministic function:               @pm.deterministic     def coupling(phi=phi, a=a,b=b,c=c):          """          Calculate the J-coupling  from Karplus forumula          """          ang=phi*3.1415/180.          return a*np.cos(ang)**2+b*np.cos(ang)+c      (a,b,c are scalars. In the future, I'll try to see what happens when they are also stochatic variables) Then my model is:               model = pm.Model([obs, phi])     mcmc = pm.MCMC(model)     mcmc.sample(20000, burn=1000, burn_till_tuned=False)      What I don't understand is why the posterior distribution of phi does not depend on the dispersion of the data. That is, I get the same std of `phi_trace` whether data_std is 0.1 or 2. What I am doing wrong?               phi_trace = mcmc.trace("phi")[:]     print "\n Phi std:", phi_trace.std()      Thanks in advance, Ramon PS. I know the nonlinearity of the system may cause the model to get trapped to one of the solutions, say close to -150 and not explore 0 or +150. I'll try to deal with that in the future... PS2. Some background: phi is a dihedral angle of a molecule and the coupling is what can be measured experimentally. In fact a uniform distribution for phi is a bad prior, but I wanted to see the difference when using more informative priors.


title: Self Study Qn on Exponential Distribution

text: I am currently working on a self-study question relating to exponential distribution. I made an attempt but was not able to make sense of the answer. When I looked up the reference guide, it offered the answer at 0.1068. Appreciate any guidance and advice on where has my method gone wrong. Many thanks. **_Question_** > When a truck arrives at a checkpoint, each truck will be checked by customs > agents. The times are exponentially distributed with a service rate of 10 > trucks per hour. What is the probability that a truck requires less than 15 > minutes but more than 10 minutes to be checked? **_Attempt_** ![enter image description here](http://i.stack.imgur.com/DaxrY.jpg)


title: String clustering and centroid computation

text: I have a text file document containing a set of words `strings` that I want to cluster. I want to use the `K-means` algorithm. As a distance I'll be using the `Levenshtein distance`. **Question** Usually if we have for example real number vectors we can just compute the mean of all the vectors and it will represent the centroid of the cluster. What approach can I use to compute the centroid of a cluster containing strings?


title: How to perform ANOVA on a modification of the dataset?

text: I just started working with R and I could really use some help with this problem. I want to perform an ANOVA, to see how my independent factors (between-subject: Group) (within-subject (Modality, Violation, Correctness) interact with my dependent variable (EndQuestion.ACC). EndQuestion.ACC means whether a participant gave an incorrect ("0") or a correct ("1") answer. ![Example dataset](http://i.stack.imgur.com/tqeas.png) I have the following function to perform the ANOVA.               fourf.aov <- ezANOVA(dat, dv=.(EndQuestion.ACC), wid=.(Subject), between=.(Group),                       within=.(Modality, Violation, Correctness), detailed=T)           However, when the times arrives to perform post hoc tests, of course, I get means between 0 and 1. But what I really want to use for the ANOVA and for the post hoc tests are the summed amount of errors per participant. Does anyone have an idea how I can solve this with code? Thank you so much in advance for helping me!


title: How to compute the weight matrix for WLS estimation of a multi-group ordinal CFA model

text: I am attempting to perform a two-group confirmatory factor analysis (CFA) of one continuous factor on six ordinal predictors in `OpenMx` (`R`) using (robust) Weighted Least Squares (WLS) estimation. While I am completely new to OpenMx (I only used `sem` and `lavaan` ) I think I have most things now: I have a way of estimating polychoric correlations/covariances and estimating the latent means in one group and a good idea on how to implement this model in `OpenMx`. I just have one problem, I need a weight matrix in WLS estimation. This should be an asymptotic covariance matrix of the elements of the observed correlations/covariances. What is this exactly and how can I compute such a matrix?


title: Growth curve analysis on orthogonal polynomial terms

text: I am conducting a study which is looking at the effect of 'Condition' (Quiet, Intelligible, Unintelligible) on the pupil(eye) response over time. Upon visual inspection of my data plots, pupil response slopes appear to change in a curvilinear way and there appear to be differences in the timing and amplitude of peaks and slopes between conditions at different timepoints. In order to describe and statistically analyse these differences I decided to use 'Growth Curve analysis' and analyse the effect of 'Condition' on the pupil response intercept, slope (linear), quadratic, cubic, quartic, and quintic terms to see if any of these terms in particular best describe any difference across conditions. So here is the base model which I constructed ('ot' refers to orthogonal time terms);               m.base <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)      To analyse, firstly, the effect of adding 'Condition' as a fixed effect on the intercept and, subsequently, the effect of 'Condition' on the slope and other time terms (quadratic, cubic, etc.) I entered the following code;               m.0 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)     m.1 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + ot1:Condition +   (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)     m.2 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + (ot1+ot2:Condition)   + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)     m.3 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + (ot1+ot2+ot3:Condition) + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)     m.4 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + (ot1+ot2+ot3+ot4:Condition) + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)     m.5 <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) + Condition + (ot1+ot2+ot3+ot4+ot5:Condition) + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)      So I guess I've got two questions relating to my analyses: (1) I'm not 100% sure that I have coded the fixed and random effects correctly, so any advice would be much appreciated. I believe that there will be variability across participants and have therefore including participants ('ParticipantNumber') as a random effect. I can't think of any reason to add a 'ParticipantNumber|Condition' random effect also although I have noticed that this is done in lots of other cases, so maybe there is a good reason for doing this (which I'm unaware of)? (2) This is my output when I run an ANOVA to assess model fit;               Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)         m.base 28 -1630.8 -1476.9 843.38  -1686.8                                  m.0    30 -1626.8 -1461.9 843.38  -1686.8  0.000      2          1         m.1    32 -1651.7 -1475.9 857.87  -1715.7 28.977      2  5.102e-07 ***     m.2    32 -1639.8 -1463.9 851.90  -1703.8  0.000      0          1         m.3    32 -1641.1 -1465.2 852.53  -1705.1  1.264      0  < 2.2e-16 ***     m.4    32 -1631.8 -1455.9 847.88  -1695.8  0.000      0          1         m.5    32 -1646.0 -1470.1 855.00  -1710.0 14.234      0  < 2.2e-16 ***      My next step is to run pairwise comparisons between the three levels of my 'Condition' variable to uncover where the difference lies for linear, cubic, and quintic terms. However, I haven't yet been able to find a straightforward way of doing this on R. Any suggestions from people who have run similar analyses would be much appreciated.


title: Why is a regression model of portfolio return giving smaller adjusted R-square (i.e., negative) than expected?

text: I have a question for adjusted $R^2$ given a specific regression model. I am doing a project on January effect and I have a model from some journal using $$R_i = a_0 + a_1D_{\mathrm{Jan}} + \varepsilon_i$$ where   * $R_i$ is daily return of portfolio/index,   * $a_0$ is non-January daily returns,    * $a_1$ is January returns over non January returns,    * $D_{\mathrm{Jan}}$ is dummy variable (1 for Jan, 0 otherwise).  With this model I tried to test whether January return is significantly greater than non-Jan returns, especially in small capitalization stocks. So I have returns for size sorted portfolio 1 to 4, where Portfolio 1 (P1) consist of smallest cap stocks, P4 consist of largest cap stocks. What I did was using regression program in Excel and SPSS, input all daily returns of the portfolio from Jan to Dec for 10 years as dependent variables and the dummy (1 for january and 0 for others) as an independent variable into the regression program and I get **negative adjusted $R^2$** for all portfolios (P1-P4), mostly about -0.05. The value of $R^2$ itself was also very small at about 0.0001. The journal I am basing my model on was using monthly returns as $R_i$, but I modified it into daily return. I still get negative adjusted $R^2$ when I use the monthly returns. Can anyone please help me point out what was wrong in the model? Did I use the wrong input based on the given model? If yes, what should be the correct input? Or if the model was wrong, what model should I use? Here are the results of my test, the p-value of the dummy variable and the adj $R^2$.   * For P1 0.54, adj $R^2$ -0.0002.    * P2 0.36, adj $R^2$ 0.0004.    * P3 0.68, Adj $R^2$ -0.0003,    * P4 0.14, adj $R^2$ 0.0005. I understand that the variables are all insignificant. But I am confused in interpreting the $R^2$


title: Visual representations of the p-value in ANOVA to assist intuitive understanding

text: I'm a statistics newbie and I would like to make sense of the p-value with ANOVA, hopefully with a visual presentation. All online visual tools I've found so far only have MSB, MSE and F values visually presented (such as this here. I'm a very visual person and I would like to see how it works before I learn more about it. Do you have any suggestions for me? It could be a Matlab program as well, something where I play with the distributions and see the p value change.


title: Experiments for performance testing software

text: I am looking to run a series of performance tests for a software system across different configurations and determine which combination works best. The software system processes a number of files, does some work with them and writes them to a database. In order to carry out these tests, I have a script that will run the software using **varying file sizes** (10, 100, 500, 1000), **varying buffer sizes** (0, 10, 100, 300, 500), and **two different "run-configurations"**. Each combination is run 30 times and the execution time is recorded for each combination for each run-configuration I am looking to test. The same files are used for each run and the database is cleared to ensure independent trials. My first question is whether or not this is a statistically sound way of conducting these experiments. Would the assumption that observations are independent hold? I would like to run factorial ANOVA across the file sizes, buffer sizes, and configurations to determine if there are interaction effects between them and the execution time. My next question is what should I do if the execution time variability across groups is roughly equal but execution times within the various groups do not follow a normal distribution?


title: Testing a 2 by 2 mixed ANOVA in R

text: I have data from an experiment that is testing how the order of two studying methods (visual or auditory) affects word recall. For analysis a multi-factor anova with a repeated measure is appropriate, but I am not sure if I am structuring my data correctly. This is the command I'm using: `aov(recalled_items~task*order)+Error(subject/task)+(order))` Here is an example of the data structure:               Subject    Task    Order    Recalled Items     A        Visual    First       13     A       Auditory   Second      22     B        Visual    First       14     B       Auditory   Second      28     C        Visual    Second      10     C       Auditory   First       15     D        Visual    Second      14     D       Auditory   First       29        * Does R know to compare Visual 1 and Visual 2 recall values and Auditory 1 and Auditory 2 recall values?  I am worried that because of the way I structured my data R is just comparing Visual 1 and Auditory 1 and as a result I am getting no effect.


title: Recommended procedure for factor analysis on dichotomous data with R

text: I have to run a factor analysis on a dataset made up of dichotomous variables (0=yes, 1= no) and I don´t know if I'm on the right track. Using `tetrachoric()` I create a correlation matrix, on which I run `fa(data,factors=1)`. The result is quite near to the results i receive when using MixFactor, but it´s not the same.   1. Is this ok or would you recommend another procedure?   2. Why does `fa()` work and `factanal()` produce an error? (`Fehler in solve.default(cv) : System ist für den Rechner singulär: reziproke Konditionszahl = 4.22612e-18`)


title: Justification for the Bootstrap Percentile Interval

text: The following is a proof of the validity of the (bootstrap) percentile interval taken from Larry Wasserman's "All of Statistics". $\theta^*_{\alpha/2}$ and $\theta^*_{1-\alpha/2}$ denote the respective quantiles of the bootstrap of our estimator $\hat {\theta}$. Could you please take a look at the proof below and explain to me where the third line of the derivation comes from? How does $U$ enter the derivation? Is there perhaps a mistake there? Thank you in advance. ![enter image description here](http://i.stack.imgur.com/e4fpL.jpg)


title: Expected number of unseen cards when drawing $2n$ cards from a deck of size $n$

text: We have a deck of $n$ cards. We draw cards from it uniformly at random with replacement. After $2n$ draws what is the expected number of cards never chosen? This question is part 2 of problem 2.12 in > M. Mitzenmacher and E. Upfal, _Probability and Computing: Randomized > Algorithms and Probabilistic Analysis_ , Cambridge University Press, 2005. Also, for what it's worth, this is not a homework problem. It's self-study and I'm just stuck. My answer thus far is: Let $X_i$ be the number of distinct cards seen after the $i$th draw. Then: $E[X_i] = \displaystyle \sum_{k=1}^{n} k (\frac{k}{n}P(X_{i-1}=k) + \frac{n-k-1}{n} P(X_{i-1}=k-1))$ The idea here is that each time we draw, we either draw a card we've seen or we draw a card we have not seen, and that we can define this recursively. Finally, the answer to the question, how many have we **not** seen after $2n$ draws, will be $n-E[X_{2n}]$. I believe this is correct, but that there must be a simpler solution. Any help would be greatly appreciated.


Can chi-square test be performed on data that does not have normal or gaussian distribution? One of the chi-square tests I am performing results in percentage points less than the 12.59 (α=0.05, v=6) using cross-tabulation. Therefore, I cannot reject the null hypothesis. Is there a test that I can perform to confirm the alternative hypothesis? Added clarification 30nov2011: I will attempt to clarify some of the issues raised in these comments. Please be gentle on me as I am not trained in statistics. In fact, thanks to the unimaginative curriculum and teaching, I abhorred statistics taught in engineering school. Just trying to help my wife with a project. So, we are trying to establish a weak (if at all) relationship/dependence between 2 different attributes of a dataset. To accomplish this, we have generated a cross-tabulation. One of the attributes matches perfectly with a lognormal distribution (acc. to Minitab) as the data points and the lines are almost coincident. We have used the k-means clustering algorithm to group the lognormal attribute into 3 broad categories (large, medium and small). Next, we are trying to perform a chi-square test on the cross-tabulation. df = 6. But, not having sufficient insight into this, I am not sure what the result of 6.8 percentage points supposed to mean for the null hypothesis. Per our understanding, the null hypothesis is "two attributes are not related". I am wondering what the next step should be. Further progress hinges upon confirming or denying the existence of a relationship. (Note, we do not have to determine the relationship. However, I wouldn't mind getting some insight into the relationship in this process.)


title: Removing human evaluator bias

text: I am working on machine translation evaluation, and am looking at ratings given by humans who are judging the quality of sentences produced by machine translation systems. The evaluators give fluency and adequacy (how well the information is preserved) scores to sentences. The data is from NTCIR7 for those that are interested. For the test data I have, there were 3 human evaluators, and they rated 100 sentences, from 5 different systems, so 500 sentences in total. The scores they gave were between 1 and 5 (which I normalised to 0.0-1.0). **Problem:** It seems that one human evaluator **A** is biased towards giving much more generous _fluency_ scores. You can see this in the rough chart below, that **A** gives many more sentences a perfect score of **1** when compared to **B** and **C**. **Question:**   * What statistical methods can I use to try to remove this bias? Is bias even the right word?   * With only 3 evaluators, is it a good idea to modify the scores at all? Does that count as 'cheating' to get the answer I think is right? **Graphs:** (sorry, the highest score is on the left) **Fluency: A** gives many more perfect **1** scores than **B** or **C** ![Fluency correlation](http://i.stack.imgur.com/j1tDJ.png) **Adequacy: C** gives more low **0.25** scores, but fewer **0** scores. I am less worried about this as it seems to average out. But it is still kind of troubling. ![Adequacy correlation](http://i.stack.imgur.com/orTft.png)


title: How do I interpret conflicting results from adf.test and pp.test in R?

text: I have difficulty in solving unit root test by using both `adf.test` and `pp.test` in R. The result p-value using `adf.test` is 0.2677 (mean that it has unit root), and the result p-value using `pp.test` is 0.01696 (mean that is has no unit root). What means that? I wonder whether I have to use diff function to remove unit root or not. (`k` argument in `adf.test`, How do I use this? what means exactly for `k`?) Raw data is below.               2935833     2622529     2719635     2625179     2311187     2101758     2552638     2883423     3128904     2959348     2759000     2233755     2560858     2548821     2625675     2326076     1662956     1772409     1797275     2639852     2799990     3133285     2438296     2583766     2610157     2493415     2094163     2174301     2283420     2505128     2873785     2339727     2985829     3037351     1828265     1038562     1474727     1523331     2122667     2571006     2252161     2422347     2155973     2294976     2809652     2436293     2561852     2199544     2674423     2551363     3110508     3177925     3046952     2850904     3002830     2910913     2809172     3136842     3355368     3604565     3013310     3125751     2548605     2646575     2231458     1962095     1958019     2143073     2305966     2620302     2356447     2427571


title: Distribution of longest streak of consecutive (non-)birthdays in birthday problem

text: Given the setting of the birthday paradox ($n$ people with birthdays randomly drawn from the uniform distribution of 365 days), let us call the days on which at least one person was born _birthdays_ and all other days _non- birthdays_. What can be said about the distribution of a) The length in days of the longest streak of birthdays; b) The length in days of the longest streak of non-birthdays? Is there an analytical solution to, for example, the expected value of either of these? (You can choose whether to break all streaks when the year changes; I think insightful analysis is more interesting than the precise formulation of the problem.)


title: Adaboosted trees and q-learning

text: Is there a solid mathematical connection between ADABOOST and Q-learning that can inform "higher performance" machine learning algorithms? Q-learning references:   * http://www.autonlab.org/tutorials/rl.html   * http://mnemstudio.org/path-finding-q-learning-tutorial.htm   * http://www.autonlab.org/tutorials/rl06.pdf (pg 20) which comes from the (beautiful) autonlab tutorials. *todo:    * _populate with ADAboost references_   * _justify the intuition behind the question_


title: Maximum Likelihood Estimation with Known Parameter Distribution

text: Consider i.i.d observation vector ${\bf x}$ from a distribution $F$ depending on vector of parameters $\boldsymbol{\theta}$ and single parameter $\alpha$. We would like to estimate parameters [$\boldsymbol{\theta}, \alpha$] using maximum likelihood estimation but we also know that $\alpha$ is a realisation from a $N(\mu, \sigma^2)$. If we didn't know the distribution of $\alpha$ we could just maximise the objective function $f({\bf x} | \boldsymbol{\theta}, \alpha)$ but knowing the distribution for $\alpha$, what does the objective function to maximise become?


title: Utility function for survey questions

text: What exactly is a utility function of a survey question? For example, suppose a question has $7$ possible answer options. We assign each answer choice a number between $0$ and $100$. What information does this tell us?


title: AR(1) on autocorrelated data that is not a time-series

text: I need to apply a regression model on observations that is not time series data but each observation presents a store and the amount of cartons that gets sent to that store. For instance `a_i = number of cartons sent to store i` and is my dependent variable. Each carton contains products of different kinds and my independent variables is the number of units for a certain product for each observation. Thus `x_ij = number of units of product j sent to store i`. So essentially i need to forecast the number of cartons given the amount of units of each product that needs to go to a store. The problem i have is that the data is autocorrelated, where smaller carton forecast tend to have bigger errors and larger carton forecast smaller errors. I would like to use AR(1) to transform the data with the following formula, `Y_t - pY_(t-1) = B_0(1-p) + B_1(X_t - pX_(t-1))` My model will have more than one independent variable, but it is correct of me to first sort the data ascending according to the dependent variable `a_i` before applying the above formula as `i` does not represent time but rather a store?


title: How to specify a contrast matrix (in R) for the difference between one level and an average of the others?

text: I have a regression model that looks like this: $$Y = \beta_0+\beta_1X_1 + \beta_2X_2 + \beta_3X_3 +\beta_{12}X_1X_2+\beta_{13}X_1X_3+\beta_{123}X_1X_2X_3$$ ...or in R notation: `y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x1:x2:x3` Let's say $X_1$ and $X_2$ are categorical variables and $X_3$ is numeric. The complication is that $X_1$ has three levels $X_{1a}, X_{1b}, X_{1c}$ and instead of standard contrasts, I need to test:   * Whether the intercept for level $X_{1a}$ significantly differs from the average intercept for levels $X_{1b}$ and $X_{1c}$.   * Whether the response of $X_2$ is significantly different between level $X_{1a}$ and the average of levels $X_{1b}$ and $X_{1c}$.   * Whether the slope of $X_3$ is significantly different between level $X_{1a}$ and the average of levels $X_{1b}$ and $X_{1c}$. Based on this post it seems like the matrix I want is...                2     -1     -1      So I do `contrasts(mydata$x1)<-t(ginv(cbind(2,-1,-1)))`. The estimate of $\beta_1$ changes, but so do the others. I can reproduce the new estimate of $beta_1$ by subtracting the predicted values of the $X_1b$ and $X_1c$ group means (when $X_3=0$ and $X_2$ is at its reference level) from twice the value of $X_1a$ at those levels. But I can't trust that I specified my contrast matrix correctly unless I can also similarly derive the other coefficients. Does anybody have any advice for how to wrap my head around the relationship between cell means and contrasts? Thanks. Is there a standard name for this type of contrast? * * * Aha! According to the link posted in Glen_b's answer, the bottom line is, you can convert ANY comparison of group means you want into an R-style contrast attribute as follows:   1. Make a square matrix. The rows represent the levels of your factor and the columns represent contrasts. Except the first one, which tells the model what the intercept should represent.   2. If you want your intercept to be the grand mean, fill the first column with all of the same non-zero value, doesn't matter what. If you want the intercept to be one of the level means, put a number in that row and fill the rest with zeros. If you want the intercept to be a mean of several levels, put numbers in those rows and zeros in the rest. If you want it to be a weighted mean, use different numbers, otherwise use the same number. _You can even put in negative values in the intercept column and that probably means something too, but it completely changes the other contrasts, so I have no idea what that's for_   3. Fill in the rest of the columns with positive and negative values indicating what levels you want compared to what others. I forget why summing to zero is important, but adjust the values so that the columns do sum to zero.   4. Transpose the matrix using the `t()` function.   5. Use `ginv()` from the `MASS` package or `solve()` to get the inverse of the transposed matrix.   6. Drop the first column, e.g. `mycontrast<-mycontrast[,-1]`. You now have a p x p-1 matrix, but the information you put in for your intercept was encoded in the matrix as a whole during step 5.   7. If you want labels in the summary output more pleasant to read than `lm()` et al.'s default output, name your matrix's columns accordingly. The intercept will always automatically get named `(Intercept)` however.   8. Make your matrix the new contrast for the factor in question, e.g. `contrasts(mydata$myfactor)<-mymatrix`   9. Run `lm()` (and probably many other functions that use formulas) as normal in standard R without having to load `glht`, `doBy`, or `contrasts`. Glen_b, thank you, and thank you UCLA Statistical Consulting Group. My applied stats prof spent several days handwaving on this topic, and I was still left with no idea how to actually write my own contrast matrix. And now, an hour of reading and playing with R, and I finally think I get it. Guess I should have applied to UCLA instead. Or University of StackExchange.


title: How to model a skewed Student's t disribution

text: I have a small number of samples (5) of a large population (~10,000). The samples are percentages and hence I know from the context that no answers are possible below 0% or above 100%. From this one set of samples I've calculated my t statistic and confidence intervals; however, for anything above 95% confidence my interval width is such that the predicted upper limit is above 100. My understanding is that is you have specific knowledge like this you can create a skewed Student's t distribution. First, is this correct? Second, how would this apply to the upper and lower limits of my confidence interval? I understand that two values are important in calculating this, the skewness and kurtosis, but I'm not sure how to use them. My skewness value: -1.82 My kurtosis value: 3.448 Sample Values:               79.465     91.905     91.096     88.144     90.599


title: R linear regression test hypothesis for zero slope

text: I'm working with some data and I used R to the a linear regression model `Y = aX + b`. The code I used was `summary(lm(Y~X))` What I got was                             Estimate  Std. Error t value Pr(>|t|)         (Intercept) -0.3884045  0.0260232  -14.93   <2e-16 ***     X            0.0062095  0.0004635   13.40   <2e-16 ***      What I want to test now is the null hypothesis `H0: a=0`, that is, the case where the slope is zero. I'm confused about how to do that. I tried using the offset parameter ( the idea was to subtract the 'a' coefficient found previously in the former fit ), but I'm not sure it is the correct way to test this hypothesis. What I did was: `summary(lm(Y~X,offset=0.0062095*X)` and I got:                            Estimate   Std. Error t value Pr(>|t|)         (Intercept) -3.884e-01  2.602e-02  -14.93   <2e-16 ***     X           -2.464e-08  4.635e-04    0.00        1        Is it right? Am I now supposed to reject H0 since the p-value found was 1?


title: Interpreting p-values

text: I am a little confused about what p-values mean under Fisher's signficance testing & Neyman-Pearson's hypothesis testing. Fisher uses p-values as a continuous measure of evidence against a null hypothesis? So a p-value of 0.06 would indicate that there is no difference and the null hypothesis is true? However, does this mean the same thing under Neyman-Pearson. I know that you have to pre-set alpha values for type 1 errors but does this affect p? Does a p-value greater than alpha indicate that there is >5% chance of a type one error occurring?


title: Check MCMC stationarity and convergence before or after thinning?

text: library(BRugs)     library(coda)     ...     output = BRugsFit("model.bug", data, inits, numChains = 1, nThin = 5,               seed=6, parametersToSave, nBurnin = 9, nIter = 50,coda = T,digits = 5)          # computes the effective sample sizes.     effectiveSize(output)          HPDinterval(output,prob=.9)          # diagnostics      geweke.diag(output)     raftery.diag(output)     gelman.diag(output)      `output` will contains the iterations after removing burn-in period (here it is 9) and thinning (here 5). So there are 50/5=10 iterations in `output`. I wonder whether to check the diagnostics (geweke, raftery, gelman, and possibly more) before or after removing burn-in and/or thinning? Thanks!


title: Response variables with interaction effects

text: Does anyone know of a statistical approach/method that lets you model linear regressions between 3-5 explanatory variables that are categorical or continuous and 3 continuous response variables that have interaction effects? Thanks in advance.


title: Why are different answers obtained for a regression F-test when using SSE and R2 formulas?

text: I cannot understand what I'm doing wrong with an F-test for the regression - I get completely different answers using SSE and R2 formulas. Unrestricted model: ![enter image description here](http://i.stack.imgur.com/aVWhX.png) Restricted model: ![enter image description here](http://i.stack.imgur.com/4dxNO.png) 1 restriction is imposed, so q=1 in both cases. Calculation with SSE: $F=\frac{SSE_r- SSE_u}{SSE_U}*\frac{N-K}1=\frac{2.252-2.142}{2.142}*48=2.465$, same result given by Stata F-test. Calculation with R2: $F=\frac{Ru^2-Rr^2}{1-Ru^2}*48$=1216.94   What am I doing wrong?!


title: With a binary Y, why are R's lowess fits so often flat?

text: I'm trying to explore (rather than formally test) relationships for various subgroups between a continuous X and a binary Y (0 or 1). A good shortcut would seem to be to graph a smoothed X-Y fit for each subgroup. But even when, at various levels of X, the mean of Y ranges from .1 to .3, R's smoothed fitline often simply hugs the horizontal where Y=0. This is not informative at all and in fact seems inaccurate. I've tried several different ways to obtain a smoother as well as many different values of f. E.g.:               x=runif(500,0,1)  #x is continuous     y=1:500     y=ifelse(y < 100,1,0)  #y is now binary with mean ~ 0.2          z=1:500     z=ifelse(z<50 | (z<300 & z>100),1,0)  #z is binary too               plot(x [z==1], y [z==1])     m=lowess(x[z==1], y[z==1])     lines(m, col='blue')  #The line hugs the points at y=0     windows()     scatter.smooth(x [z==1], y[z==1])         In this case this last line did produce the expected sort of fitline, but it's not a reliable occurrence when I work with real data.


title: Optimal rate of convergence for nonparametric estimators in Sobolev space

text: Given a regression model on interval $[0,1]$ $$ Y_{i}=f(x_{i})+\epsilon_{i},\ i=1,\ldots,N $$ with fixed design and standard error assumptions $E(\epsilon_{i})=0;\ E(\epsilon_{i}\epsilon_{j})=\delta_{i,j}\sigma^{2}$. The regression function $f$ is from Sobolev space $$ f\in W^{q}\left[0,1\right]=\left[f,\ldots,f^{(q-1)} \text{ are absolutely continuous, } \int_{0}^{1}\left| f^{(q)}(x)\right| ^{2}<\infty\right] $$ Is the optimal convergence rate with respect to $L_{2}$-norm $N^{-q/(2q+1)}$ ? If yes, could you give a reference? If not, what additional assumptions are required ?


title: Why do we use only the first component scores to construct different indices?

text: I have seen that in constructing indices by PCA the scores of only the first component is used. Although the first principle component explains the most variation contained in the data, often the VAF is too small (10% or even less). Then why do we use only the first component score instead of combining the scores of first few dimensions which would explain more variation when combined (I don't know if combining these score on different dimensions is legit)?


title: Confusion related to Singular Value Decomposition (SVD)

text: I have this problem where I am given a matrix A of size nxp. I run svd on this matrix A. Then I am supposed to plot this matrix using imagesc in matlab where rows are sorted by the first row eigenvalue and columns are sorted by the first column eigenvalue I didn't get what it means by "rows are sorted by the first row eigenvalue and columns are sorted by the first column eigenvalue" Suggestions?


title: Using a gamm4 model to predict estimates in new data

text: I have been experimenting with gamm4 to derive GAMMs of some repeated measures data. The models looks very nice and seem to give more flexibility than my LMMs. Ultimately I want to compare models not by the quality of their fit (also the reality of comparing LMM and GAMM fits seems complex?), but in the quality of their predictions in new data sets, and in simulated new data by MCMC. With LMMs I predict using the fixed effects only using:               mm <- model.matrix(terms(lmer),newdata)          newdata$predicted <- mm %*% fixef(lmer)      This is fine since we are predicting in new individuals, with new independent random effects. I cannot get this predict method to work with gamm4.               > mm <- model.matrix(terms(gamm4$mer), newdata)          Error in model.frame.default(object, data, xlev = xlev) :        variable lengths differ (found for 'X')      I think this is because the GAM process creates new variables in order to transform predictor variables. It is also complex, because I believe the transforms are stored as random effects, so I would need to extract these random effects, but not the “individual level” random effects. Does anyone know how I can:   1. Extract just the transform effect terms from the gamm4 model?   2. Make predictions into new data using gamm4?   3. Extract the model specification of the GAMM so I could implement it as a standalone algorithm?   4. General advice?


title: How to calculate 2D standard deviation, with 0 mean, bounded by limits

text: My problem is as follows: I drop 40 balls at once from a certain point, a few meters over the floor. The balls roll, and comes to a rest. Using computer vision, I calculate the center of mass in the X-Y plane. I am only interested in the distance from the center of mass to each ball, which is calculated using simple geometry. Now, I want to know the one-sided standard deviation from the center. So, I would be able to know that a certain number of balls are within one std radius, more balls within 2*std radius and so on. How do I calculate the one-sided standard deviation? A normal approach would state that half of the balls be on the "negative side" of 0 mean. This of course makes no sense in this experiment. Do I have to make sure that the balls conform to the standard distribution? Thank you for any help.


title: Linear regression sample size advice

text: I'm involved in sample size calculation for an oncology clinical trial. Our Primary outcome is quality of life (quantitative, normal). There are two treatment regimes. Therefore I started with a simple t-test sample size calculation. With an effect size of 0.55, pow=.8, alpha=.05 we need overall 106 participants. Considering lost to fup and death before QOL measurement accrual size rise to 152 participants. Trial is multicentre, with randomization stratified by centre. Patients can have different type of tumor (nsclc, cholangiocarcinoma, stomach, pancreatic) The funding institution referees requested to (freely translating) “increment sample size to account heterogeneity of patients involved". They did not request stratified randomisation center*tumor (maybe due to likely imbalance for a small sample size), eventually suggested a minimization. I'm starting to think about a linear model like the following, to fullfill requests QOL = f(Treatment.dummy, TumorType.dummies) Typically effect size relies on the R^2 with and without the covariate of interest. However given the state of the literature, i don't know how to make an educated guess about this stuff. The only idea that remains (at least to me) is to simulate with Monte Carlo. In the end it would be an F test of model with vs without Treatment.dummy. Given a certain n (growing beetween simulation steps):   * I would simulate proportion of different tumor patients recruited, given infos from our clinical database.   * Given the tumor type I would simulate proportion of patients randomized to Treatment in that strata (...thinking about Uniform with mean .5 ... eg U(0.35,0.65)   * For pts in control group, given tumor type, i would simulate QOL given normative data of the instrument    * For pts in treatment group, given tumor type, i would simulate QOL given normative data of the instrument AND effect-size of interest Then do the regression model, the test and get the power. Another way would be to start from the t-test sample, and then increment it given models sample size rule of thumbs (eg 20 pts per added variable); but I don't like it very much because you loose grip with power-analysis. Any other approach? Any suggestion (even “yes, do the Monte Carlo simulation”) would be very appreciated.


title: Time series with external regressors

text: I have monthly data for about the last 2.5 years for an economic index as well as external regressors for each month. Is it feasible to treat the months as independent and just fit a multiple regression model regressing the external data (3-5 variables) on the index? I have tried fitting an ARIMA (1,0,0) and with external regressors and the auto regressive term is not significant. Is there a test that will determine whether a time series is independent given external data? Also, given my small amount of data what is the best way to validate any models.