title: Understanding condition index used for finding multicollinearity

text: In linear models, in my book, the condition index is defined as $\sqrt{\lambda_{max} \over \lambda_{min}}$ where $\lambda_{max}$ is the maximum eigenvalue of $ZZ^*$, i.e., the correlation matrix of the independent variables, and $\lambda_{min}$ is the minimum eigenvalue. I could not find out why the condition index helps me to find out multicollinearity. So, I've come up with my own explanation: > In principal components analysis, eigenvalues explain the variance of data > in each of the eigenvectors. Thus, having a large ratio of max eigenvalue to > min eigenvalue means that the data can be explained by a smaller number of > eigenvectors than the number of independent variables, which means that > there is collinearity among independent variables. Given these, here are my questions.   1. Is my own explanation correct?   2. What's a better explanation? (Please understand that I do not have strong linear algebra knowledge. So, please explain in very easy words.)


title: In general how do you set K in K-NN?

text: As the title suggests, how should you set K in K-Nearest Neighbours? Is it just a case of lower values of K are more susceptible to over-fitting and larger values of K are likely to give a more accurate reflection (less susceptible to noise). Also the optimal value of K largely depends on the training set, but I was wondering whether there was a general 'technique' that is used?


title: How do you do time series cross-validation using python?

text: Also, any tutorials/blogs available that you are aware of?


title: Nested logistic regression in R

text: I am trying to find a way to do Nested Logistic Regression in R that fits my needs. I have a very large data set with almost 200 variables available. I have found my "best" model and it contains 12 of these variables. I am aware of the problems of using 12 variables in logistic regression and don't want to cover those here. I am using the model to predict the probability a student will drop out of high school. The dataset I'm using that includes the students I am predicting the probability for has some missing data. So for example say my model is:               Dropout ~ AttendanceRate + Enrollments + AgeDiff + Race + Absences      Now say it is not known for a particular student how many Absences he has had. I am using the predict function to find the dropout probabilities from my model. If one of the variables is missing or "NA" then there is no predicted probability as it is NA. To fix this problem I want to use Nested logistic regression. So for a student with a missing value in "Absences" the model would be               Dropout ~ AttendanceRate + Enrollments + AgeDiff + Race      I am currently using the glm function to do my logistic regression. I understand it cannot do nested logistic regression. So how do I do Nested Logistic regression from my full model so that I can predict a students dropout probablity no matter which variable or variables may be missing. I realize using Nested models will hurt the strength of my model for prediction purposes but I am willing to sacrifice that for my purposes. Any Suggestions?


title: Which test should I report, ANOVA or Mann-Whitney?

text: My data set is small, with roughly 12 in each group (two groups). In addition, the distribution of the data is not normal. As a result, I ran a Mann-Whitney U test. However, I also ran a one-way ANOVA and got the exact same results. Which test should I report and it is okay to say that I ran both with equal results but chose to report only the ANOVA?


title: Subdivisions in statistics

text: I was having a discussion with a statistics professor a while ago about the different 'flavours' of statistics (frequentist, Bayesian, ...). He posed that he would subdivide statistics in four categories: non-parametric-, robust-, frequentist- and Bayesian statistics. The subdivision is characterized by the amount of assumption the methods make about underlying distributions (non- parametric statistics makes none, while Bayesian makes those assumptions very explicit). I was going to to ask if CrossValidated agrees with this subdivision, but since that is a subjective question I'll ask: 1) Is this subdivision widely recognized in statistics; 2) do 'real world' problem usually require one particular method? Ie, given some problem, is there a method most suitable for solving that problem or can multiple methods work for a given problem? Thanks in advance.


title: Are truncated numbers from a random number generator still 'random'?

text: Here 'truncating' implies reducing precision of the random numbers and not truncating the series of random numbers. For example, if I have $n$ truly random numbers (drawn from any distribution, e.g., normal, uniform, etc.) with arbitrary precision and I truncate all the numbers so that finally I end up with a set of $n$ numbers, each with exactly 2 digits after the decimal. Can I call this new set of numbers 'random'? I came up with this question when I was reading about **hardware generated random numbers**. The Wikipedia article says _they generate random numbers by measuring a physical process._ But since this measurement has its limitations (measurement error, finite precision, etc) can we call these hardware generated numbers random?


title: How to interpret relatively good -2LL value & Nagelkerke $R ^2$ without any significant IVs

text: My model's chi square says my 'full' model isn't significantly better than the null model. But my -2LL is 80 (isn't this pretty _low_?) and my Nagelkerke $R^2$ is .20. On the other hand, there isn't a single significant IV in my model... So how should I interpret the lower -2LL and the not that small .20 Nagelkerke $R^2$?


title: General Non-linear Regularized Models

text: Had a general questions. Are there any good non-linear models with regularization? I've heard of some linear models with regularization but not too many non-linear ones. I understand that you can use SVMs, which have built in regularization, but just curious what else might be available. Thanks


title: One-way repeated measures ANOVA and missing data

text: I have following data for the one-way repeated measures ANOVA test and I don't know how to deal with missing values. The data looks something like this:               users  | type 1 | type 2 | type 3 | type 4 | type 5          user 1 |    x   |    x   |    -   |    -   |   x     user 2 |    x   |    x   |    x   |    -   |   -     user 3 |    x   |    x   |    x   |    x   |   -     user 4 |    x   |    x   |    -   |    x   |   x     ...       Can anyone help me out?


title: Reporting the use of different likert-scales in survey for structural equation modelling?

text: Suppose you are doing a survey and from the data you want to build a structural equation model. The research model consists of 5 constructs and each of them has between 2 and 5 items. The survey uses Likert scales but the scales differ between the constructs: 2 constructs use 7-point likert scales, one uses a 6-point likert scale and two use 4-point likert scales.   1. Is it a problem from a methodological point of view to mix these different scales when doing structural equation modelling? (perhaps especially when odd and even number scales are mixed?)   2. Should I report the use of different likert-scales somewhere in the body of the paper, e.g. in the limitations section?


title: Correlation coefficient for use with nonlinear finite sets

text: I need to find the correlation between two data sets that are nonlinear and only have values between 0 and 10. Both sets contained the same number of values. I've tried using Pearson's R and introducing some linearity by adding n+10 to each data point. That however introduces strong correlation where there really shouldn't be any. I'm completely new to statistics, so here is a layman's explanation of what I'm trying to achieve. I have a questionnaire that asks customers for an overall score between 0-10. I then ask several more specific questions for things like Price, Service, etc. These questions are also 0-10. I want to find the degree of correlation between the overall score and the scores for the specific questions. For example, if customers that give a high score for Price also give a high overall score, then I can see that Price is important to the customer. Whereas if their scores for Service do not closely match the overall score, then the quality of Service is not as important to the customer. A perfect correlation in my data would look like this: Price: 10, 10, 10, 10 Overall score: 10, 10, 10, 10 Could someone suggest a suitable coefficient to achieve this? Cheers


title: How to analyze this dataset with 2 between-, two within- factors using lme4?

text:
This post is a follow up on my previous post (which was interested in `lme`) and uses the same dataset. Now I would like to know how to analyze it using `lme4`. ## The data The data is from a behavioral experiment in which participants in 6 groups (based on two crossed factors) worked on 16 trials (two crossed 4-level factors). That is, we have a dataset `d` with two between-subject factors, `group` and `condition`, and two within-subject factors (i.e., repeated- measures factors), `topic` and `problem` (I uploaded the data to pastebin, so everybody should be able to obtain it), the participant id is `code`, the dv is `mean`:               > d <- read.table(url("http://pastebin.com/raw.php?i=4hRFyaRj"), colClasses = c(rep("factor", 6), "numeric"))     > str(d)     'data.frame':   2928 obs. of  6 variables:       $ code     : Factor w/ 183 levels "A03U","A08C",..: 1 1 1 1 1 1 1 1 1 1 ...
       $ group    : Factor w/ 2 levels "control","experimental": 2 2 2 2 2 2 2 2 2 2 ...       $ condition: Factor w/ 3 levels "alternatives",..: 3 3 3 3 3 3 3 3 3 3 ...
       $ topic    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...       $ problem  : Factor w/ 4 levels "AC","DA","MP",..: 3 4 1 2 3 4 1 2 3 4 ...
       $ mean     : num  94.5 94.5 86.5 84.5 80 46.5 73.5 43.5 51 39 ...      The usual way to analyze this data would to fit an ANOVA on this data (note how the error term is constructed for the within-subject factors):               m1 <- aov(mean ~ (condition*group*problem*topic) + Error(code/(problem*topic)), d)      ## The Question My main interest in the data is the following:   * Is there an effect of the group factor on any level (i.e., main effect or interaction)? I hope there is not.    * Is there an interaction of condition with problem ? Or even an interaction of condition with problem and topic? I have two questions regarding the analysis in `lme4`:   1. How can I specify these questions using `lme4`?   2. As `lme4` does not provide p-values, how do I determine whether a variable (e.g., group) has any effect (I imagine using some kind of likelihood ratio test) and what is the critical value above which I need to accept effect to be 'significant'? As is probably obvious from the above description I am no expert in `lme4` neither a statistician, so both Venables & Ripley and the `lme4` Book by Bates gave me a hard time. Leaving me kind of clueless as before.


title: how to fit a michaelis-menten function with a random effect using the nlme package in R?

text: I am working on fitting a model using the nlme package in R. `y` is a saturating function of `x`, similar in form to a michaelis-menten function where               y = Vm*x / (K+x)      I would like to parameterize this model using my data set, and control for a third variable, `z`, as a random effect. I cannot seem to implement this so far in the nlme framework. Any tips on how to do so, or how to do so in a different package, would be greatly appreciated.


title: SAS: When using "by" statement, what does "var1 = ." mean in output?

text: Let's say I have some simple code like what is below.                proc sort; by  var1 var2;      run;           proc gplot;      plot y*x ;      by var1 var2;      run;      In the output, there are a bunch of graphs. In addition to all the combinations of values of var1 and var2, I also get charts where `var1=.` in combination with each value of var2 (etc). Does the `var1=.` mean the plot is including all values of var1? Or data where var1 is somehow blank? Or something else? Is there a way I can get these to not print out (i.e. if I want exactly as many plots as there are combinations of var1 and var2)? It's not a big deal for only two variables, but creates a lot of extra plots to scroll through when using more variables. Thanks!


title: Variance of plugin estimator

text: This question related to my previous question. Let $$X_1,\dots,X_n$$ are i.i.d. with distribution function $F$ and $$Y_1,\dots,Y_n$$ are i.i.d. with distribution function $G$. Suppose that there exists an unknown increasing function $\psi:\mathbb{R}\mapsto\mathbb{R}$ such that $\psi(X_i)\sim N(0,1)$ and $\psi(Y_j)\sim N(0,\sigma^2)$ for all $i,j=1,\dots,n$. I wish to evaluate the variance of plug-in estimator of $\sigma^2$ in this problem. I have obtained the following facts: Note that $F(x)=P(X\le x)=P(\psi(X)\le \psi(x))=\Phi(\psi(x))$ where $\Phi$ is the cumulative standard normal distribution. This implies $\psi(x)=\Phi^{-1}(F(x))$. Recall that $\psi(Y_j)\sim N(0,\sigma^2)$. So, $$\check\sigma^2=\frac1n\sum_{j=1}^n\psi^2(Y_j)$$ is an optimal estimator for $\sigma^2$. Since $F$ is unknown then I replace it with its empirical distribution function $\hat F_n$ based on $X_1,\dots,X_n$. Hence, it is natural to replace $\psi$ with $\hat\psi=\Phi^{-1}(\hat F_n)$. Therefore, I conjecture that $$\hat\sigma^2=\frac1n\sum_{j=1}^n\hat\psi^2(Y_j)$$ is the estimator that I wanted. The problem is I don't know what the variance of $\hat\sigma^2$ is. Any suggestion? * * *


title: Kernel density estimation for a variable with lots of zeros

text: I am trying to estimate the kernel density for number of days a child is sick. Around 73% of children report not being sick, i.e. zero. How do I estimate a kernel density for this censored variable using Stata?


Sequential testing is common for Time series. Think SupF kind of tests, or a chow-test for structural break in the parameters over a grid of data points. So, we have many (as many as data points modulo few points from the start and end) statistics, and their correspondent P-values. FDR is good for the independent case, what I just described yields highly correlated test statistics, so the FDR is too restrictive. That is, a rejection is a clear sign for a break, but no-rejection.. might be because we have many data points (so FDR_pv is small), I can search over more coarse grid to reduce the number of tests, so a max over a grid = {1,2,3,...,T} may not be rejected but a max over a grid = {1,3,5,...,T} may be rejected, which does not make too much sense. Is anyone familiar, and may be able to give some references for FDR extension to the Time Series domain?


title: Calculating standard deviation of the output of a numerical experiment

text: I am simulating a game of Roulette (in `R`), where I simulate a gambler that starts with £100, spins the wheel 10 times according to some strategy (e.g. bet £5 each spin on red). At the end of this simulation the bank balance will be either higher or lower than £100. If higher, the spin is considered a win, otherwise a loss. By running the above 100 times, I can calculate the proportion of times where a win occurred, and present this as a number between 0 and 1. My question is how can I apply a standard deviation to this number? Do I need to simulate the above (simulations of simulations) many more times to do so? Because that would take a while...


title: Pitfalls in time series analysis

text: I am just starting out self-learning in time series analysis. I've noticed that there are a number of potential pitfalls that aren't applicable to general statistics. So, building on What are common statistical sins?, I'd like to ask: **What are common pitfalls or statistical sins in time series analysis?** This is intended as a community wiki, one concept per answer, and please, no repetition of more general statistical pitfalls that are (or should be) listed at What are common statistical sins?


title: Cut-off point in a ROC curve. Is there a simple function?

text: I want to find the cut-off point for gender based on an anthropological measurement. I can draw the curves and I know that in case sensitivity and specificity are both similarly important, the point closest to the upper left corner of the frame (or if the curve is negative, the closest point to the lower right corner) should be determined as cut-off. However, I don't know whether there is any already implemented function in R or any other programs for determining this, simply? I just know SPSS doesn't have such a function. Do you know any programs, or R which can do it via an already implemented function?


title: How to find the distribution of the result of a compound experiment

text: I'm trying to find the distribution of data collected from a die-roll and coin-toss experiment. The experiment is as follows: 1)Roll a fair die so that you get a number $D \in (1,...,6)$   2)Flip a fair coin $D$ number of times (i.e. the value on the dice) and record the number of tails $T$ Repeat this experiment for 10000 runs and define the random variable $W$ as the number of runs needed until an outcome of $T \geq 5$ occurs for the first time. Am I right in thinking that $T$ would be geometrically or negative binomially distributed? Just a note, I'm doing all my modelling and simulations in R.


title: KM curve: log rank test significant but HR crosses 1

text: Hello and thanks for this source about statistics. My Problem is about a situation where Kaplan Meier curves differ in meaning of a significant log rank test, but the hazard ratio crosses 1 when looking on the 95% CI. How to interpret this? Could the problem be too much censored data? I have some individuals with data for > 3000 days and most with data for only about 500 days. Thanks in advance an greetings from switzerland. Martin


title: SEM or path analysis: Latent vars with individual items as indicators, or composite scores as observed vars?

text: I am unsure which analysis is more appropriate out of path analysis or SEM (Structural equation modelling). I am looking at the relationships between a number of variables which have been measured by scales (i.e. Extraversion & Neuroticism subscales of the BFI, Self-Liking subscale of the SLSC-R, etc.) I believe I have two choices:   1. Conduct a SEM using each individual question/item as an indicator of the factor which the subscale represents. The scales I'm using vary from 4 items to 10 items. Measurement error is controlled but model fit might be more difficult due to greater complexity.   2. Conduct a path analysis by using the composite scores of the subscales as the observed variables and conduct a path anaysis only. Assumption of no measurement error Thanks for your advice.


title: Probability model for repeated survey data

text: I'd like to find a suitable probability model but I'm not sure if the one I've thought of makes sense. The intention is to apply a Bayesian analysis to survey data which has been sampled twice. The setting:   * Individuals are asked a series of nine questions (out of many others)   * These nine questions are yes/no   * The same individuals are asked the same nine questions one year later   * I'd like to know if there has been a change in response over time Modeling the number of yes's as $X \sim$ Binomial(9, $\theta$) seems obvious, but how should I incorporate the two sampling times? It makes sense to take the number of responses from the second sampling and subtract them from the first sampling and do inference on $\theta_{diff} = \theta_{after} - \theta_{before}$. However, the difference of Binomial r.v.s is ugly and it would be nice to have conjugacy. I could just do $X \sim$ Binomial(9, $\theta_{before}$) and $Y \sim$ Binomial(9, $\theta_{after}$) but the questions were answered by the same people. In other words, why would one perform a two-sample t-test when the information necessary for a paired t-test is available? Since $X + Y \sim$ Binomial($n_1 + n_2, \theta$) if $X \sim$ Binomial($n_1, \theta$) and $Y \sim$ Binomial($n_2, \theta$), does it make sense to model the the first batch of questions as $X \sim$ Binomial(9, $\theta_{first}$) then add the second batch to the first and model $X + Y \sim$ Binomial(18, $\theta_{both}$) and look at the difference, overlap, etc. of the posteriors of $\theta_{first}$ and $\theta_{both}$? Am I being too picky about the difference of Binomials? Am I missing something obvious? If anyone has suggestions or useful references they would be much appreciated.


title: What's the distribution of $\bar{X}^{-1}$?

text: What's the distribution of $\bar{X}^{-1}$ with X being a continuous iid random variable that is uniformly distributed? Can I use the CLT here?


title: Cosine of a uniform random variable

text: I have to find pdf of Y = cosine(X) where X is a random variable distributed uniformly in [-pi,pi]. I solved this using distribution function method, and the result was: $f_{Y}(y) = \dfrac{1}{\pi sin(cos^{-1}y)}, y\ \epsilon\ [-1,1]$ I am not able to understand this result intuitively. At first thought, it occured that cos(X) should be distributed uniformly as well, but that is clearly wrong. I completely fail to understand why the pdf of cos(X) goes to infinity at -1 and 1 as well. Please help me form an intuitive understanding of this result, I am not from a statistics background.


title: Given that one can sample $X \sim f(x)$, is there an easy way to sample $Y \sim k \cdot f(g(y))$ (such as $k \cdot f(e^y)$)?

text: Say I'm able to sample an RV $X$ from a PDF $f(x)$, can I exploit this to efficiently sample another RV $Y \sim k \cdot f(g(y))$ (where $k$ is a normalizing constant)? I'm interested in something for a Gibbs sampler that would be better than using Metropolis, slice, inverse transform sampling, etc.


title: Confusion related to Kulldorff's scan statistics

text: I was reading this paper related to Bayesian spatial scan statistics where I came across the Kulldorff's scan statistics. I have attached the screenshot of the paper. My objective is to find a location of potential disease outbreak. ![enter image description here](http://i.stack.imgur.com/ml9XB.png) I have a little doubt with the formula given in the paper. In the paper it is mentioned that the Fscore of any location is ![enter image description here](http://i.stack.imgur.com/FALCe.png) My question is when they calculate the F score. The number of sick people in current location(Cin) = 5 The baseline population is around(Bin) 50 The number of sick people outside current location(Cout) = 2245 The number of sick people outside current location(Bout) = 62475 Call = 5+2245 = 2250 Ball = 50+62475 = 62525 Now if I calculate the Fscore it comes out to be Nan. Is it the correct formula? Any insights. In the next link I found something like this ![enter image description here](http://i.stack.imgur.com/qyqQE.png) Now even though I use this formula, the value is NaN. It is because like in the case of denominator the prob < 1 and I am raising it to a high power, so it comes out to be zero. So my question is, is it the problem with my data itself? If I plug it into the formula it will definitely come out to be like that.


title: Finding most informative feature subsets given dataset, clustering algorithm and gold standard partition

text: I have an $n \times m$ matrix of data $\mathbf{D}$ as well as a $k$-partition $P$ of $n$ indices each representing a row in $\mathbf{D}$. Assuming an arbitrary clustering algorithm $A$, I would like to find a subset $F$ of $\\{1,\ldots, m\\}$, representing indices of the columns of $\mathbf{D}$, such that applying $A$ to only the columns indexed by $F$, an $n \times |F|$ matrix, minimizes the average variation of information (VI) between $P$ and the output of $A$. Another way of saying this is that I would like to be able to find the subspace on which $\mathbf{D}$ is most informative about $P$. (That the metric is VI or that the statistic taken over multiple runs is an average is not necessarily important. I'm looking for a general solution, but if a solution that specifies either the metric or the clustering algorithm exists, I'd be interested.) Attempting to build the distribution of (normalized) variation of information values by brute force is not practical since the number of possible subsets of $\\{1,\ldots, m\\}$ is the Stirling number (of the second kind) $S(m, 2)$, which obviously gets monstrous as $m$ gets to any interesting size.


title: Estimating the area between two ecdfs

text: I wish to calculate the area between two ecdfs in R (see below):               x<- 1     y <- c(0.329, 0.635, 0.558, 0.397, 0.719, 0.657, 0.356, 0.685, 0.618)     plot(ecdf(x), do.points = TRUE, verticals=T, xlim=c(0,1))     lines(ecdf(y), lty=3, verticals=T)      A hint in the right direction will be much appreciated. The plot would look like this: ![ecdf](http://i.stack.imgur.com/jL3Pr.jpg) I wish to calculate the area between the dotted line (ecdf(x)) and the solid line (ecdf(y)).


title: Simulate/Predict soccer standings after X rounds, based on current results/standing

text: How can I simulate/predict the soccer standings after a certain number of rounds, based on the currently played maches/standings. I would like to do this with R. After the first round, the result will not be precise, but after a few rounds, the simulation would probably become closer and closer to the real result. Data example:               Round    Home team    Away team    Home score    Away score     1        Man Utd      Arsenal      3             1     1        Liverpool    Man City     0             1     1        Blackpool    Newcastle    0             0     1        QPR          Leeds        2             1     1        Chelsea      Wigan        4             0     2        Arsenal      Chelsea      2             2     2        Liverpool    Wigan        0             1     2        Blackpool    Leeds        2             2     .     .     .     .     And so on


title: Is it valid to use only some of the components as regressors from a PCA?

text: I was trying a 4 components PCA where the 1st and 2nd components have emerged exactly as I expected, but the other two components do not commensurate with the theory. I mean, I expected a little different kind of loads. Some of the variables expected to belong to component 3 have actually loaded highly on component 4 and some expected to belong to component 4 have loaded on component 3. I have tried in different ways, but things don't really improve and theoretically component 3 and component 4 do not make much sense to me. You can have a view of the loadings table here- ![enter image description here](http://i.stack.imgur.com/BoOIH.jpg) Now, I am planning to use the 1st and 2nd components as regressors (as they have clearly emerged) along with the 9 individual variables that have loaded highly on the other two components (because their linear combinations don't hold to be meaningful) in a further regression on a dependent variable. My questions are- 1) Is it statistically valid to use only component 1 and 2 (component scores, basically) along with the 9 individual variables of component 3 and 4 as regressors in a regression? 2) Is it possible to include a few other variables as regressors too which were not included in the PCA?


title: Stepwise Regression Models in JMP

text: In JMP, I am building a regression model by using "Analyze"->"Fit Model" and choosing "Stepwise" for the personality. Once I click "Run" in the "Model Specifications" window, I get the "Fit Stepwise" window that allows me to specify how I would like my model to be built. So I hit "Go", and it builds my model. My question comes in understanding the difference between what is displayed to me in "Current Estimates" and the values I get if I click on "Run Model" and look under "Parameter Estimates". It seems like these should be the same, but clearly JMP is doing something inbetween the Stepwise window and the Fit Model window. Could someone please explain to me the difference between "Current Estimates" (in the Stepwise window) and "Parameter Estimates" (in the Fit Model window)?


title: Correlated features produce strange weights in Logistic Regression

text: I have a data set with highly positively correlated features that I'm classifying with LR. AFAIK correlated weights are not a problem in the same way they are in Naive Bayes - overcounting will not occur with LR. The strange things that I'm seeing is that some of the highly correlated features assume opposite weights: feature A might be highly positive and feature B highly negative, though not as much. Is this a symptom of something going wrong with optimization or is this expected (a priori I expect A and B to be positive class indicators)


title: Can ITML overfit?

text: When training using Information Theoretic Metric Learning (ITML), should I be concerned that the algorithm will overfit the data? If yes, what is the best way to avoid it? http://www.cs.utexas.edu/~pjain/itml/ Thanks. P.S. Can someone with permission add the tags ITML and Overfitting to this question?


title: In message-passing methods, what is the actual content of the messages?

text: In message-passing methods, factors and random variables exchange messages that typically encode marginals, but as much as I look at their formulas, I still don't understand what those messages actually look like. For example, in belief propagation, when using the factor graph representation, we have the following scheme:   * Messages **from a variable** node $v$ **to a factor** node $u$ $\forall x_v\in Dom(v),\; \mu_{v \to u} (x_v) = \prod_{u^* \in N(v)\setminus\\{u\\} } \mu_{u^* \to v} (x_v).$   * Messages **from a factor** node $u$ **to a variable** node $v$ are the product of the factor with messages from all other nodes, marginalised over all variables except the one associated with $v$: $\forall x_v\in Dom(v),\; \mu_{u \to v} (x_v) = \sum_{\mathbf{x}'_u:x'_v = x_v } f_u (\mathbf{x}'_u) \prod_{v^* \in N(u) \setminus \\{v\\}} \mu_{v^* \to u} (x'_{v^*}).$ I understand that these messages represent marginals and that the eventually we end up with the marginals in every node/factor, but my question is:   * Are these messages **vectors** holding real numbers?    * What is the **dimensionality** of those vectors?   * What does this dimensionality **represent**? To make things more concrete let's assume a distribution over **real** random variables (e.g. gaussian), although I am also interested in the discrete case (e.g. categorical random variables). Also, how are these messages initialized?


title: Assumptions in linear discriminant analysis

text: Suppose we have 4 groups of animals (pig, dog, rabbit, cat). There are 10 animals in each group. We measure the weight and pulse of each animal. Does linear discriminant analysis assume that the data within **each** group is multivariate normally distributed? Or that the data as a whole is multivariate normally distributed?


title: Contribution of residuals to the mean error

text: I have data for 1000 students' performance over 10 different tests on a scale of 0 -100 (a 1000 rows X 10 col matrix). I calculated the mean score and the associated std. deviation for each student. I wish to find out if the student(s) are performing low on a particular test. I assume that the if such a test exists, then its weightage in the overall mean score will be higher relative to other tests. The question is how do I determine the influence of individual test on the mean score of each student ? Anybody ? Thanks.


title: Why is $R^2$ poor for AR model selection used for forecasting?

text: There is a related question here, about how to calculate the R-squared on a regression with ARIMA errors. I found the answer quite useful, and hoped for some elaboration, particularly on Rob's closing remark: > Once you move away from the usual regression set up with an intercept and > iid errors, $R^2$ ceases to be uniquely defined **and is not particularly > useful**. Suppose I am trying to select a model with the objective of making reasonable k-step ahead forecasts. I also have my choice of exogenous variables (including data for k-steps ahead). For example, I have observations on $t \in [0,50]$ for Y, X1, X2... I also have a scenario analysis with forecasted data for on $t \in [51,60]$ for X1, X2... The response and exogenous variables are clearly integrated. As a type of experiment, I examine two different models. One is given by `Y ~ L(Y,1)`, and another by `Y~X1`. It so happens that the `Y~L(Y,1)` model has insanely high $R^2$, but I know that this model has much worse predictive power than the `Y~X1` model. How do I justify that $R^2$ should not be used to rank these models? The best I can think of is just holding out some data, and performing out-of- sample tests to demonstrate that `Y~X1` is better than `Y~L(Y,1)`, but this feels anecdotal to me. Is there some more rigorous proof? Thank you!


title: Effects Coding in R

text: I am in interested in how do effect coding in R. I know that someone else has asked this question (i.e How to do regression with effect coding instead of dummy coding in R?). Here is the lm() model on its own:                 Fixed.Full<-lm(sqrtClay_Tot ~ Physio.Code+ROCKTYPE1, data=PedonsTx_Reread,                  contrasts =  list(Factor.0.Physio = contr.sum))          > summary(Fixed.Full)          Call:     lm(formula = sqrtClay_Tot ~ Physio.Code + ROCKTYPE1, data = PedonsTx_Reread)          Residuals:         Min      1Q  Median      3Q     Max      -4.5719 -1.1543  0.0102  1.1366  6.0146           Coefficients:                                         Estimate Std. Error t value Pr(>|t|)         (Intercept)                          5.18961    0.54194   9.576  < 2e-16 ***     Physio.Code2                        -0.90588    0.15752  -5.751 1.11e-08 ***     Physio.Code3                         0.70532    0.21141   3.336 0.000874 ***     Physio.Code4                         0.56462    0.31295   1.804 0.071444 .       Physio.Code7                        -0.04656    0.26410  -0.176 0.860075         Physio.Code8                        -0.41242    0.24330  -1.695 0.090294 .       Physio.Code9                        -0.81833    0.32048  -2.553 0.010782 *       Physio.Code11                       -0.04619    0.16420  -0.281 0.778515         ROCKTYPE1chert                       1.43390    1.80191   0.796 0.426319         ROCKTYPE1clay or mud                 0.35031    0.53861   0.650 0.515562         ROCKTYPE1evaporite                  -0.95719    0.65049  -1.471 0.141408         ROCKTYPE1fine-grained mixed clastic -1.02193    0.59167  -1.727 0.084377 .       ROCKTYPE1gravel                     -0.16708    0.56139  -0.298 0.766041         ROCKTYPE1limestone                   0.36667    0.55323   0.663 0.507585         ROCKTYPE1mudstone                    0.18934    0.58553   0.323 0.746476         ROCKTYPE1rhyolite                   -0.25743    0.83684  -0.308 0.758426         ROCKTYPE1sand                       -0.92281    0.53193  -1.735 0.083014 .       ROCKTYPE1sandstone                  -1.57669    0.57599  -2.737 0.006281 **      ROCKTYPE1serpentinite               -0.18647    1.84384  -0.101 0.919464         ROCKTYPE1shale                       0.34408    0.62577   0.550 0.582518         ROCKTYPE1silt                       -0.82415    0.65554  -1.257 0.208912         ROCKTYPE1siltstone                  -1.21094    0.79293  -1.527 0.126970         ROCKTYPE1terrace                    -0.64082    0.60119  -1.066 0.286663         ROCKTYPE1tuff                        1.10034    1.31108   0.839 0.401482         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 1.748 on 1261 degrees of freedom     Multiple R-squared:  0.2072,    Adjusted R-squared:  0.1928      F-statistic: 14.33 on 23 and 1261 DF,  p-value: < 2.2e-16      Now here is the model using effecting coding within lm():               > Fixed.Full<-lm(sqrtClay_Tot ~ Physio.Code+ROCKTYPE1, data=PedonsTx_Reread,     + contrasts =  list(Physio.Code = contr.sum))     > summary(Fixed.Full)          Call:     lm(formula = sqrtClay_Tot ~ Physio.Code + ROCKTYPE1, data = PedonsTx_Reread,          contrasts = list(Physio.Code = contr.sum))          Residuals:         Min      1Q  Median      3Q     Max      -4.5719 -1.1543  0.0102  1.1366  6.0146           Coefficients:                                         Estimate Std. Error t value Pr(>|t|)         (Intercept)                          5.06968    0.50933   9.954  < 2e-16 ***     Physio.Code1                         0.11993    0.13291   0.902  0.36703         Physio.Code2                        -0.78595    0.11689  -6.724 2.68e-11 ***     Physio.Code3                         0.82525    0.15844   5.208 2.22e-07 ***     Physio.Code4                         0.68455    0.24897   2.750  0.00605 **      Physio.Code5                         0.07337    0.21084   0.348  0.72792         Physio.Code6                        -0.29249    0.18905  -1.547  0.12208         Physio.Code7                        -0.69840    0.26157  -2.670  0.00768 **      ROCKTYPE1chert                       1.43390    1.80191   0.796  0.42632         ROCKTYPE1clay or mud                 0.35031    0.53861   0.650  0.51556         ROCKTYPE1evaporite                  -0.95719    0.65049  -1.471  0.14141         ROCKTYPE1fine-grained mixed clastic -1.02193    0.59167  -1.727  0.08438 .       ROCKTYPE1gravel                     -0.16708    0.56139  -0.298  0.76604         ROCKTYPE1limestone                   0.36667    0.55323   0.663  0.50759         ROCKTYPE1mudstone                    0.18934    0.58553   0.323  0.74648         ROCKTYPE1rhyolite                   -0.25743    0.83684  -0.308  0.75843         ROCKTYPE1sand                       -0.92281    0.53193  -1.735  0.08301 .       ROCKTYPE1sandstone                  -1.57669    0.57599  -2.737  0.00628 **      ROCKTYPE1serpentinite               -0.18647    1.84384  -0.101  0.91946         ROCKTYPE1shale                       0.34408    0.62577   0.550  0.58252         ROCKTYPE1silt                       -0.82415    0.65554  -1.257  0.20891         ROCKTYPE1siltstone                  -1.21094    0.79293  -1.527  0.12697         ROCKTYPE1terrace                    -0.64082    0.60119  -1.066  0.28666         ROCKTYPE1tuff                        1.10034    1.31108   0.839  0.40148         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 1.748 on 1261 degrees of freedom     Multiple R-squared:  0.2072,    Adjusted R-squared:  0.1928      F-statistic: 14.33 on 23 and 1261 DF,  p-value: < 2.2e-16      The code above is functional within the lm() function as seen above,and what I am trying to do in the code above is to predict clay content (square-root transformed) from the categorical variables 'Physio.Code' and 'ROCKTYPE1'. What I would like to know, however, is if there is any way to do effect coding using the gls() or glm() functions in R, where method="REML" or REML estimation can be used. I am interested in comparing the levels within 'Physio.Code' specifically to to the grand sample mean (i.e. sum contrasts). The other question I had was how to interpret the numbers after the categorical variable and whether they correspond to the same levels. With respect to 'Physio.Code', in the two models above '1' was used as the baseline - the levels were originally '1,2,3,4,7,8,9,11" but when '1' was used as the baseline group, it switched to "2,3,4,7,8,9" - I wasn't sure if levels are the same (i.e. Physio.Code2 in model is the same as Physio.Code1 in the 2nd model). Thank you.


title: Comparing migration statistics of the same city

text: I've got a city with 12 districts and I've got a table of the number the people who moved within the city to other districts, which looks roughly like this, but the numbers are just for illustration:               to\from 1   2   3   4   ... 12     1       -   100 200 100     300     2       200 -   100 300     100     3       250 100 -   100     300     ...     12      200 100 200 100     -      Also, I have a different source which also maps the migration process somehow, but they have measured a smaller number of overall mitigations, because they didn't extrapolate their measurements to the whole population, e.g. like this:               to\from 1   2   3   4   ... 12     1       -   1   2   2       2     2       2   -   1   2       1     3       2   1   -   1       2     ...     12      2   2   1   1       -      Which tests are useful to compare these kind of measurements? I think that looking at a PCA of these tables (and their transpositions) might give an idea if these measurements are different, but I'm not sure which kind of test should be used to compare these kind of measurements. Would a Chi-Squared test work in this setting?


title: Data tranformation in panel data

text: I am doing panel data analysis and some of my variables have high kurtosis. I am not sure whether I have to transform these variables. I have tried to delete outliers but one of the variables still not normal unless I delete many observations which then change the results to insignificant model note: this variables has a maximum of around 17,000,000 and min of - 5000,000 which i couldn't use log transform any help is appreciated. thanks


title: Restricting model parameters in logistic models in R

text: Is there any function in `R` that can solve the problem like this example from the SAS website: > Beginning in SAS 9.3, PROC FMM can be used as an alternative to the LOGISTIC > and GENMOD procedures for fitting generalized linear models such as logistic > and poisson models. You can fit the model in PROC FMM and use its RESTRICT > statement to impose equality or inequality constraints on the model > parameters. > > For example, in the following logistic model suppose you want to constrain > the parameters for X1 and X2 to be equal. >   >   >      proc logistic; >         model y = x1 x2 x3 x4; >         run; >   > > The following statements fit the model in PROC FMM and impose the > restriction. >   >   >      proc fmm; >         model y = x1 x2 x3 x4 / dist=binary link=logit; >         restrict x1 1 x2 -1; >         run; >   > > To restrict the parameter on X1 to exceed that of X2, use the following > RESTRICT statement. >   >   >      restrict x1 1 x2 -1 > 0; >


title: Scikit linear model predictions not matching

text: I'm using scikits linear model for modelling two input variables and one output variable. I suspect that the two input variables have a quadratic relationship with the output variable. So I created the data set ($a^2$, $b^2$, $ab$, $a$, $b$, $1$) for the two original input variables $a$ and $b$, then used the scikits linear model for this data which gave the coefficients:               [ 0.21746023  0.3457555  -0.15043191 -0.15318758 -0.27072946  0.        ]      the .predict() method produces very good predictions however, when I run the code               for i in data:         print i[0]*coeffs[0]+i[1]*coeffs[1]+i[2]*coeffs[2]+i[3]*coeffs[3]+i[4]*coeffs[4]+i[5]*coeffs[5]      The outputs do not match the predictions at all. Here is my full code:               def import_data():         data=np.genfromtxt("input_file")         y=data[:,2]         all_data=[(a**2,b**2,a*b,a,b,1) for (a,b,c) in data]         return np.array(all_data),y          def main():         data,y=import_data()              regr = linear_model.LinearRegression()         regr.fit(data, y)         coeffs=regr.coef_         print coeffs         for i in data:             print i[0]*coeffs[0]+i[1]*coeffs[1]+i[2]*coeffs[2]+i[3]*coeffs[3]+i[4]*coeffs[4]+i[5]*coeffs[5]         print regr.score(data, y)         print regr.predict(data)      Could someone tell me where my code is wrong?


title: Moment generating functions question

text: > Suppose that $X_1, X_2, ..., X_n$ are independent, where each $X_i$ has > probability (mass) function $p_i(x_i)$ given as $p_i(x_i) = > \frac{e^{-\lambda} \lambda_i^{x_i}}{x_i!}$ (only the parameter $\lambda_i$ > differs int he distribution of each $X_i$ for $x_i = 0, 1, ...$. What is the > distribution of their sum: $\sum\limits_{i=1}^n X_i$ ? Prove it (perhaps > with moment generating functions). Answer: The moment generating function of Poisson (sum of $\lambda$) I don't get how to do this question and I don't really understand the question. Can someone please help me? Thanks in advance.


title: Do low-discrepancy sequences work in discrete spaces?

text: Low-discrepancy sequences in a real space ($[0,1]^n$) seem like a really excellent tool for evenly sampling a sample space. As far as I can tell, they generalise well to any real space, if you use an appropriate map (eg. $[0,1]\to[a,b]$ linear map). Do such sequences generalise to discrete spaces? eg. if I have a space that has only two elements in each dimension (eg. boolean switches), can I just map $[0,0.5]\to 0;\ (0.5,1]\to 1$? What about for dimensions with more elements (eg. a 4-state switch?). And for spaces with different number of states in each dimension? My intuition says that this might work ok, especially if the sub sequence is longish, but that it might work better for some sequences than for others, depending on the number of states (eg. a Halton sequence might have odd interactions with a dimension with a prime number of states, or a Sobol' sequence might only work for dimensions with $2^n$ elements). But I have done no testing. If this doesn't work, why not?


title: "Moderation" versus "interaction"?

text: I have come across these two terms which are used interchangeably in many contexts. Basically, a moderator (M) is a factor that impacts on the relationship between X and Y. Moderation analysis is usually done using a regression model. For example, gender (M) can affect the relationship between "product research" (X) and "product purchase" (Y). In interaction, X1 and X2 interact to influence Y. The same example here is that "product research" (X1) is affected by "gender" (X2) and together they affect "product purchase" (Y). I can see that in moderation, M affacts the X-Y relationship but in interaction, M (which is gender in this case) affects the other IV. **Question** : If the aim of my project is to see how gender affects the relationship between X and Y, should I use moderation or interaction? **Note:** My project is about the correlation between X and Y, not causal relationship between X and Y.


title: How to test for association between proportions with multiple variables?

text: I have tested 100 bacterial isolates for resistance to 6 antibiotics (tetracycline, erythromycin, kanamycin, streptomycin, gentamicin, chloramphenicol). The result for any antibiotic could be either sensitive (1) or resistant (0). How can I determine if there is any significant association between resistance to various antibiotics (I can frequently find isolates resistant to three or more antibiotics but how to know if this is significant)? I could use Chi- square to compare between any two antibiotics but how to make multiple comparisons and will this make any difference except saving time. Scientists in the papers used Chi-square with Bonferroni adjustment. I am using Statistica for Windows software. thanks


title: Is there any way to fit an exponential model with a discrete predictor in R?

text: I'm trying to estimate the growth rates of a number of ant colonies. I have measured colony mass at three different time points. I was thinking about measuring growth rate by fitting an exponential growth model to each colony's masses across time. It was drawn to my attention that I may not be able to perform non-linear least squares (nls) as I had hoped, because the predictor is a discrete variable (i.e., it's t1, t2, or t3). I was wondering:   1. whether this is indeed problematic &   2. whether anyone had any recommendations for how else to estimate growth rates in this case?


title: Simulation of KS-test with estimated parameters

text: We know that the original KS-test has a limitation, requiring the Null hypothesis testing underlying distribution to be fully specified, rather than estimated. But in practice, we usually need to test goodness of fit for the fitted distribution. It is suggested to conduct a Monte Carlo Simulation by Wikipedia for KS-test. But there are very few limited literature review about such simulation. There are several ideas come up in my mind. * * * ## Approach 1 It seems that the reason why original KS-test is no longer valid for testing estimated distribution, is that we need a smaller critical threshold KS Distance (d value) to reject the null hypothesis underlying distribution, as it is fitted around the data. The Lilliefors-test solves such problem. Although it is only for Normal distribution, we can extend its method.   1. Find the KS Distance ($D_{0}$) between the data (n data points) and the estimated underlying distribution, like original KS-test.   2. Generate n data points by the estimated underlying distribution for M times (around millions), find their KS Distance $D_{m}$ set = {$D_{1}$, $D_{2}$, $D_{3}$, ..., $D_{M}$} for each time.   3. Compare the $D_{0}$ value and the percentile of $D_{m}$ set as the p-value. $P(>D_{0}|Null Hypothesis Distribution)$ * * * ## Approach 2 Some post said that the original KS-test can not be applied to estimated distribution, because it does not take into account the standard error of estimation (or fitting). So they suggest bootstrap from raw data as a simulation, but they have not described the details. I think the approach can be as follow.   1. Find the KS Distance ($D_{0}$) between the data (n data points) and the estimated underlying distribution A, like original KS-test.   2. Bootstrap from raw data and fit distribution again, getting a new fitted distribution B. Find the KS Distance ($D_{1}$) between this new raw sample and new distribution B.   3. Repeat step 2 for M time, getting their KS Distance $D_{m}$ set = {$D_{1}$, $D_{2}$, $D_{3}$, ..., $D_{M}$} for each time.   4. Compare the $D_{0}$ value and the percentile of $D_{m}$ set, to see whether it is out of the 95% empirical interval. It seems that this takes into account the estimated error (or confidence interval) about the estimation, which affecting the original KS-test. But it seems that this approach does not follow the p-value definition, or hypothesis definition, assuming Null hypothesis and find the sample statistics probability. * * * ## Approach 3 My fellow also comes up an idea, like the approach 2.   1. Find the KS Distance ($D_{0}$) between the data (n data points) and the estimated underlying distribution A, like original KS-test.   2. Set the Null hypothesis that the data is from the distribution in question (A).   3. Generate n points from the distribution in question (A) and fit the same distribution to this new sample which gives parameter set B. Measure the KS distance $D_{1}$ between this sample and the distribution in question with parameter set B.    4. Repeat step 3 a large number of times (M) which gives a distribution of KS distances $D_{m}$ assuming the actual data is from the distribution in question using our fitting method.   5. If distance $D_{0}$ is greater than 95% of the distances calculated in 4 then reject the Null hypothesis that the data is from the distribution in question with parameter set A. It seems that it is something parametric bootstrap from approach 2. But it appear strange to generate such sample and fit data again. And it seems still not following the hypothesis definition. * * * Do you have any idea about these three approaches? or any modification? Or do you find any formal literature saying simulation details for KS-test for estimated distribution? I only find Durbin (1973) with details about this. But sorry I can not find the full text about this book. Does anyone have any other electronic version of similar details about simulation KS-test? Does Anderson-Darling test also suffer from this estimated parameter issue? Any ideas or discussion are highly appreciated. Thank you.


title: The most common extracted features for image recognition

text: I want to build a neural network, but because I have high resolution pictures, I rejected the idea of passing the entire image to the NN. I was wondering what are the _most common extracted features_ for use as input to the NN. I know that features depends on the kind of images. But I guess I'm wondering what are the _standard features_. I'd also like to know about papers or books with information and/or examples about this topic.


title: How to obtain the inverse of the F cumulative distribution based on the F cumulative distribution?

text: I am currently implementing a application that needs to obtain the inverse of the F (Fisher-Snedecor) cumulative distribution. I already have a library that contains the F distribution and I can easily obtain the cumulative distribution based on a confidence interval and degrees of freedom. How can I easly adapt this library to return the INVERSE of the F cumulative distribution?


title: Why it is good to take log on Finance data? Does it have nice properties?

text: Just like what I am asking in the title. I see nearly all the financial datas take logs before the data analysing step, Why? Dose it have nice properties?


title: Building a mixed model for multi-location split plot over 2 years in SAS

text: For the last 3 weeks I have been reading books, papers and internet sources incl. cross validated but could not find something that would work in my case. I have data from a split-plot field trial. The field trail was laid out in the same way at 2 locations over 2 years. The literature said that it would be best to do the statistics in SAS with PROC MIXED because GLM would give me some false results. So my question would be what would be the best mixed model for 2 years and 2 locations. Fixed factors are:   1. tillage as the main plot factor   2. plant as the subplot factor   3. location Random factors are:   1. block    2. year For one year and one location I have managed to get a model as it can be found in the "hitchhiker's guide to mixed models for randomized experiments". The model for one year one location would look like this.               Model yield = tillage plant tillage*plant;     Random = block block*tillage block*tillage*plant      The next step would be the model for 2 years at one location, but for this case I could not find any literature that would work for a split plot. Because every year there was a new set up for the experiment I don't think that I could consider it as repeated measurements. What should the model look like? Step 3 would be the model over 2 locations. I have read that locations can be seen as 'superblocks' but don't know how I should build the model for the 2 year 2 location split plot. I would be so happy if somebody could help or point me in the right direction where I could find some models like this. Thanks in advance :)


title: LDA vs. labeled LDA

text: I have gone through the techniques and understood the basic ideas. But I want to know which one usually is expected to work better, LDA or Labeled LDA? What are the features of the dataset that help decide amongst the two?


title: What is Mahalanobis distance, & how is it used in pattern recognition?

text: Can someone explain to me the concept of Mahalanobis distance? For example, what is the Mahalanobis distance between two points x and y, and especially, how is it interpreted for pattern recognition?


title: Mahalanobis distance on singular data

text: I have an issue which I could not solve, although I tried and I got some help on R help forums too. I am trying to calculate Mahalanobis distances on a data frame, where I have several hundreds of groups and several hundreds of variables. Whatever I do, however I subset it I get the "system is computationally singular: reciprocal condition number" error. It is clear that it is singular, but is there any way to get rid of it and run Maha? Should I forget Maha? Than what to use? I have uploaded the data file to my ftp: http://mkk.szie.hu/dep/talt/lv/CentInpDuplNoHeader.txt It is a tab delimited txt file with no headers. I was working with the (R) StatMatch Mahalanobis (also tried stats Mahalanobis) function. I have a deadline for this project (not a homework:)), and I could always use this funct, so I thought I will be able to quit the calculations short, but now I am just lost. link to previous question, sorry for crossposting, I have no idea how to migrate the previous one. http://stackoverflow.com/questions/13078909/system-is-computationally- singular-reciprocal-condition-number I would really appreciate any help. Thanks for any help


title: Specifying add and drop thresholds for stepwise regression in R

text: I am running a stepwise regression using the F test as the criterion. Is there a way to explicitly set the add and drop thresholds (alpha levels) in R? The documentation does not make it clear.


title: Justification for low/high or tertiary splits in ANOVA

text: I have two groups of participants and they all took the same measurements. When I perform a median split on one of the measurements to test my hypothesis about an interaction, the dependent variable fails Levene's test for equality of variances. When I perform a teritary split, it passes Levene's test - both when I am comparing only high and low thirds and when I compare all three thirds too. There are main effects and interaction effects in all of the ANOVAs that I have run, but only the 2x2 ANOVA fails Levene's test for significance, neither of the 3x2 ANOVAs fail. How is this possible?


title: How to decide the number of Hidden layers and number of neurons in an Artificial neural network?

text: How to decide the number of Hidden layers and number of neurons in an Artificial neural network ? I am having an artificial neural network for a rainfall prediction application which is having four neurons in input layer and two in output layer.Which methods and heuristics should I apply to choose the number of Hidden layers and number of neurons in my network?


title: A question on sampling distributions

text: I once read the following interesting distribution related question: > An insurance company identified the 10 safest and 10 least-safe cities from > 200 cities in the US based on the mean number of years drivers went between > automobile accidents. It turns out that the cities on both lists were all > smaller than the 10 largest cities in the list of 200 cities. > > How do we use the sampling distribution model of the mean to explain this > fact?


title: Standard deviation of lift

text: I am trying to calculate the standard deviation of lift given treatment and control groups, and the respective response values, which are binary, for each entity in the groups. Here is the problem in detail: My data looks like this: > Control:   >  D1, 0   >  D2, 0   >  D3, 1 > > Treatment:   >  D4, 0   >  D5, 1 Here D_i, is a data point, and the adjoining value, 0 or 1, is the response. The response rate of the control group is 33% while that of the treatment group is 50%. Hence, the lift is 50/33 = 1.5. Now, I want to know whether calculating the lift on the whole of the control/treatment population is painting things with too fat a brush? i.e.if I had divided the control and treatment groups into an equal number of smaller populations, paired these smaller control and treatment groups, and calculated the lifts on each of these paired populations, how much dispersion would the lift exhibit? I guess, given the above explanation, the naive way to go about it is as described: randomly partition the original groups into an equal number of populations, create pairs of control and treatment subsets, calculate the lifts on each of these paired populations, and then compute the standard deviation of these lifts from the mean lift. Is there a better/principled way to do this?


title: Regression test coefficient on variable is greater than coefficient on interaction term

text: With a model such as: $y \approx B_0 + B_1\cdot \log(x) + B_i\cdot \log(x):\text{group}_i +B_j\cdot group_i$, where group can take on several values ($i = 2$ to $15$, let's say): In an OLS regression, a statistically significant coefficient on any one of the interaction terms $B_2$ through $B_{15}$ means I know up to an acceptable level of Type 1 error that the coefficient is not 0. In this particular model, the theory suggests that the slope between $y$ and $\log(x)$ in every group should be negative. So I want to test the hypothesis that $B1 + B2 < 0$, for instance, indicating that I'm confident the relationship between $y$ and $\log(x)$ in group 2 is negative. How can I accomplish such a test? It would be even more helpful if you have specific guidance on how to do so in Python's statsmodels package.


title: Explain a continuous variable by other continuous variables

text: I'm trying to analyse the correlation between data extracted from a running tracking app. From this dataset, I would like first to study how some continuous variables (amongst elevation, heart rate, session hour, pace, etc) are influencing the speed/pace. Which method should I use to analyse this relation ? Linear regression with each variable and compute the correlation coefficient ? Multiple linear regression ?


title: Is Student's $t$-test the right choice?

text: I have 2 patient populations taken from the same time period that underwent 2 different surgical laparoscopic procedures. I want to compare the rate of conversion to an open surgical procedure. In my data this is represented as dummy code under the variable OpenConversion. Can I do a _t_ -test to compare the rates, or am I choosing the wrong test? Thanks!


title: Number of baseline hazard functions

text: Suppose we have a Cox proportional hazard model and we stratify by group. Note that there are three groups. Does this mean that we would estimate three baseline hazard functions? If we just adjusted for group instead of stratifying would we just estimate one baseline hazard function?


title: Is there a difference between administrative censoring and right censoring?

text: It seems that administrative censoring is when something happens that ends the study (e.g. lack of money). People who have not experienced the event during the follow up period are right censored. Is there a difference between the two?


title: How to 'read' (understand ) an expected value equation (example inside)

text: I have just come across expected values and they are giving me a bit of grief trying to understand them. e.g. for covariance the equation is $\text{E}\left((x - \bar{x})(y - \bar{y})\right)$ To me this means the expected value is take my x value , subtract the mean, and then multiply it by the y value minus the y mean. What I don't get is what in this equation says to do this for every pair of xy values i have? I cant see a symbol in it that says to sum everything up ... $\Sigma$. How do we know to sum everything up ? Also looking at the using a 4 sided dice, lets look at the expected value of our roll. E = 1/4 * 1 + 1/4 * 2 + 1/4 * 3 + 1/4 * 4 = 2.5 Here I have multiplied each possible value by its probability of occurring. This is what I understood an Expected value to be, but I don't get how this relates to $\text{E}((x - \bar x)(y -\bar y))$ ...??


title: Bayesian hurdle model and predictions

text: I fitted a Bayesian hurdle model assuming a Poisson distribution: $ P(Y_i=y_i) = \begin{cases} 1-\pi_{i} & \mbox{if }y_i=0\\\\[0.5em] \pi_{i} \frac{\lambda_{i}^{y_i} e^{-\lambda_i}}{y_i!(1-e^{-\lambda_i})} &\mbox{if }y_i > 0 \end{cases} $ Because the log-likelihood is separable with respect to $\pi$ and $\lambda$, I fitted two separate regressions (i.e., logistic and zero-truncated Poisson regression) using a set of predictors. In other words: $ \lambda = \log(\beta_{1}^T X) \\\ \pi = \mbox{logit}(\beta_{2}^T X) $ My issue arises in trying to predict new values. In a Bayesian context, I have posterior distributions for the coefficients from the logistic and zero- truncated Poisson regressions. For a new predicted value, I will have a posterior distribution for $\pi$ and $\lambda$. To get the $E(\pi \cdot \lambda)$ is easy, I can just take the product of the posterior means. But what about getting the 2.5 and 97.5% percentiles for the credible interval around for the new predicted value?


title: Relationship between momentum and minibatch size, should they be used together?

text: It seems to me the both momentum and using minibatchs achieves almost the same thing. They both add a form of relating the changes to the weights from other test cases. In minibatch, the weight changes are summed over the batch, then applied before the next. In momentum, the change from this test 1 is averaged (generally with a strong weighting to more recent) with past changes, then applied. They don't quiet do the same thing, but it seems similar.   **Do they actually have fundamentally different effect?   ** Should they be used together?   Should the minibatch size be smaller if using both? * * * 1 Or minibatch of tests if using both


title: Simple Probability of an event

text: I have 1000 computers and 2 have failed with a specific code bug. They all are exposed to the bug, because they all have the same code running. Just using those numbers, is the probability of the failure event 2/1000 ?


title: $X_{n}$ and $Z$ be random variables if $X_{n} \ge Z$ then $ E[\liminf_{n\to \infty} X_{n}] \le \liminf_{n\to \infty} E[X_{n}] $

text: Let $X_{n}$ and $Z$ be random variables on probability space $(\Omega ,\mathcal F,P)$ and $Z$ be integrable. show that $$X_{n} \ge Z \qquad \Longrightarrow \qquad E[\liminf_{n\to \infty} X_{n}] \le \liminf_{n\to \infty} E[X_{n}] $$


title: When to use Student's or Normal distribution in linear regression?

text: I am looking at some problems, and in some, to test the coefficients, sometimes I see people using Student's distribution, and sometimes I see Normal distribution. What is the rule?


title: Testing a binomial model

text: I have a computational experiment that takes a size parameter $s$ and returns true or false. It has been suggested that the result of the experiment is independent of $s$, but I don't think so. I expect a general decrease in the numbers, though only logarithmically dependent on $s$. (This makes it hard to detect and is the essential reason for my posting -- a larger effect would be easier to detect.) For s = 5..30 I have generated 1000 samples for which I give the number of occurrences, respectively: 108, 108, 99, 99, 98, 87, 82, 85, 84, 94, 97, 87, 86, 102, 82, 84, 83, 87, 75, 69, 66, 79, 91, 81, 74, 74 Another wrinkle: the time needed to generate the samples depends on the size parameter. Making 1000 tests at size 10 took less than a second; at size 30 it takes nearly a minute. So I could generate many more examples at small sizes if that would be helpful. ($s$ relates to $10^s$, so I can't pick too many more than $10^s$ for a given $s$. But I could run a billion tests at s = 10 without trouble [in a week or two].) The flip side, of course, is that it's difficult to generate large examples, even though those are the ones that would be most useful. With a bit more effort (half a dozen processor-hours) I was able to do 1000 tests at higher levels:   s = 40: 75   s = 45: 81   s = 50: 74   s = 55: 73   s = 60: 74   s = 65: 68 So my question: How should I design an experiment to decide if these are identically distributed? Once I have my data, how can I decide (for some $\alpha$) whether to reject or not reject the hypothesis that the distributions are identical, i.e., independent of $s$? For example, I could continue to generate data for 1000 tests at each value of $s$ from 31 to 70. But given the increasing time needed as the parameter increases I'm not sure that this is an appropriate way to sample the data. But once I introduce non-uniformity I'm not sure what testing technique would be appropriate. (I'm a statistical neophyte.)


title: Using k-fold cross-validation for time-series model selection

text: **Question:** I want to be sure of something, is the use of k-fold cross- validation with time-series is straightforward, or one need special attention before using it? **Background:** I'm modeling a time series of 6 year (with semi-markov chain), with a data every 5 min. To compare several models I'm using a 6-fold cross- validation by separating the data in 6 year, so my training sets (to calculate the parameters) have a length of 5 years, and the test sets have a length of 1 year. I'm not taking into account the time order, so my different sets are :   * fold 1 : training [1 2 3 4 5], test [6]   * fold 2 : training [1 2 3 4 6], test [5]   * fold 3 : training [1 2 3 5 6], test [4]   * fold 4 : training [1 2 4 5 6], test [3]   * fold 5 : training [1 3 4 5 5], test [2]   * fold 6 : training [2 3 4 5 6], test [1]. I'm making the hypothesis that each year are independant from eachother. How can I verify that? Is there any reference showing the applicability of k-fold cross-validation with time series.


title: Confidence intervals and central estimates for a functional of an estimated function with uncertain parameters

text: I've got a problem that is leading me to dip my toes into Bayesian stats, and I've got a question about confidence (or, I suppose, credible) intervals: Say you want to know how $X$ maps to $y$. You fit a model $y=f(X)+\epsilon$. Then, you want to optimize $X$ to get the best $y$: $$y_{max} = argmax_X(\hat{f}(X),s.t. \text{whatever constraints}) $$ This gives you the model's best estimate of the optimal $X$ for getting the biggest $y$. But obviously $\hat{f}(X)$ is uncertain. If you take a Bayesian standpoint that $\beta$ is distributed multivariate normal, you can take samples from it, which gives new coefficients (see, for example, this). Taking many samples, using them to pick new optimal values of $X$, one gets a distribution of $y_{max}$ that reflects uncertainty in $\hat{f}(X)$. Here is the problem: the central estimate of $y_{max}$ (i.e.: optimizing based on the parameter estimates of the fitted model) is not necessarily the mean or the median of the $y_{max}$ distribution that one gets when optimizing the functions based on the posterior draws. So what should I do with the "central estimate"? Which estimate should I consider to be my "best guess" of the value of $y_{max}$? Should it be $y_{max}$ at the (ML) parameter estimates? Should it be the mean or median of the posterior simulations of $y_{max}$? I don't know whether there is a right answer here: maybe this is a somewhat of a philosophical question? (Or am I making some relatively fundamental mistake, which makes my whole question moot? If so, I'd be grateful for replies that point it out.)


title: Dependent variable selection for loglinear segmented regression in time-series analysis of rare events

text: My office is going to implement a bundle of infection control measures in hospital and see if it can effectively reduce the infection rate of some pathogen. The unit of measurement will be "case per thousand patient bed days". We have chosen 4 wards for implementing the control measures for 12 months, and do the measurement monthly, but by looking at the current infection data, even thses wards are considered as having the highest count of infections, they are still considered as relatively clean, as they have several months with zero measurement in rates. I have created a regression model showing if they managed to reduce 50% of infection, and the betas in the model (with three variables inside), are not statistically significant. My colleagues are worried that with all these hard work, giving non-significant result will be very frustrating for the front line staff. Is there any alernative outcome measure, or even alternative statistical methods, for time series of rare events? Thanks!


title: Binary classification with too few positive samples

text: I met a problem of doing a binary classification with quite few positive samples. For example:   * Binary classification with either labelled 1 samples (positive) or labelled 0 samples (negative samples).   * 110 cases: 10 positive samples labelled with 1 and 100 negative samples labelled with 0.   * Classifier like SVM classified 5 positive as positive but mis-classified the other 5.    * The accuracy is very high, but low AUC. Is there a better way than doing leave-one-out cross validation? Can I split the 100 negative samples to 10*10 sets and do 10 leave-one-out experiment with 10 neg+ 10 pos samples and average the performance? Or if you have a better idea please tell me. Thanks very much. A.


title: Timing functions in R

text: 1. I would like to measure the time that it takes to repeat the running of a function. Are `replicate()` and using for-loops equivalent? For example:              system.time(replicate(1000, f()));     system.time(for(i in 1:1000){f()});      Which is the prefered method.   2. In the output of `system.time()`, is `sys+user` the actual CPU time for running the program? Is `elapsed` a good measure of time performance of the program?


title: Time series 'clustering' in R

text: I have a set of time series data. Each series covers the same period, although the actual dates in each time series may not all 'line up' exactly. That is to say, if the Time series were to be read into a 2D matrix, it would look something like this:               date     T1   T2   T3 .... TN     1/1/01   100  59   42      N/A     2/1/01   120  29   N/A     42.5     3/1/01   110  N/A  12      36.82     4/1/01   N/A  59   40      61.82     5/1/01    05  99   42      23.68     ...     31/12/01  100  59   42     N/A          etc       I want to write an R script that will segregate the time series {T1, T2, ... TN} into 'families' where a family is defined as a set of series which "tend to move in sympathy" with each other. For the 'clustering' part, I will need to select/define a kind of distance measure. I am not quite sure how to go about this, since I am dealing with time series, and a pair of series that may move in sympathy over one interval, may not do so in a subsequent interval. I am sure there are far more experienced/clever people than me on here, so I would be grateful for any suggestions, ideas on what algorithm/heuristic to use for the distance measure and how to use that in clustering the time series. My guess is that there is NOT an established robust statistic method for doing this, so I would be very interested to see how people approach/solve this problem - thinking like a statistician.


title: What statistic to use

text: My 8 year old daughter is conducting an experiment to assess memory. She has 3 groups. `Condition 1`, `Condition 2` and `Control`. Does she simply run a t-test for `group 1` against `Control` and then another t-test for `group 2` against `Control`, then a t-test for `group 1` against `group 2`?


title: How to do a repeated measures analysis with binary DVs and a single binary IV?

text: I ran an experiment where participants were randomly allocated to 2 conditions. They then completed 3 seperate tests with a binary pass/fail outcome for each test. So I have 3 binary (within-subject) test scores for each participant and 2 (between-subject) groups of participants. I need to determine:   1. Were pass rates were equivalent for the 3 tests?   2. Did the participant's condition influence pass rates overall?   3. Did condition interact with any particular test outcome? If the DVs were continous this would be a straightforward repeated measures, mixed model ANOVA. I have been researching what I should do. I came across Cochran's Q test which seems a perfect way of determining if pass rates varied between the 3 tests but there does not seem to be any way of accounting for the two experimental conditions. Is there a similar test that allows for an additional IV? I'm using SPSS and have also looked at Generalized Estimating Equations (GEE), which seems like it might be helpful but it appears that there it is only possible to enter a single DV. If I use this procedure how do I enter the 3 repeated test scores?


title: Can I apply a t-test to normalized data?

text: I'm conducting biological experiments and have to normalize my raw data to the value of the untreated control in order to analyze them. This means that I have for the untreated sample something like {1,1,1,1,1} and for the treated sample something like {1.12,2.1,2.05,1.85,1.78}. How do I test if both sets are significantly different? I obviously cannot enter 1,1,1,1,1 into a t-test, can I? In case I want to test for a significant difference in two treated samples, lets say treated with two different doses which have, on their own, also been normalized to the untreated control (Dataset A {1.12,2.1,2.05,1.85,1.78}, Dataset B {2.2,3.1,2.5,2.9,1.8}), how do I do this the right way? Thanks a lot!


title: Article recommendation dataset

text: I am experimenting with a recommendation systems for articles (news, magazines etc). Is there any dataset available for the same which includes the demographic information of users as well? I am actually trying to incorporate demographic information as well in recommending, hence looking for such a dataset. I really don't want to simulate such data as I think real data would give much more interesting insights than simulated one. Anyone aware of such? I found one http://contest.plista.com/ but unfortunately its has german dataset and demographic information seems to be missing. There is Yahoo Front Page but its only for academic research. Would appreciate any information here.


title: Hierarchical models for multiple comparisons - multiple outcomes context

text: I've just been (re-)reading Gelman's _Why we (usually) don't have to worry about multiple comparisons_. In particular the section **"Multiple outcomes and other challenges"** mentions using a hierarchical model for situations when there are multiple related measures from the same person/unit at different times/conditions. It appears to have a number of desirable properties. I understand that this is not necessarily a Bayesian thing. Could somebody show me how to properly construct a multivariate multilevel model using rjags and/or lmer (regular JAGS and BUGS should be fine too, as well as other mixed model libraries e.g., MCMCglmm) so that I can play around with it to compare and contrast results? The type of situation I would like a model for is reflected in the toy data below (multivariate, repeated measures):               set.seed(69)     id     <- factor(rep(1:20, 2))                # subject identifier     dv1    <- c(rnorm(20), rnorm(20,  0.8, 0.3))  # dependent variable 1 data for 2 conditions     dv2    <- c(rnorm(20), rnorm(20,  0.3, 0.6))     dv3    <- c(rnorm(20), rnorm(20, -0.3, 0.8))     dv4    <- c(rnorm(20), rnorm(20,  0.2, 1  ))     dv5    <- c(rnorm(20), rnorm(20,  0.5, 4  ))     rmFac  <- factor(rep(c(1, 2), each=20))       # repeated measures factor     dvFac  <- factor(rep(1:5, each=40))           # dependent variable indicator          dfwide <- data.frame(id, dv1, dv2, dv3, dv4, dv5, rmFac)     dflong <- data.frame(id, dv = c(dv1, dv2, dv3, dv4, dv5), rmFac, dvFac) # just in case it's easier?


title: Are those two models equals?

text: > **Possible Duplicate:**   >  R-squared: X "explains" the percentage of variation of the Y values. Does > axis order matter? Excuse me about the strange question but I have a doubt regarding the linear regression. Take a look at the code below:                 x = rnorm(100)       y = rnorm(100)            mod1 = lm(x~y+0)       mod2 = lm(y~x+0)      As You can see I change the regressand and the regressors on those two linear models. My question regarding the residuals of those two models, why they should be different if I use the same variables?


title: Target and Non-target scores, BOSARIS Toolkit

text: I am working on a binary classification problem using Artificial Neural Networks. My questions are: 1) What are "target scores" and "non-target scores" in a classifier, and how can I make my classifier such that I can obtain scores instead of hard decision in its output. (I am using Artificial Neural Networks toolbox of Matlab to perform the classification.) 2) The output of my classifier are "Estimated outputs", "True Positive Ratio", "True Negative Ratio", and "Threshold". I have also "actual targets". How can I relate these matrices to the "target scores" and "non-target scores" to be able to use them in BOSARIS Toolkit to calculate Equal Error Rate of the system?


title: How to check if the results of two methodologies are the same?

text: I have two models that essentially are supposed to produce the same results. Given that the first model is correct, what tests are available that I can use to check the validity of the second model by comparing the results they both produce? Ideally the comparison of the tests will be done through R, but I am also interested in the theory. EDIT about the results: The results of both models are raster images For each pixel in each raster each module produces a values. I can create a sample with n samples eg:               PixelId | modelA | modelB     1       |    223 | 222     2       |     1  |  1      3       |     3  |  3      and so on...


title: Interpreting a mediation model with binary outcome

text: I have a mediation model with two continuous mediators (m1; m2), a continuous input variable (x) and a dichotomous output variable (y). The two mediators are different mechanisms of the input variable. I understand mediation should be useful to determine whether one, both or none of the mediators increases the chances of the outcome variable to be 1. I have tried to calculate the model with two different methods: a) the lavaan package for R, and b) Baron and Kenny's steps. One difference in the calculations is that the lavaan modeling approach calculates the input variable x based on regression coefficients of 4 predictors, whereas for the Baron and Kenny approach the variable x is just the mean of these 4 predictors (I just didn't know how to do it). The outcomes of the two models differ quite a lot as can be seen in the attached figures. However the main problem I have at this moment of time, is the interpretation of the coefficients. In figure 3 and 4 the estimates as they are returned from R are on the respective arrows. ![Lavaan model with estimates](http://i.stack.imgur.com/wMAjT.png) ## Lavaan model with estimates ![Baron and Kenny model with estimates](http://i.stack.imgur.com/htDtg.png) ## Baron and Kenny model with estimates To make some sense out of these estimates, I tried to transform then into a percentage. A percentage increase of the odds ratio of the dependent variable in case of increasing the respective predictor by 1. These are in figures 5 and 6. ![Lavaan model with standardized estimates](http://i.stack.imgur.com/3uZ8x.png) ## Lavaan model with standardized estimates ![Baron and Kenny model with standardized estimates](http://i.stack.imgur.com/DXRTL.png) ## Baron and Kenny model with standardized estimates To do so I did the following for the regression coefficients with binary outcomes:               round((exp(COEFF)-1)*100,2)      Now if I look at the lavaan model, could I say: (1) The total effect of x on y is high (11.52%), the direct effect of x on y is rather low (-4.78%). (2) An increase of 1 in m1, increases the odds of Y to be 1 by 90%. Generally, is it a bad idea to try some kind of transformation here? I can't really understand how to make sense out of the different coefficients in the indirect effects paths.               #### Minimum Working Example ####          library(QuantPsyc)     library(psych)     library(lavaan)          data <- read.csv("http://dl.dropboxusercontent.com/u/63191123/data.csv", header = TRUE, sep ="\t", stringsAsFactors=TRUE)     sapply(data, mode);sapply(data, class);attach(data);head(data);psych::describe(data)               model <- '      # measurement model     m1 =~ m11 + m12 + m13 + m14 + m15 + m16     m2 =~ m21 + m22 + m23 + m24 + m25     m24 ~~ m25     m13 ~~ m14     m16 ~~ m21     m25 ~~ m22     m23 ~~ m24     m13 ~~ m15     m12 ~~ m13     m24 ~~ m22     '          # fit model     fit <- cfa(model, mimic="Mplus", data=data)          m1 <- (parameterEstimates(fit)[1,4] * m11 + parameterEstimates(fit)[2,4] * m12 + parameterEstimates(fit)[3,4] * m13 + parameterEstimates(fit)[4,4] * m14 + parameterEstimates(fit)[5,4] * m15 + parameterEstimates(fit)[6,4] * m16)/6     m2 <- (parameterEstimates(fit)[7,4] * m21 + parameterEstimates(fit)[8,4] * m22 + parameterEstimates(fit)[9,4] * m23 + parameterEstimates(fit)[10,4] * m24 + parameterEstimates(fit)[11,4] * m25)/5               #### Baron and Kenny ####          #<<<<<< Baron and Kenny - Step 1 >>>>>>#     Baron_c.mod <- glm(Y ~ X, family = binomial(link = "logit"))          #<<<<<< Baron and Kenny - Step 2 >>>>>>#     Baron_a1.mod <- lm(m1 ~ X)     Baron_a2.mod <- lm(m2 ~ X)          #<<<<<< Baron and Kenny - Step 3 >>>>>>#     Baron_b.mod <- (glm(Y ~ X + m1 + m2, family = binomial(link = "logit")))          # Baron and Kenny Summary:     Path_a1 <- summary(Baron_a1.mod)$coefficients[2,]          Path_a2 <- summary(Baron_a2.mod)$coefficients[2,]     Path_b1 <- summary(Baron_b.mod)$coefficients[3,]         Path_b2 <- summary(Baron_b.mod)$coefficients[4,]     Path_c <- summary(Baron_c.mod)$coefficients[2,]         Path_c_ <- summary(Baron_b.mod)$coefficients[2,]          Baron <- rbind(Path_a1, Path_a2, Path_b1, Path_b2, Path_c, Path_c_)     print(Baron, digits=3)          #Baron.std <- rbind(lm.beta(Baron_a1.mod), lm.beta(Baron_a2.mod), lm.beta(Baron_b.mod)[2], lm.beta(Baron_b.mod)[3], lm.beta(Baron_c.mod)[1], lm.beta(Baron_b.mod)[1])     #Baron.std <- Baron.std[,1]     #Baron <- data.frame(Baron, Baron.std)               ################     #### Lavaan ####     ################                data$m1 <- m1         data$m2 <- m2          mediation <- '     X =~ x1 + x2 + x3 + x4     Y ~ b1*m1 + b2*m2 + c*X     m1 ~ a1*X     m2 ~ a2*X     indirect_m1 := a1*b1     indirect_m2 := a2*b2     total := c + (a1*b1) + (a2*b2)     '          mediation.fit <- sem(mediation, std.lv=FALSE, ordered = "Y", data=data) #      summary(mediation.fit)               ###### Standardized coefficients for logit regressions:     # e.g. for the logistic regressions:     # Baron and Kenny Model     round((exp(summary(Baron_b.mod)$coefficients[3,1])-1)*100,2) # b1         round((exp(summary(Baron_b.mod)$coefficients[4,1])-1)*100,2) # b2     round((exp(summary(Baron_c.mod)$coefficients[2,1])-1)*100,2) # c         round((exp(summary(Baron_b.mod)$coefficients[2,1])-1)*100,2) # c'          # Lavaan Model     round((exp(parameterEstimates(mediation.fit)[5,5])-1)*100,2) #b1     round((exp(parameterEstimates(mediation.fit)[6,5])-1)*100,2) #b2     round((exp(parameterEstimates(mediation.fit)[28,5])-1)*100,2) #c     round((exp(parameterEstimates(mediation.fit)[7,5])-1)*100,2) #c'


title: Extract BIC and AICc from arima() object

text: **Problem:** I would like to extract the BIC and AICc from an arima() object in R. **Background:** The arima() function produces an output of results, which includes the estimated coefficients, standard errors, AIC, BIC, and AICc. Let's run some sample code to see what this looks like:               # Load the sunspots dataset     data(sunspots)     # Build an ARIMA(2,0,2) model and store as an object     model <- arima(x=sunspots, order=c(2,0,2), method="ML")     # Show a summary of the model     model       The output of results for the model appears like this:               Series: sunspots      ARIMA(2,0,2) with non-zero mean           Coefficients:              ar1     ar2      ma1      ma2  intercept           0.9822  0.0004  -0.3997  -0.1135    51.2652     s.e.  0.1221  0.1196   0.1206   0.0574     8.1441          sigma^2 estimated as 247.9:  log likelihood=-11775.69     AIC=23563.39   AICc=23563.42   BIC=23599.05      On the bottom line, we can see values for AIC, BIC, and AICc. (Note: this is the output shown by arima() when the forecast package has been loaded, i.e. library(forecast)) Accessing the AIC value is quite easy. One can simply type:               > model$aic     [1] 23563.39      Access to the AIC value in this manner is made possible due to the fact that it's stored as one of the model's attributes. The following code and output will make this clear:               > attributes(model)     $names      [1] "coef"      "sigma2"    "var.coef"  "mask"      "loglik"         [6] "aic"       "arma"      "residuals" "call"      "series"        [11] "code"      "n.cond"    "model"              $class     [1] "Arima"      Notice, however, that bic and aicc are not model attributes, so the following code is no use to us:               > model$bic     NULL     > model$aicc     NULL      The BIC and AICc values are, indeed, calculated by the arima() function, but the object that it returns does not give us direct access to their values. This is inconvenient and I've come across others who've raised the issue. Unfortunately, I've not found a solution to the problem. Can anyone out there help? Which method can I use to access the BIC and AICc from the Arima class of object. **Note:** I've suggested an answer below, but would like to hear improvements and suggestions. Edit (Version details as requested):               > R.Version()     $platform     [1] "i686-pc-linux-gnu"          $arch     [1] "i686"          $os     [1] "linux-gnu"          $system     [1] "i686, linux-gnu"          $status     [1] ""          $major     [1] "3"          $minor     [1] "0.2"          $year     [1] "2013"          $month     [1] "09"          $day     [1] "25"          $`svn rev`     [1] "63987"          $language     [1] "R"          $version.string     [1] "R version 3.0.2 (2013-09-25)"          $nickname     [1] "Frisbee Sailing"


title: Few Qs about contingency tables: Which Distribution?

text: I've collected data about the number of visible symptoms of a disease from few hospitals (5 classes) as well as some other data such as sex (2 classes), age class (3 classes). The data were collected on 3 different years (i.e., 3 classes). I am trying to use log linear models on contingency table to analyze the data with the aim of knowing if the prevalence of the number of visible symptoms is different between years, hospitals, sexes and age classes. What I am not sure about is which distribution should I use and how to test for the best distribution for my data? I would appreciate any clues on what to do, thanks


title: paired t-test advantages - literature reference

text: I am having trouble finding any web-based resource or (even better) publication which details (or even names) the kind of problems which one avoids by running a paired t-test as opposed to an unpaired test when testing the same population twice. I am surprised resources like wikipedia also do not detail this. Will an unpaired test overestimate or underestimate my _p_? Or does this depend on other considerations? Could you help me out?


title: What are features that distinguish clustering, blind signal separation and dimensionality reduction?

text: In terms of `input -> [process] -> output` what are features that distinguish clustering, blind signal separation and dimensionality reduction? From this Wikipedia article, the implication is that there are two types of unsupervised learning:   * clustering; and    * blind signal separation I have never heard of the term _blind signal separation_ before. How does it differ from clustering and how does dimensionality reduction underpin this?


title: Show that expected Fisher information for location families is a constant function

text: We have a location family, with location parameter $\theta$, and densities defined over the real numbers $\mathbb{R}$ such that $p_\theta(x) = p(x-\theta)$. I have to prove that the expected Fisher information $J(\theta) = \int_\mathbb{R} \frac{(p'(x))^2}{p(x)}\,\mathrm{d}x$ implying that $J(\theta)$ is a constant function. Is it just about solving the integral? If yes, please a hint, my weak math background does not allow me to solve this...


title: How do I establish shape and scale parameters?

text: I am trying to forecast customer LTV using an exponential gamma distribution suggested in a Journal of Forecasting article (Empirical Comparison of New Product Trial Forecasting Models; authored by Bruce Hardie, Peter Fader, and Michael Wisniewski). The specific formula that I am looking to solve is $P(t)=1-(α/(α+t))^r$, where $t$ is the period, $P(t)$ is the probability of a customer still being a customer at time $t$, $α$ is the scale parameter, and $r$ is the shape parameter. I have some initial data on the per period attrition / retention of customers over 12 periods but I don't know how to use this data to calculate $α$ or $r$ to enable me to forecast future periods attrition/retention and ultimately LTV. **Can anybody explain how I can use the data I have to calculate $α$ and $r$?**


title: SVM low accuracy but significant p value

text: I'm also new at doing SVM... We have a data set of about 50 (25 in each group). We are trying to identify babies at a high risk of developing a disorder from babies that are low risk for developing a disorder. When I run SVM (LOOCV, linear kernal, ttest filtering), I'm getting really low accuracy levels (i.e. 24% in one case), but with statistically significant p-values. Can anyone help me understand why this is happening/what this could mean? Thanks!


title: Papers on the risks of using regular subsets

text: I'm working on my Masters thesis and am using exhaustive search (also called regular subsets) for linear models. There was an answer on CrossValidated, I believe by whuber, which talked about the risks of using regular subsets and had a link back to a page with several papers discussing this risk. Unfortunately, I have been searching for this post and I can't find it. If anyone could point me toward this post, or papers on the limitations/risks/problems with using exhaustive search for linear models, I would be very grateful.


title: Overdispersion in GLMM in SPSS

text: I am using a generalized linear mixed model to analyze poisson and binary data. Animals were observed on several moments, so my model should account for that, which is why I am using the GLMM. I am convinced that the model is adequate, however, I do need to correct for overdispersion. However, I cannot find the way to do this in SPSS 20. I am not an SPSS-expert, so I have not found it. Can anyone tell me if it is possible to address overdispersion in SPSS GLMM? Thank you!


title: Whether and how to calculate adjusted r-squared in factor analysis?

text: I am using the `factanal` function in `R` to conduct factor analysis. Somebody told me that I need to use Adjusted $R^2$ to get more accurate results in factor analysis. But when I looked at the option of factor analysis, I couldn't find the way to do that. Any ideas?


title: Why does a model perform worse after reintroducing observations with missing data imputed?

text: My dataset contains about 12% missing data, and much of the missing data is grouped along observations (not randomly scattered or along columns). I optimized a regression method after removing the observations with large amounts of missing data (and used KNN-impute to infer the remaining missing data). I am now trying knn-impute on the entire data set, without removing any data prior, and my previous optimal regression method is now performing worse. Why might this be? Note: optimal regression method == optimizing parameters. So is adding in this new data worse i) as it contains lots of missing data and knn-impute is the wrong method, or ii) I just need to search for new optimal parameters.


title: What test should i perform?

text: I have a questionnaire with a five point scale. Differences in means is obtained. As my data is not normal I am using NON-PARAMETRIC tests for my comparisons. Population: 1 group (Pharmacists only) Question example:               Statement                                                Min   Max  Mean Std.D      Are important reference tool while  dispensing medication  1    5   3.4  .829     Are important reference tool while responding to symptoms  1    5   3.7  .876      How do you test if a significance exists between these two statements?


title: Calculating error/accuracy in Logistic Regression

text: I am using Logistic Regression to create a binary classifier in R. Wrote a predictor for test data, and calculated the mean of my predicted result (I guess it can be translated to the efficiency of my method). I did find some references to training and test errors calculation, but can't seem to understand how to calculate them from the predicted fit. Any reference will help!                   glm.fit=glm(salary~age+workclass+fnlwgt+education+education.num+marital.status             +occupation+relationship+race+sex+capital.gain+capital.loss             +hours.per.week+native.country,             data=trainAdults,family=binomial)         summary(glm.fit)         glm.probs=predict(glm.fit,newdata=testAdults, type="response")          glm.probs[1:10]         glm.pred=ifelse(glm.probs>0.5,"2","1")         setattr(testAdults$salary,"levels",c(1, 2))         attach(testAdults)         table(glm.pred,salary)         mean(glm.pred==salary)         plot(glm.probs)         boxplot(glm.probs)      Result of summary:               Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.6016  -0.6337  -0.3539  -0.0786   3.3131            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)    -1.433467   0.029795 -48.111  < 2e-16 ***     age             0.437952   0.026454  16.555  < 2e-16 ***     workclass      -0.104634   0.021628  -4.838 1.31e-06 ***     fnlwgt          0.052655   0.023348   2.255 0.024117 *       education       0.055044   0.028739   1.915 0.055455 .       education.num   0.851533   0.026941  31.607  < 2e-16 ***     marital.status -0.342673   0.026837 -12.769  < 2e-16 ***     occupation     -0.027290   0.023192  -1.177 0.239308         relationship   -0.223131   0.033368  -6.687 2.28e-11 ***     race            0.100934   0.026330   3.833 0.000126 ***     sex             0.404090   0.034456  11.728  < 2e-16 ***     capital.gain    2.392861   0.108381  22.078  < 2e-16 ***     capital.loss    0.259485   0.019611  13.232  < 2e-16 ***     hours.per.week  0.335252   0.025005  13.407  < 2e-16 ***     native.country -0.005913   0.024925  -0.237 0.812480         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 16832  on 15080  degrees of freedom     Residual deviance: 11829  on 15066  degrees of freedom     AIC: 11859          Number of Fisher Scoring iterations: 7