title: Good clustering Java library

text: I'm looking for a good Java library implementing several clustering algorithms. I'll have to cluster some programs execution traces and I still don't know which algorithms I am going to need, so I'd like to use a library providing lot of them and that makes it easy to swap algorithms. So far I had a look at Weka but I don't know whether there is a more complete library available I'm missing.


title: Why are the residuals in this model so linearly skewed?

text: I am fitting a simple linear regression model with 4 predictors: `lm(Outcome ~ Predictor1 + Predictor2 + Predictor3 + Predictor4, data=dat.s)` I'm finding that the model predictions are consistently off as shown in this graph: ![scatterplot of predictions and residuals](http://i.stack.imgur.com/CNLJz.png) The model clearly overestimates the low values and underestimates the high values, but the miss-estimation is very linear -- it seems like the model should be able to just adjust the slope and fit the data better. Why is that not happening? In case it helps, here are scatterplots of the the Outcome against each of the four Predictors: ![enter image description here](http://i.stack.imgur.com/uc55e.png) Using the `car` package `outlierTest` function did not identify any outliers.


title: PCA proof needed for proportion of variance explained by L PCs = mean R-square from regression on PC scores

text: I observed the following relation and would like to know where I can find a general proof for this: Assume a data matrix $A = [a_{ij}]_{t x k}$. 1) Perform principal component analysis (PCA) using the correlation matrix $corr(A)$ 2) Calculate the proportion od variance explained by the first $L$ principal components $: \frac{\sum_{l=1}^{L} Eigenvalue_l}{\sum_{l=1}^k Eigenvalue_l}$ 3) Take the first $L$ principal component scores, $PCS_{l}$, perform the following regression and store the respective coefficient of determination $R^2$: $a_1 = b_{10} + b_{11}PCS_{1} + \dots + b_{1L} PCS_{L} + e_{1}$ => $R^2_1$ $\vdots$ $a_k = b_{k0} + b_{k1}PCS_{1} + \dots + b_{kL} PCS_{L} + e_{k}$ => $R^2_k$ 4) Take the average: $\bar{R}^2 = \frac{1}{k}\sum_{i=1}^k R^2_i$ 5) Then, I observe that $\bar{R}^2 = \frac{\sum_{l=1}^{L} Eigenvalue_l}{\sum_{l=1}^{k}Eigenvalue_l}$ The relation in 5) is an interesting result and I would like to find a general proof for this. Does anyone know a good statistics textbook or a paper that shows this result? Thanks H


title: What to look for in a pocket calculator for a Graduate level statistics course midterm / final

text: Soon I'll be sitting for a midterm (and then the dreaded final) for a graduate class in Statistics. It's open book / notes but I'll be needing a calculator -- something I've not used for years. Any thoughts on what features I should look for when shopping for one. I've always wanted to own and use an HP 10 series calculator but I was struck with sticker shock when I looked them up on eBay. Since a tag does not exist for this topic, I'll be tagging it with topics that I'll be examined on. Thanks in advance.


title: Weights in Cumulative Link Models - package Ordinal

text: I am running a cumulative link model on an ordered response variable that represents an assigned preference value between 0 and 4. The setup of my research is a questionnaire in which respondents are asked to rate a collaboration scenario (higher score means a more preferable scenario. I already know that specific scenario elements (e.g. funding) are highly influential in the determination of preference scores. Christensen (2013) suggests in the appendix of his "analysis of ordinal data" that adding `weights` using the frequency of an important variable in the CLM procedure can be sensible (Or at least I infer that from some examples). Hence I calculate weights based on the frequency of different types of funding (4 types in 1,515 observations:               tab2 <- data.frame(with(dt, tapply(dt$dum,funding,sum)))         tab2$z <- rownames(tab2)     colnames(tab2) <- c("fundweight","group")     dt$fundweight <- tab2$fundweight[match(dt$funding,tab2$group)]      The following happens when adding these weights to a previous model (output anova). So the only difference is the presence or absence of weights: AIC goes from 4,510 to 1,683,600, log Likelihood from -2,223 to -841,768. While these seem to clearly indicate that weights are highly problematic, I do not really understand why this happens. Additionally almost all covariates in my model become highly significant when adding the weights which also somehow increases my number of observations from n = 1,515 to n = 574,303. I don't find any information in the package information nor in the CLM tutorial about what these weights actually represent and why they have such a massive impact on the results. Any pointers would be very welcome! Thanks in advance, Simon


title: Machine Learning and training time: is it really relevant?

text: I have a question regarding the time needed for training a classifier. I am facing the specific problem of Sentiment Analysis (classification of text as pos/neg/neu). (Excepting online learning algorithms) Since the training stage of a classifier is usually conducted only once and offline, in the end does it really matter how much time would it take for training ? If the classifier is to be trained over a mobile or a resource-constrained platform, reducing such training time would still make more sense, though... Thanks in advance.


title: Estimation problem in econometric model

text: I am trying to estimate the following equation: $$Y_{i,t+1}=(λβ)X_{i,t}^{'}+(1-λ)Y_{i,t}+\epsilon_{i,t+1}$$ but I don't understand how to find the $\lambda$ estimate especially since $X'$ consists of multiple variables. In this equation, $\lambda$ is a scalar and $\beta$ is a vector. Can anyone point me in the right direction in terms of model type? I am using Stata 12 for this problem and my data is an unbalanced panel. Thanks.


title: How can I estimate a "rank" dependent variable in a multivariate dataset?

text: I have a dataset where $x_{1}$, $x_{2}$, ..., $x_{k}$ are of mixed type and $y$ is a ranking from $1$ to $n$ where $n$ is the number of observations in the dataset. I'd like to predict $y$ based on the $x$s. What is a good way to proceed here? My intuition is that a GLM approach won't be particularly useful given that $y$ is not independently distributed (although the $x$s are iid).


title: Estimating n in coupon collector's problem

text: In a variation on the coupon collector's problem, you don't know the number of coupons and must determine this based on data. I will refer to this as the fortune cookie problem: > Given an unknown number of distinct fortune cookie messages $n$, estimate > $n$ by sampling cookies one at a time and counting how many times each > fortune appears. Also determine the number of samples necessary to get a > desired confidence interval on this estimate. Basically I need an algorithm that samples just enough data to reach a given confidence interval, say $n \pm 5$ with $95\%$ confidence. For simplicity, we can assume that all fortunes appear with equal probability/frequency, but this is not true for a more general problem, and a solution to that is also welcome. This seems similar to the german tank problem, but in this instance, fortune cookies are not labeled sequentially, and thus have no ordering.


title: How to test uniformity in several dimensions?

text: Testing for uniformity is something common, however I wonder what are the methods to do it for a multidimensional cloud of points.


title: Testing count variables with size as a covariate

text: I think linear regression could handle this but I am asking myself if that method has good statistical power. Let's think of circles with different radii and let there be points inside those circles.                         radius points     Circle A      20      8     Circle 1      10      7     Circle 2      25     12     Circle 3       8     10      How can I test the null hypothesis that the number of points in circle A is the same as in the circles 1,2 and 3 given that the area has an linear influence? ## My idea Here is a very basic idea: What do we expect for each Circle given the simplification that the area of the circle is given and not the radius? In total we have an area of $63$ with $37$ points. Thus $0.59$ points per area. So we would expect $11.75$ points for circle A. As I am interested in the aspect that A differs from the other circles which can be seen as an homogene area i will sum them up to an area of $43$ with $25.25$ expected points. This leads to the Chi-Square-Test Statistic: $$ \chi^2 = \frac{(20-11.74603)^2}{11.74603} + \frac{(43-25.25397)^2}{25.25397} \sim \chi_1 $$ After calculation we obtain the p-value: $0.186$ for this example ## Improvement As we have a sample of of point numbers $x=(7,12,10)$ we can calculate the variance $Var(x) = 6.3$ That's why the formula will change to: $$ \chi^2 = \frac{(20-11.74603)^2}{11.74603} + \frac{(43-25.25397)^2}{3 \cdot 6.333} \sim \chi_1 $$ And the p-value will be: $0.0164$ Is this "safe" ?


title: How to code and interpret a regression with 4-level categorical variable and continuous covariates

text: I have a problem in which 4 groups of people are compared in their performance on a certain test. One way anova shows that they are in fact different. But I also know that my groups are different in some other covariate (e.g. age). I'd like to see that if those differences are still significant if I "adjust" for age etc. When I run a regression it shows me that one of the levels is "significant" but not 2 others (1 is a reference, I guess). Intercept is significant also, but it always is anyway. If I run the model withot the constant term, and with 3 dummy variables then all 3 are "significant". How do I interpret these results? Thanks.


title: How to predict success rates using historical data over time

text: Let's say I have a data set which shows me basketball players: list of their games which contains date, number of shots they took and number of shots they scored. For example: 1.10.1992 Michael Jordan took 12 shots - scored 6 I would like to predict what is going to be their scoring success rates in the next games using the historical data set. What is the best way to do it? The real model I'm doing it for has very low success rates < 1%.


title: Dimensionality reduction on binary matrix

text: Is there a way to do dimensionality reduction on a **binary matrix**. From what I know, Principal Components Analysis does not work. Any suggestions on this matter. Edit: It's not a symmetric matrix.


title: Understanding proxy variable

text: Can somebody help me to understand a topic? Say in my analysis I have the true variable _education_ ($E$) which is measured by a proxy variable _number of years in school_ ($S$). Here, the error is measurement is $E-S$. I was told that, this error is independent of $S$ not of $E$. I can not understand how it can be the case? Should not be that error dependent on both $S$ and $E$? Thanks,


title: How to verify that simulated data is normally distributed?

text: I have a program generating purportedly normal distributions and I would like to test it. I have a number of issues; perhaps the experts here will help me sort out the essential from the inessential and answer most of these.   1. I'm looking for a simple test, ideally -- one I can implement without too much trouble.   2. There may be correlation between adjacent values. Some tests may not be sensitive to this failure if the data are 'otherwise' normally distributed.   3. Ideally, I'd like to allow a (small!) amount of non-normality. Most tests I've seen allow slightly non-normal data to pass only because a small number of values are tested (where "small" might still mean millions, depending on the size of the deviation). This is reminiscent of this question about the value of normality testing.


title: Multivariate regression for spatial dataset

text: I am trying to decide on the correct technique for a multivariate regression with spatial data. I would like to run a regression where the dependent variable is the current snow depth and the independent variables include physiographic parameters (slope, aspect, elevation, etc.) and snow depth for the same site in past years using daily data. The goal is to produce a statistical model with which I can interpolate snow depth across a whole basin based on the physiographic parameters. Initially I was going to use a standard MVR but came across geographically weighted regression (GWR), which I think is more appropriate since snow depth is very spatially correlated. The third step, after establishing a model and interpolating, would be to distribute the residuals that I'll have at points where I know the snow depth; a common approach in the literature is elevation detrended inverse distance weighting.   1. Would it be incorrect to use MVR instead of GWR?   2. If I use GWR, would it still makes sense to distribute the residuals? From what I read, GWR already includes some correction for the inevitable residuals unlike MVR.  Please correct me if I'm wrong or seem to have misunderstand anything. I'm quite new to spatial statistics. Most of my GWR knowledge comes from Geographically Weighted Regression.


title: Calculate p value for Pearson's R

text: I have a bunch of bivariate correlation tests for which $R$ has already been calculated. Is it possible to figure out the p-value for each $R$ knowing only (in addition) the number of data points $n$? If not, then what else do I need to know and how do I calculate it? (I also have results for bivariate regressions, including F, sigma, coefficient, coefficient sigma, intercept, intercept sigma if they are of use).


title: Power of FDR vs. FWER approaches in multiple comparisons

text: Regarding multiple comparisons, could someone please explain me why the power of false discovery rate (FDR) is greater than the power of family-wise error rate (FWER)?


title: How to get ellipse region from bivariate normal distributed data?

text: I have data which looks like: ![Figure](http://i.stack.imgur.com/lB6Yh.png) I tried to apply normal distribution (kernel density estimation works better, but I don't need such great precision) on it and it works quite well. Density plot makes a ellipse. I need to get that ellipse function to decide if a point lies within the ellipse's region or not. How to do that? R or Mathematica code are welcomed.


title: Expected value of a multinomial distribution

text: A multinomial distribution can be given as $ M(m_1,\dots,m_K|N,P) = {N \choose m_1\dots m_K}\prod_k p_k^{m_k} $ The expected value is $Np_k$. How can I prove it?


title: Binomial distribution confidence interval for log plot

text: In simulating iterative decoding of low-density parity-check codes there may be (for a certain signal-to-noise-ratio of a noisy channel) for example 10 decoding failures out of $10^6$ trials. The log of probability of decoding failure is plotted against signal-to-noise-ratio. I quote from the Appendix of MacKay's paper Good Error-Correcting Codes based on Very Sparse Matrices to show how he deals with creating error bars on such plots: > The experiments result in a certain number of block decoding failures $r$ > out of a number of trials $n$. We report the maximum likelihood estimate of > the block error probability, $\hat{p}=r/n$, and a confidence interval > $[p_-,p_+]$, defined thusly: if $r \geq 1$ then $p_{\pm}=\hat{p}\exp(\pm > 2\sigma_{\log p})$ where $\sigma_{\log p}=\sqrt{(n-r)/(rn)}$; else if $r=0$ > then $p_+=1-\exp(-2/n)$ and $p_-=0$. I've read the wikipedia page on Binomial proportion confidence interval and a number of helpful answers on this website. Mackay's formula can be derived by approximating the normal approximation to a 2 sigma confidence interval: $\hat{p}\pm z_{1-\alpha/2}\sqrt{\hat{p}(1-\hat{p})/n} = r/n \pm 2\sqrt{r(n-r)/(n^3)}=(r/n)(1\pm2\sqrt{(n-r)/(rn)})$ which is approximately equal to $(r/n)\exp(\pm2\sqrt{(n-r)/(rn)})$ when the argument of the exponential is small. His notation also suggests he is really doing statistics with $\log p$. His formulas give symmetric error bars on log plots. I don't know whether this is a reasonable goal. So I have a couple of questions:   1. Is my derivation the best way to understand what's going on with this formula?   2. Is it the most appropriate confidence interval?


title: Why use Lasso estimates over OLS estimates on variable subset?

text: For Lasso regression $L(\beta)=(X\beta-y)'(X\beta-y)+\lambda*norm(\beta,1)$, suppose the best solution (minimum testing error for example) selects $k$ features, so that $\hat{\beta}^{lasso}=\left(\hat{\beta}_1^{lasso},\hat{\beta}_2^{lasso},...,\hat{\beta}_k^{lasso},0,...0\right)$. We know $\left(\hat{\beta}_1^{lasso},\hat{\beta}_2^{lasso},...,\hat{\beta}_k^{lasso}\right)$ is biased estimate of $\left(\beta_1,\beta_2,...,\beta_k\right)$, why we still take $\hat{\beta}^{lasso}$ as our final solution, but not the more 'reasonable' one $\hat{\beta}^{new}=\left(\hat{\beta}_{1:k}^{new},0,...,0\right)$, where $\hat{\beta}_{1:k}^{new}$ is the LS estimate from partial model $L^{new}(\beta_{1:k})=(X_{1:k}*\beta-y)'(X_{1:k}*\beta-y)$. ($X_{1:k}$: keep the columns of $X$ corresponding to the $k$ selected features). In brief, why we use Lasso for both feature election and estimation, instead of only for variable selection and leaving the estimation by other models on the selected features? More over, what does 'Lasso can select at most $n$ features'? $n$ is the sample size.


title: Increased Type I error - GLM

text: Some of you might have read this nice paper: O’Hara RB, Kotze DJ (2010) Do not log-transform count data. Methods in Ecology and Evolution 1:118–122. klick. In my field of research (ecotoxicology) we're dealing with poorly replicated experiments and GLMs are not widely used. So I did a similar simulation as O’Hara & Kotze (2010), but mimicked ecotoxicological data. **Power simulations** : I simulated data from a factorial design with one control group ($\mu_c$) and 5 treatment groups ($\mu_{1-5}$). Abundance in treatment 1 was identical to the control ($\mu_1 = \mu_c$), abundances in treatments 2-5 was half of the abundance in the control ($\mu_{2-5} = 0.5 \mu_c$). For the simulations I varied the sample size (3,6,9,12) and the abundance in the control group (2, 4 ,8, ... , 1024). Abundances were drawn from a negative binomial distributions with fixed dispersion parameter ($\theta = 3.91$). 100 datasets were generated and analysed using a negative binomial GLM and a gaussian GLM + log- transformed data. The results are as expected: the GLM has greater power, especially when not many animals were sampled. ![enter image description here](http://i.stack.imgur.com/3fkVr.jpg) Code is here. **Type I Error** : Next I looked at type one error. Simulations were done as above, however all groups had same abundance ($\mu_c = \mu_{1-5}$). However, the results are not as expected: ![enter image description here](http://i.stack.imgur.com/bUJ2b.jpg) Negative binomial GLM showed a greater Type-I error compared to LM + transformation. As expected the difference vanished with increasing sample size. Code is here. **Question:** **Why is there an increased Type-I Error compared to lm+transformation?** If we have poor data (small sample size, low abundance (many zeros)), should we then use lm+transformation? Small sample sizes (2-4 per treatment) are typical for such experiments and cannot be increased easily. Although, the neg. bin. GLM can be justified as being appropriate for this data, lm + transformation may prevent us from type 1 errors.


title: Distribution of complex results for pool of ten-sided dice, specific to custom tabletop RPG system

text: I'm no statistician, but I do find the area of mathematics fun and fascinating. I have a friend who has designed a "roll-mid" dice system for his custom tabletop RPG and as part of the game balancing, I'm interesting in finding out certain result possibilities. An individual die has the following results:               1        : critical failure, total of -1 "success"     2-3, 8-9 : no "successes"     4-7      : +1 "success"     0        : critical success, total of +2 "successes"      A player will roll from one to ten individual ten-sided dice (d10s). I haven't any idea how to approach an issue like this with the additional layer of complexity. The concept, however, is just really cool to me. My question: What is the distribution of "success" totals for any number of dice 1-10?


title: How to capture competitive spatial interactions between multiple stores and customers

text: I am estimating sales with data on customer and store locations and attributes using a Huff-style model, where sales decay with drive time and increase with attractiveness of the store. One hypothetical instance of the problem is illustrated below. The 3 pushpins represent stores and the white flag represents the site where the customers live.** The size of the pins represents the attractiveness of the store (like floor space). The problem I am having is that I would expect red store's sales at the site to be lower than the the green store's even though they are just as far away and have identical attractiveness, because the purple store is somewhat in between the the red store and the site. I would like to translate this intuition into rigorous (but tractable) math so that I add it into my statistical model. I am also having trouble figuring out what this is called in the literature (other than the n-body gravity problem in physics). ![Example](http://i.stack.imgur.com/xCZTH.jpg) ** For the daltonic folks, the red store is in the upper left corner. The green store is in the bottom right. Purple store is to the left of the site.


title: How to generate a non-integer amount of consecutive Bernoulli successes?

text: Given:   1. A coin with **unknown** bias $p$ (Head).   2. A strictly positive **real** $a > 0$. Problem: Generate a random Bernoulli variate with bias $p^{a}$. Does anyone know how to do this? For instance, when $a$ is a positive integer, then one can flip the coin $a$ times and see whether all the outcomes were Heads: if they are then issue '0', otherwise issue '1'. The difficulty lies in the fact that $a$ is not necessarily an integer. Also, if I knew the bias $p$, I could just build another coin with the desired bias. Thanks in advance!


title: Correlation or Regression Analysis

text: I am proposing a study for a dissertation paper. The paper has 1 dependent variable, and 4 independent variables. The variables are all measured over a 4 year period. The data is quantitative, ratio data. What statistical analysis can be run that uses such a small sample. I have received several different responses from my mentors, but not the same, and with different arguments for each. Any advise here would be much appreciated.


title: Minimize the standard deviation of the total values of groups of items (optimization?) in r

text: I am a total novice, but I think what I'm looking for is help with an optimization problem. I have a set of about 26 items, each with a certain value. I want to distribute these items into 8 groups as equally as possible, having granted that a perfectly equal distribution is not possible because of the individual item values. I think, one way of measuring how equally the items have been divided into groups is by taking the standard deviation of the sums of the values of the items in each group. Is there a package in R that could help me to minimize the standard deviation of these sums by altering the allotment of items to groups? FWIW, this is an 'estate division' problem, in which each item has a known value, and no bidding/auctioning can take place.


title: Difference Between Matlabs classify and ClassificationDiscriminant functions

text: I recently came upon a discrepancy between Matlab's classify and ClassificationDiscriminant functions. By default, they should both be doing Linear DiscriminantAnalysis, but when given unbalanced classes, they give very different results. I have created a small test script here at Gist which creates some data that is unbalanced but easily separable. I test my classifiers on the training data for optimal results. For a linear classifier, the result should be unique, but it appears this is not the case.


title: Average and standard deviation of timestamps (time wraps around at midnight)

text: I have lots of sensor data with timestamps like "2014-09-09 16:10:45" and accompanying sensor readings. To get some insight into these I want to find "unusual" events by looking at the average and standard deviation of the time part of the timestamp. How can I handle the wrap-around of time on midnight? A made up example: Imagine power readings being influenced by people turning on machines in the morning (a power sensor would notice increasing values) and turning them off in the evening (decreasing values). I want to find sensor readings that are unusual. This would be decreasing sensor readings in time periods where readings usually increase and increasing readings in time periods where readings usually decrease. My idea was to extract the time part of the timestamps (eg. 12:55:10), convert that to seconds (the day has 86400) and then, split by the tendency of readings (eg only looking at increasing readings) calculate the average and standard deviation. If I then take the time window from "average second of the day minus standard deviation" to "average second of the day plus standard deviation" (maybe using twice the standard deviation) I would have the typical periods and every increasing reading outside this time window would be unusual. **The problem** : Time wraps around at midnight! A reading at 00:15:00 would actually be really close to 23:50:00 in reality, but "far away" in calculation. This surely skews the statistics unless everything happens mid- day. Is there a standard practice to handle this? Can you give me ideas? I am totally stumped at the moment. I would love to stay in PostgresQL but as that is not a requirement I did not tag that. Anything helps! Below is some example data, I have about 200-300 readings per sensor. You can see that in this example the increases happen in the morning.               "Timestamp as %Y-%m-%d %H:%M:%S";"Day of the year";"Second of the day";"Tendency of reading"     "2014-03-01 14:45:00";60;53100;-0.030     "2014-03-03 08:18:00";62;29880;0.150     "2014-03-03 14:17:00";62;51420;-0.120     "2014-03-03 16:37:00";62;59820;-0.030     "2014-03-04 08:11:00";63;29460;0.150     "2014-03-04 10:21:00";63;37260;-0.150     "2014-03-04 16:12:00";63;58320;-0.030     "2014-03-05 08:04:00";64;29040;0.150     "2014-03-05 14:42:00";64;52920;-0.060     "2014-03-05 17:27:00";64;62820;-0.030     "2014-03-06 08:29:00";65;30540;0.090     "2014-03-06 12:06:00";65;43560;-0.030     "2014-03-06 13:49:00";65;49740;-0.120     "2014-03-07 08:21:00";66;30060;0.150     "2014-03-07 10:27:00";66;37620;-0.030     "2014-03-07 11:27:00";66;41220;0.030     "2014-03-07 13:46:00";66;49560;-0.060     "2014-03-07 16:59:00";66;61140;-0.030     "2014-03-07 18:52:00";66;67920;-0.030     "2014-03-08 08:47:00";67;31620;0.120


title: probability of a pair in cards

text: I have a simple problem involving probability of drawing at least 1 pair of cards in a four card hand. I am not getting the right answer but I dont understand the flaw in my logic. Can anyone explain to me why my approach is wrong? The problem: Bill has a small deck of 12 playing cards made up of only 2 suits of 6 cards each. Each of the 6 cards within a suit has a different value from 1 to 6; thus, for each value from 1 to 6, there are two cards in the deck with that value. Bill likes to play a game in which he shuffles the deck, turns over 4 cards, and looks for pairs of cards that have the same value. What is the chance that Bill finds at least one pair of cards that have the same value? My solution: Probability of 1 or more pair = ((6) * (10 choose 2))/(12 choose 4). i.e. the # of hands with at least 1 pair are 6 (the number of ways to create a pair) * (10 choose 2) (the # of ways to select 2 cards from the remaining 10 cards after the pair). This simplifies to 6/11, however, the correct answer is 17/33. Any help understanding would be greatly appreciated.


title: Expected value of a product of random variables

text: Knowing the distribution, \begin{array}{c c c c} & \quad &X=1\quad &X=-1\quad &X=2 \\\ &Y=1\quad &0.25\quad &0.25\quad &0 \\\ &Y=2\quad &0\quad &0\quad &0.5 \end{array} How does one calculate E(XY)?


title: Making sense of scatter plots

text: I had a recent question which I probably should re-formulate to a more general one. I came across this: Using scatter plots to understand multiple values of Y for a given X and thought the accepted answer was very good, but what's unclear to me is: given a nasty scatter plot, how would you visually get the idea what sort of relationship the x:es and y:s have? In my mind, there's no way of telling if a linear, quadratic, etc regression is appropriate. If we look at the plot that is currently confounding me: ![scatterplot](http://i.stack.imgur.com/2S3uQ.png) My idea was to somehow plot the AVERAGE of f(x), with x clumped together in intervals, instead of each observation. Is this how you would go about it? If not, what other way do you visually make sense of this data?


title: Are world cup predictions testable?

text: As of today, dozens of soccer world cup predictions exist, some more complex, some more elegant, and most of them predict every nation's "chance" of winning a particular match/ the cup. As I am writing a blog post for a newspaper and want to make the statement that these predictions are nice but essentially worthless as they can not be tested, I would like to bulletproof my assumptions and enhance my limited statistical knowledge with this question. **First: Chance**. Weirdly, for instance, although Nate Silver explains his model pretty well, he does never explain the actual notion of "chance". If 538's model "predicts" that Brazil has a 45% chance of winning the tournament, this does indeed mean that if the exact same cup was played a thousand times, and the model _would be correct_ , Brazil would win 450 times on average. And _not win_ in 550 cups, on average. Am I right? **Second: Evaluation of the predictive power of a model** : My even more limited knowledge of data mining, machine learning and predictive modelling tells me that usually, in order to assess the strength of a model, a procedure such as x-validation or validation with a dedicated test set is used. For this to be possible, a _test set must exist_. In the case of the world cup, this would actually mean that a model such as 538's would have to be applied to either the past 20+ wold cups (not possible because of lacking historical data and extremely tedious) or to countless "instances" of the current world cup (not possible because of.. well yeah). Just saying that a model is a good predictor because it correctly predicted (most of the matches of) _this_ world cup 2014 is not strictly valid, is it? In the same vein, when Nate Silver states that SPI worked pretty well for 2010, this actually tells us nothing as it could have been just as well _chance_ , in the sense of a random outcome? So, we have _no method_ of telling whether one of these so prominent world cup predictions are actually good models, given we are not applying these model to the next 100 world cups and evaluating them afterwards? What do you mean?


title: Expected value calculation

text: I have great difficulty in calculating $\mathrm{E}\left(\sqrt{T_{1} T_{2}}\right)$ where $T_{i}$ is a random variable distributed according to the Birnbaum-Saunders distribution. Are there any suggestions on how to calculate the expression? The approximation is in Kundu et. al 2010 but when programming in R the approximation diverges, I'm thinking that the error is in the approach and maybe can meet other approach.


title: How to calculate a mean and standard deviation for a lognormal distribution using 2 percentiles

text: I am trying to calculate a mean and standard deviation from 2 percentiles for a lognormal distribution. I was successful in performing the calculation for a normal distribution using `X = mean + sd * Z` and solving for mean and sd. I think I am missing an equation when I try to do the same thing for a lognormal distribution. I looked at wikipedia and trying to use `ln(X) = mean + sd * Z` but I'm getting confused whether the mean and sd in this case are for the normal distribution or the lognormal. Which equations should I be using? and will I need more than 2 percentile to solve the calculations?


title: (Multiple) Correspondence Analysis for count data entered as binary variables

text: I have a data set, 1014 cases and 55 variables which are binary and is in the form of               Brand1Attribute1, Brand1Attribute2, ... , Brand1Attribute11, Brand2Attribute1, ... Brand5Attribute11      with 0 = "No, I do not associate this attribute with this brand" and 1 = "Yes, I associate this attribute with this brand". Basically, I have been trying to use multiple correspondence analysis in SPSS to plot a perceptual map of how the brands and attributes relate, but I have been getting nowhere. I am imagining a 2-dimension plot, with no meaningful axes, and the distance between points is proportional to the correlation between brands, or attributes, or between brand and attribute. Am I going about this the right way, or should I be using a different method?


title: How can I use Akaike's Information Criterion to compare two models of multi-peaked emission spectra?

text: I have several photoluminescence emission spectra that I am trying to fit curves to. The spectra each have a slight baseline and four peaks. The independent variable $x_i$ is wavelength (converted to energy for analysis) and the response variable $y_i$ is intensity of light recorded at that wavelength. The experimental errors are assumed to be $\sigma_i = \sqrt{y_i}$ . I'm fitting these data with two different models. Model 1 has four Gaussian peaks and a baseline of "$y = mx + c$ type", i.e. 14 parameters in total. Model 2 is again four peaks + linear baseline, but this time the each peak is fitted with a pseudo-Voigt function whose parameters are area, "mixing", mean, and FWHM, meaning a total of 18 parameters for the second function. The fits are quite similar but there are slight differences, with Model 2 seeming to be a slightly better fit to every spectrum in the dataset. I would like to use the AICc to compare the models, because I wish to ascertain whether the Model 2 fit is genuinely a better representation of the physical processes that generate the spectra, or if it is just overparameterized. But I don't know the likelihood functions or how to calculate or maximize them. I looked here and thought that perhaps I could use the $\chi^2$ values from each fit to calculate the AIC in the absence of a maximum likelihood estimate. But I am not sure that the assumption "that the underlying errors are normally distributed (with mean zero) and independent" is met here. If "underlying errors" refers to the residuals after each fit, they do seem to have normal distributions (but not with mean zero); if it actually refers to the experimental errors, then there is not a normal distribution (see above). I don't really understand what "likelihood function" means in the context of photoluminescence emission spectra. After all, the spectra are not records of events following a binomial probability distribution like the coin toss, or a Poissonian probability distribution like the electronic shot noise, that seem to be the bases of all the examples I could find. I am at a loss how to understand the probability mass function in this case. Is it justifiable to use the model functions as the probability mass functions and derive the likelihood functions from them, or is that just circular logic? I'd be grateful for any advice.


title: Machine Learning System Design: a practical advice

text: Recently, I stared working on a machine learning competition hosted on Kagge.com. As the first step, a quick and dirty system was developed using Logistic Regression (LR). After running the system with increasing the number of training examples step by step and plotting the learning curve, I realized that my system is suffering from the high-biased issue. In order to overcome this issue, I increased the number of features step by step and measure the training and cross validation errors. Unfortunately, this did not show a significant improvement and both training and cross validation accuracies remained in the 76% - 79% range. At this time I’m considering to follow one of the following avenues.   1. Try few more learning algorithms (single layer NN, SVM or decision tree) and take the majority vote.   2. Since this is a text (actually web pages) analysis problem, instead of unigram, try bigram (and trigram) with LR and check whether I can achieve some significant improvement. Your expertise advices are highly appreciated.


title: German credit data: neural network, svm, logistic regression : input variables

text: I'm using the following data set on some credit scoring models: https://archive.ics.uci.edu/ml/datasets/Statlog+(German+Credit+Data) My teacher told me that it's best to use the same data set for all the different techniques, but how do you handle the different restrictions? The data set consists of 7 numeric and 13 categorical variables, how do you use those categorical variables for the sum? Doesn't the support vector machine only accept 1 or 0 as input? Or is it values ranging between 0 and 1 ?


title: Discrete Choice Model with Some Levels Combined

text: I am new to discrete choice modeling, and I am half way through reading Train's Discrete Choice Methods with Simulation. My question is: what do I do if some of the choice levels are combined? For instance, if there are choices 1,2,3,4,5, but I only know whether an individual selected 1, or 2, or one of 3,4,5, but not exactly which one of 3,4,5. I can see how the likelihood function can be constructed in this case -- it would be a marginal likelihood if I am not mistaken. I would like to know:   1. What is the name of this sort of model. I have not been able to find anything on google on my own.    2. What are some practical implemnations (i.e. software packages) that are already out there? Thank you very much!


title: Correct use of Chi-Square?

text: I am trying to test whether one group is significantly different from another in two different collections of data. I have seen that people tend to use chi- squared or fisher's exact for this calculation, but am I doing it correctly? Here is the data:                                               TOTAL   NUMBER   PROPORTION          Group 1:       WGD + SSD         19590     14184  72%                    WGD               19590     7045   36%                    SSD               19590     7139   36%                    Singleton         19590     5406   28%          Group 2:       WGD + SSD         110       102    93%                    WGD               110       75     68%                    SSD               110       27     25%                    Singleton         110       8      7%      If I want to, for example, check whether the number of WGD + SSD is significantly different in group 2 verses group 1, would I set up my contingency table like so?               19590    110     14184    102      This returns a P value of 0.0839 using a two tailed Chi-square with Yates' correction, but the proportions are 93% versus 72%. Have I done something wrong here, or are the numbers in the second group too low compared to the first? Thanks in advance!


title: Between subjects repeated measures ANOVA in R - 5 brands, 11 measurements

text: In my research I have performed a series of measurements on 5 different brands of blocks. Each block has been inspected for deformation under incremental forces (20, 30, 40, 50, 60, 70, 80, 90, 100, 110 and 120 N). The deformation for each force was measured 3 times and the mean values were assigned to each brand for a specific amount of force. I was successful in creating linear regression graphs for these 5 different brands. Now my wish is to see whether a brand makes a significant difference in deformation values and to perform a post-hoc analysis to compare brands among themselves. In other words to compare the linear regression lines. Sorry if what I am saying makes no sense. So far, I have tried the following commands:               anova(lm(Deformation~Force*Brand, data=Data))     lm(Deformation~Force, data=Data))          # and     aov.data = aov(Deformation~Force*Brand, Data)      I have gotten suspiciously low p-values ( _*_ ) which clearly indicates that I might be doing something wrong. I would be grateful if you could help me with this issue.               Force   Brand   Deformation       20  Brand1  0.65       30  Brand1  1.23       40  Brand1  1.25       50  Brand1  2.39       60  Brand1  2.45       70  Brand1  2.93       80  Brand1  3.13       90  Brand1  3.57       100 Brand1  4.68       110 Brand1  4.84       120 Brand1  5.33       20  Brand2  1.24       30  Brand2  1.11       40  Brand2  1.6       50  Brand2  2.13       60  Brand2  2.69       70  Brand2  3.60       80  Brand2  3.90       90  Brand2  3.99       100 Brand2  4.51       110 Brand2  4.74       120 Brand2  5.98       20  Brand3  1.21       30  Brand3  1.37       40  Brand3  2.56       50  Brand3  2.49       60  Brand3  3.17       70  Brand3  3.33       80  Brand3  3.38       90  Brand3  4.2       100 Brand3  4.22       110 Brand3  5.22       120 Brand3  6.28       20  Brand4  0.92       30  Brand4  0.89       40  Brand4  1.2       50  Brand4  1.67       60  Brand4  1.98       70  Brand4  2.25       80  Brand4  3.8       90  Brand4  4.17       100 Brand4  4.94       110 Brand4  5.4       120 Brand4  5.76       20  Brand5  0.69       30  Brand5  1.26       40  Brand5  1.61       50  Brand5  2.17       60  Brand5  2.07       70  Brand5  3.35       80  Brand5  3.27       90  Brand5  4.13       100 Brand5  4.25       110 Brand5  4.59       120 Brand5  5


title: Global warming testing

text: I'm very interested in weather. I use data from weather-stations. Here are 2 examples of changing _average annual temperature_ in Vienna and Moscow.   (I've got more, but only 2 links allowed ((   _LinearModelFit_ method from Mathematica have been applied (blue lines).   YearTemp1   YearTemp2   Well we can see, that in model $y=a\cdot x+b$, $a$ is small, but always positive number.   Is it so called global warming, I wonder? There is a great scatter of data, from year to year. So the question to community is:   _are there additional criteria_ (exept just the fact that $a>0$) for someone can say:   "Yes, with probability 50% (20?...90?) we have global warming".   Mathematica has a lot of stats methods so the problem is just "what of them to use". UPD There is similar question,   What is the best test for validating trends in time-series data? but I do not see any real answers.


title: Help me figure out the profile likehihood given this covariance function

text: I'm taking a spatial stats class, and I'm on the road so I can't ask the prof for help. Would appreciate help understanding what is going on here. The problem is set up with $Y = X_s'\beta + e$ (with $s$ indexing coordinates). $e$ is assumed to be second-order stationary, with the variogram $$ \gamma(h) = \left\\{\begin{array}{l r} 0 & h=0\\\ \tau^2 + \sigma^2\left(\frac{3h}{2\rho}-\frac{h^3}{2\rho^3}\right) & 0 < h \leq \rho\\\ \tau^2 + \sigma^2 & h > \rho\\\ \end{array}\right. $$ where $h$ is distance between two points. For second-order stationary processes, apparently it is true that $C(h) = lim_{h' \rightarrow \infty} \gamma(h') - \gamma(h)$, where $C(h)$ is a covariance function, which is only a function of distance, because the processes is second-order stationary. This means that the covariance function is $$ C(h) = \left\\{\begin{array}{l r} \tau^2 + \sigma^2 & h=0\\\ \sigma^2\left(1-\left(\frac{3h}{2\rho}-\frac{h^3}{2\rho^3}\right)\right) & 0 < h \leq \rho\\\ 0 & h > \rho\\\ \end{array}\right. $$ which makes intuitive sense. Now for the trickier part (for me anyway): we are told to assume that $e$ is a gaussian process. We are asked if we might calculate a profile likelihood as a function only of $\rho$. The likelihood function is that of a standard multivariate normal. My question: how do I go about the algebra to (either) find a profile likelihood that is only a function of $\rho$, or perhaps a function of $\rho$ and $\tau$? If the latter, how do I deal with the fact that the parameter $\tau$ will show up on the diagonals of the covariance matrix, but not on the off diagonals? It seems like it makes the algebra intractable. Hints appreciated.


title: Clustering of set of matrices

text: I have 25 matrices 19x19 containing coherence measure between EEG electrodes. I want to divided them into some groups by clustering or any other method. I know how to deal with vectors, but I can't find anything about clustering of set of matrices. If it can help - I think we can use coherence as a distance between cells into matrix. So may be there are some method for clustering of distance matrices? * * * **EDIT** Traditional clustering methods cluster vectors. In the vector space, the distance metric and other distance functions are well defined. The Euclidean distance between vectors $x_1$ and $x_2$ is $|x_1-x_2|_2$, the 2-norm of $x_1-x_2$. Analogously, to compare two matrices $M_1$ and $M_2$, we may want to compute $|| M_1-M_2||_p$, the p-norm of $M_1-M_2$. * * * **EDIT** While there are two solution:   1. Compute distance between matrices as sum of squares of the difference matrix   2. Unwrap a triangle of each matrix into a vector **Is there something else?**


title: Mutual Information really invariant to invertible transformations?

text: "Estimating mutual information" [A Kraskov, H Stögbauer, P Grassberger - Physical Review E, 2004] states that > Mutual information is invariant under reparametrization of the marginal > variables. If $X ′ = F(X)$ and $Y′ =G(Y)$ are homeomorphisms [ie. smooth > uniquely invertible maps], then $$I(X,Y ) = I(X′,Y′)$$ This paper is also used on Wikipedia to justify the same claim. But if that is true, doesn't that imply that if $F=G$ and $X=Y$ we get $$H(X) = I(X,X) = I(X',X') = I(F(X),F(X)) = H(F(X))$$ where $H$ is the differential entropy? Wouldn't this be "proof" that differential entropy is also invariant wrt. such transformations (which it obviously isn't because e.g. for a constant $a$, $H(aX) = H(X) + \log|a| \neq H(X)$)? In particular I'm wondering if $I(X,F(X))=H(X)=H(F(X))$ for all homeomorphisms $F$. Can someone help me reduce my _uncertainty_? PS: I'm talking about the continous case, ie. differential entropy and differential mutual information


title: audit sample prediction

text: I have a finite population of transaction data. The data is broken down into various sub-categories. A subset of these transactions has been audited and determined to either pass or fail. I want to undertake a logit (and will try probit) regression to make predictions about the un-audited transactions. How do I compute the aggregated category-level and overall population prediction interval? And, separately, how should I know if the process and the number of samples selected was appropriate?


title: Interpreting $R^2$, F-statistic & p-value of a model

text: From p. 277 of _R Cookbook_ : > Let's say I have a R model `lm(formula = y ~ u + v + w)` and the `Summary()` > shows: >   >   >     Multiple R-Squared: 0.4981, Adjust R-Squared: 0.4402 F-statistic: >     8.603 on 3 and 26 DF, p-value: 0.0003915 >   > > Using Adjusted `r-Squared` I can say that my model explains 44.02% of the > variance of `y` with the remaining 55.98 unexplained. Question: **Does the associated F-statistic (with the p-value being < .05) tell me:**   1. the model, in general, is significant (not taking into account other values from Summary)   2. the model is significant in explaining the 44.02% variance (adjusted r-squared)


title: Interpreting significance in hierarchical regression

text: I've run a Hierarchical Multiple Regression to understand the unique predictive value of three different maternal variables (maternal PTSD, maternal depression, and quality of the parent-child relationship) in child behavior following trauma exposure (measured by two different subscales of Externalizing behavior problems and Internalizing behavior problems). Based on theory, I decided to enter the variables in the following order and ran the regression twice (once with externalizing behavior as the DV and once with the internalizing behavior as the DV):   1. Child's age & sex   2. Number of traumatic events child experienced    3. Maternal PTSD symptoms score   4. Maternal Depressive symptoms score    5. Quality of parent-child relationship I also included some interaction terms in a sixth step but they were not significant so I dropped them from the model. In the final model, both maternal depression and the quality of the parent- child relationship had significant standardized beta coefficients, and change in R squared was also significant for both the externalizing and internalizing subscales. This would allow me to conclude that the quality of the parent- child relationship uniquely predicts child behavior following trauma exposure, controlling for child age, sex, trauma history, maternal PTSD and maternal depressive symptoms. But how do I interpret the fact that Maternal Depressive symptoms remain significant after the addition of the last step variable (parent-child relationship)? (How would I explain that in plain English?)


title: Contrasts for one factor in a repeated-measures ANOVA

text: I am trying to set custom orthogonal contrasts for a repeated-measures design using R. Is it possible to set contrasts just for one factor? The model looks like this:               lme(percent.diff.mean ~ factor1*factor2*factor3, random = ~ 1|id/factor1/factor2/factor3, data=data, method="ML")      Factor2 and factor3 both have 2 levels. Factor 1 has 4 levels, and it makes sense to divide this into two contrasts: First, levels 1 and 2 are compared against 3 and 4. Second, level 1 is compared against level 2, as follows:               # Orthogonal contrasts for factor1:     Level1and2Vs3and4 <- c(-1, -1, 1, 1)     Level1Vs2 <- c(-1, 1, 0, 0)     contrasts(data$factor1) <- cbind(Level1and2Vs3and4, Level1Vs2)     cbind(Level1and2Vs3and4, Level1Vs2)     #      Level1and2Vs3and4 Level1Vs2     # [1,]                -1        -1     # [2,]                -1         1     # [3,]                 1         0     # [4,]                 1         0      Does one have to make contrasts for the other factors as well? If so, how can one make them orthogonal, when factor1 and factor2 have just two levels? Many thanks


title: Perform feature normalization before or within model validation?

text: A common good practice in Machine Learning is to do feature normalization or data standardization of the predictor variables, that's it, center the data substracting the mean and normalize it dividing by the variance (or standard deviation too). For self containment and to my understanding we do this to achieve two main things:   1. Avoid extra small model weights for the purpose of numerical stability.   2. Ensure quick convergence of optimization algorithms like e.g. Conjugate Gradient so that the large magnitude of one predictor dimension w.r.t. the others doesn't lead to slow convergence. We usually split the data into training, validation and testing sets. In the literature we usually see that to do feature normalization they take the mean and variance (or standard deviation) over the whole set of predictor variables. The big flaw I see here is that if you do that, you are in fact introducing future information into the training predictor variables namely the future information contained in the mean and variance. Therefore, I do feature normalization over the training data and save the mean and variance. Then I apply feature normalization to the predictor variables of the validation and test data sets using the training mean and variances. Are there any fundamental flaws with this? can anyone recommend a better alternative?


title: Calculating correlation coefficient

text: a <- c(10,11)     b <- c(2,3)      The covariance (`C`) between`a` and `b` =               sum ( (a-mean(a)) * (b-mean(b)) ) # equal to  cov(a,b) in this case as n=2      The sample covariance (`Csamp`) between`a` and `b` =               sum ( (a-mean(a)) * (b-mean(b)) ) / n-1 # equal to  cov(a,b)      Now lets say the sum of square differences for `a` = `SSa`. This means the Variance(`V`) for `a` = `SSa/n` and the sample variance(`Vsamp`) is `SSa/n-1`. The sample standard deviation (`SDsamp`) for `a` would be `sqrt(SSa/n-1)` The correlation coeffictent equations as I have found it are:               (1) C  /  sqrt(SSa) * sqrt(SSb)          # or          (2)  C  /  SDsamp_a * SDsamp_b          # (1) uses sqrt(SSa), not V or V samp and (2) uses SDsamp  not V and does not use Csamp      I am confused over how to calculate the correlation coefficient (`CF`) and the sample correlation coefficient (`CFsamp`). The two formulas above do not seem the same to me. can someone explain how to calculate `CF` and `CFsamp`?


title: Entropy of (Sum of Gaussians) versus Sum of (Entropy of Gaussians)

text: **Short version:** How can the joint entropy of two independent variables be less than the sum of those independent variables? The joint entropy should encode all information that a scalar function can, right? **Long version:** Assume there are 2 independent normal random variables $X, Y$ both with mean $0$ and variance $\sigma^2$.   1. We know that entropy of $X$ and $Y$ is $H(X)= H(Y)= \ln(2\pi e \sigma^2)/2$ (derivation)   2. The variance of the random variable $SUM = X + Y$ is $2 \sigma^2$   3. 1 and 2 mean that $H(SUM)= \ln(2 \pi e (2 \sigma^2))/2$   4. The sum of the entropies of 2 independent random variables is the entropy of their joint distribution, i.e. $H(X, Y) = H(X) + H(Y)$ . This implies that in this particular case $$H(X, Y) = (\ln(2\pi e \sigma^2)/2) \cdot 2.$$   5. Now note that if $\sigma^2=(\pi e)^{-1}$, then from 3 and 4 $$H(X,Y)= (\ln(2)/2) \cdot 2 = H(SUM)= ln(2 \cdot 2)/2 = \ln(2).$$ And if you increase $\sigma$ then $H(SUM) > H(X,Y)$. It seems quite fantastic that $(\pi e)^{-1}$ is an entropy tipping point for Gaussians. Do you know of any papers or books that make this observation? Is $N(0, (\pi e)^{-1})$ discussed as an alternative to $N(0,1)$ because of its neutrality in this context? And why is this happening at all? Shouldn't joint entropy be greater than entropy of any scalar function, it seems to be more general than any scalar function?


title: Comparing shape of time series data whist ignoring scale

text: I have a number of sets of time series data. Specifically each set of time series is associated with a single study subject. Each series measures the area between the vocal chords over a breath period (averaged over 10 breaths). Series are recorded for each patient under different dosages of a drug. I would like to, for each subject, quantitatively compare the curves produced under each drug dosage to see how similar they are. By visual inspection it appears that the graphs are very similar, however scale does vary somewhat, thus I would like to ignore scale if possible. I've been doing my statistics so far in R and Ruby, however Python would be fine also.


title: What is the difference between "likelihood" and "probability"?

text: The wikipedia page claims that likelihood and probability are distinct concepts. > In non-technical parlance, "likelihood" is usually a synonym for > "probability," but in statistical usage there is a clear distinction in > perspective: the number that is the probability of some observed outcomes > given a set of parameter values is regarded as the likelihood of the set of > parameter values given the observed outcomes. Can someone give a more down-to-earth description of what this means? In addition, some examples of how "probability" and "likelihood" disagree would be nice.


title: Posterior and Likelihood probabilities meaning

text: I am a computer scientist, so I have a background at maths (however limited). I am reading about posterior distribution from here http://en.wikipedia.org/wiki/Posterior_distribution . It says there: The posterior probability is the probability of the parameters θ given the evidence X: p(θ|x). It contrasts with the likelihood function, which is the probability of the evidence given the parameters: p(x|θ). My question is firstly can you provide me with a very simple example to understand better these concepts? And in addition in machine learning, what we want isn't the probability p(Class|x)? Thanks a lot


title: What causes perfect prediction but no significant predictors in logistic regression?

text: I want to do a logistic regression with R. I have 18 continuous covariates and a sample consisting of 100 observations. When I enter all of the covariates into the `glm()` model, none of them are significant, but the model predicts the outcome perfectly on the test data! My questions are:   1. Is the sample size large enough for running `glm()` with this many covariates?   2. What might be other causes of the problem?   3. How can I run such a model properly?


title: Using self organizing maps for dimensionality reduction

text: Over the past few days, I have been conducting some research on self organizing maps for a project at school. I have come to understand that self organizing maps can be used to reduce the dimesionality of your data. However, I do not understand how this works. For example, say you have a 10x10 network of neurons in a SOM, and your input is 25-dimensional. So, by my understanding, you would create a feature vector for each neuron that is also 25D. By the time training is done, you end up with 100 25D vectors. How is this exactly reducing the dimensions of the data? Am I supposed to be concerned with the location of the neurons? Thanks! EDIT: I've already read the question at Dimensionality reduction using self- organizing map but I dont feel it answers the question that I have.


title: Applying Kaplan-Meier survival function estimate to get expected number of events

text: I have the Kaplan-Meier estimate for the survival function which I obtained using R's survival package:               data.surv <- Surv(data$days_to_event, data$event_flag)     data.fit <- survfit(data.surv ~ 1, data = data)      The data are censored at 90 days after the initial signal so the survival function has 90 points and its value at the last point is equal to the ratio of the total number of non-events at 90 days after the signal to total number of observations in the data set. This value is far from 0 since the signal causes the event to happen rather infrequently. Is that to be expected? Now I have another data set where each observation has had the signal in the past (within 90 days) but no event so far. What is the proper way to apply the survival function estimate to the data to get expected number of events for some horizon? I thought $1-S(t_{current} + t_{horizon})/S(t_{current})$ for each observation is supposed to give me the probability, then adding them up would give the total count of the events expected from now to the horizon but even applying this to the original data set gives result which is way too low. It feels like there must be some mistake but I cannot find it. Could it be because of censoring (the majority of observations is censored - if the event does not occur within 90 days, the signal is thought to have lost its ability to cause the event)? Thanks!


title: Standard error of the residuals for a non-linear model

text: Hi I am new to R and statistics and used to linear models. Can you please explain the output? I used it to make a growth curve.               Formula: length ~ a * (1 - exp(-c * est_age))          Parameters:        Estimate Std. Error t value Pr(>|t|)         a 1.097e+03  1.026e+01 106.966  < 2e-16 ***     c 1.539e-01  1.982e-02   7.765 2.33e-09 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 41.74 on 38 degrees of freedom     Number of iterations to convergence: 6      Achieved convergence tolerance: 7.32e-07


title: Sample size to detect exponential distribution

text: I'm not a statistician, so I hope this question makes sense. I'm an EMT and we have patients known to be "frequent flyers". Example, we have a diabetic that we have seen 100+ times in 3 years. I have looked over his pattern of calls and it looks like they follow an exponential distribution, which makes sense to me. However, some times we will see him in clusters. Medically, this does make sense too. If he has a fever he will burn through his glucose more quickly and he may deviate from his baseline "lambda". I want to test the hypothesis that it is possible to 1) identify a new patient as a "frequent flyer" when N calls approximate an exponential distribution. and 2) identify when an established FF deviates from a baseline (indicating some change in patient status) A common FF might be that grandma starts falling out of bed more frequently, so this process might not be stationary. Basically, I think I'm asking what statistical test do I need to take a set of data and show that it is exponentially distributed? I've looked at Fisher Exact, chi Square, and K-S. Also is there a way to determine a "good" value for N? In my case N is the number of times I see the patient, so I would like that to be small if possible, like under 10 so I can flag them for additional non-EMS interventions as soon as possible.


title: How to analyse environmental temperature datasets

text: I need some advice regarding my statistical analyses with SPSS. I have a dataset of temperatures recorded automatically every 30 minutes in three different habitat types for the same time frame. I would like to: a) assess the differences between these 3 habitat types. b) compare these 3 habitat types with another dataset comprising temperatures recorded with data loggers that were attached to study 9 animals for the same time frame. Can you please give me some advice which tests I need to perform?


title: Negative values in predictions for an always-positive response variable in linear regression

text: I'm trying to predict a response variable in linear regression that should be always positive (cost per click). It's a monetary amount. In adwords, you pay google for clicks on your ads, and a negative number would mean that google pays you when people clicked :P The predictors are all continuous values. The Rsquared and RMSE are decent when compared to other models, even out-of-sample:                 RMSE        Rsquared      1.4141477     0.8207303      I cannot rescale the predictions, because it's money, so even a small rescaling factor could change costs significantly. As far as I understand, for the regression model there's nothing special about zero and negative numbers, so it finds the best regression hyperplane no matter whether the output is partly negative. This is a very first attempt, using all variables I have. So there's room for refinement. Is there any way to tell the model that the output cannot be negative?


title: Using Resampling to understand a large table

text: I have a data set that is very large. The attributes (columns) are several thousand. Some are sparse others are not. Some are ordinal, others interval, nominal or ratio. The row size is 10s of millions of rows. The data set is evolving -- that is the business is not sure which rows should finally end up in the set -- so they are developing their criteria -- what that means is that a record could be in the data today -- not in it tomorrow and back in the next day -- this is a work in progress. I have access to the data over a limited connection -- Limited as in limited tools I can use, limited computation power, limited time, etc. My first problem is that I want to be able to estimate for each column:   1. How many records have the field filled in or not (so null or populated)   2. And then try to estimate the field the value of the field. What I'm thinking of is deriving subsets of the data using re-sampling and then estimating the counts (item "1"). Is there a formula that can help me understand how how many samples I must take and what confidence the counts will have. Thanks.


title: Model matrix of a mixed model

text: Suppose I have a mixed model like this:               set.seed(123)     require(lme4)     df<-data.frame(id=rep(LETTERS[1:3],each=4),                    days=as.integer(rnorm(12,100,10)),                    event=rbinom(12,1,0.3))     df<-df[with(df,order(id,days)),]     df$event[sample(1:12,3,T)]<-NA        id days event     1   A   94     0     2   A   97     1     4   A  100    NA     3   A  115    NA     8   B   87     1     5   B  101     0     7   B  104     1     6   B  117     0     9   C   93     0     10  C   95    NA     12  C  103     0     11  C  112     0          fit<-glmer(event~days+(1|id), data=df, family="binomial")          mm<-model.matrix(fit)        (Intercept) days     1            1   94     2            1   97     8            1   87     5            1  101     7            1  104     6            1  117     9            1   93     12           1  103     11           1  112     attr(,"assign")     [1] 0 1      What I want is the full model matrix including the ones with missing responses, namely dim(mm) is (12 , 2). I have tried mm<-model.matrix(fit, na.action=na.pass), but it still returns the same thing. Would somebody here know how to realize it?


title: What is your favorite "data analysis" cartoon?

text: This is one of my favorites: ![alt text](http://imgs.xkcd.com/comics/correlation.png) One entry per answer. This is in the vein of the Stack Overflow question _What’s your favorite “programmer” cartoon?_. P.S. Do not hotlink the cartoon without the site's permission please.


title: Given a sample of random variables, and n, how do I find the ecdf of the sum of n Xs?

text: I can't fit $X$ to a common distribution so currently I just have $X \sim \text{ecdf(sample_data)}$. How do I calculate the empirical distribution of $\sum(X_1 + ... + X_n)$, given $n$? $X_1$ to $X_n$ are iid. My sample size is around $6\times10^7$. $n$ needs to be around 10-20. To give an example of what I am doing: Suppose $X$ is the age of a person. If I put 10-20 people in a room, what is the distribution of the sum of ges in the room? I have lots of peoples ages, and no reason to fit the distribution of an age to any of the common distributions. (previously on StackOverflow - suggested to move here).


title: Form of $\chi^2$ test

text: This article uses the following form of the chi-squared test: $\chi^2=\sum\limits_{i}\frac{(m_i-n_i)^2}{m_i+n_i}$ where $m_i$ is the number of counts of bin $i$ in a **modeled histogram** and $n_i$ the equivalent for the **observed histogram**.   1. How does one obtain this form of the $\chi^2$ test?   2. Is this a non-parametric statistic or am I assuming a Gaussian distribution of the data? * * * **Edit** : so, based on @Michael's answer, am I to assume that the statistics used in that article are in fact **wrong**?


title: Nonparametric test to compare link prediction algorithms

text: I have some experimental data, and I want to use a non-parametric test to verify my results. I am considering link-prediction on complex networks. I have a complex network $N$ and an algorithm/model $A$ which finds structure in $N$. I then remove some of the edges in $N$, first $10\%$, then $30\%$ up to $90\%$. For each of these states of $N$, I use $A$ to infer the structure, and then compute and AUC-score. I repeat this procedure 3 times, so I get a total of 15 scores. Now I take a new algorithm $A'$, which is a modification of $A$, and perform the same tests for $A'$. I now have 30 scores, 15 for $A$, and 15 for $A'$.                    10%   30%   50%   70%   90%     A1    x     x     x     x     x     A2    x     x     x     x     x     A3    x     x     x     x     x     A'1   x     x     x     x     x     A'2   x     x     x     x     x     A'3   x     x     x     x     x      Now I want to use a non-parametric test (probably a rank-sum test), to verify that $A$ is better at inferring structure in the networks. How do I go about doing this? I am using MATLAB for my experiments, so a call to `ranksum($x,y$)` is the solution to "how to do the test". What I am looking for here, is the arguments $x,y$, and if there are any things I need to consider about independence? There seem to be two obvious possibilities for the input: $$x = [A_1 (10\%), A_2 (10\%), A_3 (10\%), A_1 (30\%), \ldots , A_3 (90\%)]$$ $$y = [A'_1 (10\%), A'_2 (10\%), A'_3 (10\%), A'_1 (30\%), \ldots , A'_3 (90\%)]$$ or $$x = [A_1 (10\%), A_1 (30\%), \ldots, A_1 (90\%), A_2 (10\%), \ldots, A_3 (90\%)]$$ $$y = [A'_1 (10\%), A'_1 (30\%), \ldots, A'_1 (90\%), A'_2 (10\%), \ldots, A'_3 (90\%)]$$ Are any of these correct, and if only one of them, why?


title: Maximum entropy inference for k-means clustering

text: I am taking a course on Maximum Entropy Inference (MEI), where its application to k-Means was discussed. I am confused about the problem setting. As far as I understand, our goal is to find clustering solutions $c(\cdot)$ with small expected cost $\mathbb{E}[R^{km}(c)] \le R_{max}$. We are therefore interested in sampling from a distribution $p^*(c)$ over these solutions where good solutions have high probability and bad solutions have low probability. In order to be as 'agnostic' as possible, we ask for the distribution which maximizes the entropy, and we thus get: $$ p^*(c) = \arg \max_{p(c)} {H(c)} \quad \text{ s.t. } \mathbb{E}_p[R^{km}(c)] \le R_{max} $$ Obviously, $c(\cdot)$ are random variables in this scenario. But what is the underlying proability space? If $p^*(c')$ is large, it tells that $c$ is a good clustering solution, but... **for what problem?** Is this for a specific dataset $\mathcal{X}$ or for a large variety of such datasets? In either case, where does this or do these datasets play into the above equations? **EDIT:** I hope I am not being too vague here. Unfortunately I can't express myself very clearly as I don't understand the topic that well. However, I might be able to answer specific questions less ambiguously.


title: Overlap Probability of Empirical Distributions

text: I have a load versus capacity problem, and I'm trying to determine the likelihood of failure. Here's a simple example of what I mean. I have load and capacity discrete pdfs that were determined through a series of computer calculations that analyzed different input values (Monte Carlo simulations to address uncertainties). The figures below are a simplified version of the output: ![Load Distribution](http://i.stack.imgur.com/LoVIt.png) ![Capacity Distribution](http://i.stack.imgur.com/KRGZa.png) As you can see, there is some amount of overlap between the two. I am trying to determine the probability that the system will fail (i.e., the load with exceed the capacity). Is this as simple as integrating the overlapping area of the empirical cdfs and multiplying, or am I looking for something more complicated, like a Bhattacharyya coefficient? Thanks in advance for any advice!


title: Chi squared test with expected frequencies coming from another observation

text: I have done an experiment to test the performance of a system. In this experiment I collected answers from people. These answers are categorical (they were able to select on a 5 points scale). Each subject participated to 2 conditions: in one they are giving evaluations for observations under an ideal situation and in the other they are giving evaluations for observations coming from the real results of my system. My hypothesis is that the distribution of the ideal condition fits the distribution well under the real condition, so I want to demonstrate that the results of my system are perceived equally good as in the ideal condition by users. My results for the 5 points of the scale are the following: Observed frequencies (real condition): {16, 112, 42, 308, 100} Expected frequencies (ideal condition): {9, 81, 53, 340, 95} As I saw that most of the time the chi squared goodness of fit test is used with an expected frequency that comes from laws and not from empirical observations, I was wondering if:   1. My approach of using the chi squared test could be good for my aim, or if there is another better tool;   2. If it is okay to use empirical frequencies given by the observation of the same population under another condition.


title: Testing threshold cointegration in vector error-correction models

text: In Hansen and Seo's paper on Testing two regime threshold cointegration in VECM (J. Econometrics, 2002; 110:293), the authors proposed a test based on Lagrange Multiplier for testing treshold in model. However, suddenly the model denotes heteroskedasticity. Why the model denote heteroskedasticity?


title: Regression to obtain autocorrelation measure (AR(1))

text: This is not homework. I am a frequent user on math.stackexchange, but I am learning a bit about time series models and came across this example. Any ideas would be greatly appreciated. A linear regression model ￼was fit to some time-series data by ordinary least squares. The residuals￼ from the fit were then used to create two new variables, namely $E$ with values $\hat{e}_2,...,\hat{e}_n$￼ and $E_1$ with values￼ $\hat{e}_1,...,\hat{e}_{n-1}$. A linear “regression through the origin” was then run with E as the dependent variable and $E_1$ as the predictor. The slope estimate was 0.412 with a standard error of 0.133. Assume that the $e_t$ follow a standard AR(1) model.   1. Estimate the first order autocorrelation $\rho$ of the AR model.    2. Can the output be used to obtain a valid standard error for the estimate in 1? Is the answer to #1 above just the slope of the regression `E~E_1`?


title: Tests needed for repeated measures within groups 2*2 design?

text: I've recently conducted an experiment with two independent variables: musictype and wordtype, measured by response time. 20 participants. I've done a one-way repeated measures ANOVA along with the usual tests of normality, means etc. Do I need to include a t-test? Am I missing any other important tests? Many thanks in advance!


title: Pearson correlation coefficient test: low r and low p

text: I am running a correlation coefficient test, and the results were: $r = 0.382, p = 2.76 \times 10^{-13}$. So the $r$ value is not that impressive (usually we see $r>.5$), but the $p$-value is still significant. Usually I would think a low $r$ value would mean high $p$-value (no significant correlation), or vice versa (low $p$-value would mean a high $r$ value). Could anyone please explain what this means?


title: When to use generalized estimating equations vs. mixed effects models?

text: I have been quite happily using mixed effects models for a while now with longitudinal data. I wish I could fit AR relationships in lmer (I think I'm right that I can't do this?) but I don't think it's desperately important so I don't worry too much. I've just come across generalized estimating equations (GEE), and they seem to offer a lot more flexibility than ME models. At the risk of asking an over-general question, is there any advice as to which is better for different tasks? I've seen some papers comparing them, and they tend to be of the form: "In this highly specialised area, don't use GEEs for X, don't use ME models for Y". I haven't found any more general advice. Can anyone enlighten me? Thank you!


title: Comparing logistic coefficients on different logit models?

text: I've developed a logit model to be applied to six different sets of cross- sectional data. All the models have same independent variables but the sample size is different for each model. What I'm trying to test is the impact of each independent variable on dependent variable across all the models say for example income in model 1 vs. income in model 2 vs. income of model 3 etc. I have observed from the Odds Ratios, that these variables have different impact on dependent variable for different household categories and regions. My question is: How do I statistically test this change? **Update:** Further information regarding my question: I'm trying to look at the effects of IVs on DV at different disaggregated levels. The models are: $$ y_1 = \beta_{01} + \beta_{11}x_1 + \beta_{21}x_2 + \beta_{31}x_3 + ... \\\ y_2 = \beta_{02} + \beta_{12}x_1 + \beta_{22}x_2 + \beta_{32}x_3 + ... \\\ y_3 = \beta_{03} + \beta_{13}x_1 + \beta_{23}x_2 + \beta_{33}x_3 + ... \\\ y_4 = \beta_{04} + \beta_{14}x_1 + \beta_{24}x_2 + \beta_{34}x_3 + ... \\\ y_5 = \beta_{05} + \beta_{15}x_1 + \beta_{25}x_2 + \beta_{35}x_3 + ... \\\ y_6 = \beta_{06} + \beta_{16}x_1 + \beta_{26}x_2 + \beta_{36}x_3 + ... $$ The question is can I test whether $\beta_{11}$ is statistically different from $\beta_{12}, \beta_{13}, \beta_{14}, \beta_{15}, \text{ & }\beta_{16}$ etc.?


title: Conditional Probability equal to zero

text: Let $X_1$ and $X_2$ be two independent Normal variables with mean $\theta$ and variance $1$. Let $S=X_1+X_2$ and consider the following conditional probability: $$P \left[ X_1 \leq x_1; X_2 \leq x_2 |S=s \right]$$ If we assume that $x_1 \leq x_2$ and $x_2 \leq s$, then why does this conditional probability equal zero? I have been struggling to notice the impossibility here. Thanks.


title: Classification of biological samples with binary variables

text: **I have a question about a very basic classification problem.** The objective is to discriminate between feature A and B. Where A is a specific biological sample, and B all others. And the goal is to classify new unknown samples. To do so I have 12 variables with binary data, representing the presence or absence of certain traits. For example: 6 (3 known A and 3 known B) samples with 4 variables.                      Feature   Var1. Var2. Var3. Var4.       Sample1    *A*       1     1     1     1       Sample2    *A*       1     0     1     1       Sample3    *A*       0     1     1     0       Sample4    *B*       1     0     0     0       Sample5    *B*       0     0     1     0       sample6    *B*       0     1     1     0        What kind of classification algorithm can I use to predict a new unknown sample?


title: A buggy but effective feature?

text: Through error analysis I found that a quite effective feature actually has bugs in its implementation. Correcting the bugs actually decreased the classifier's performance. What do you do? Correct the bugs, or leave them there?


title: Is it possible to have a variable significant in multiple regression but not significant in stepwise regression?

text: I have run a stepwise regression and found that some of the selected variables are not significant yet in a multiple regression with all variables included in the model those variables were significant. What can be the cause of this and which results would you recommend that I use?


title: Checking for differences in time trends

text: I have two waterbodies that were sampled for a single water quality variable over time. Both waterbodies were sampled in the fall of each year from 2009-2012. In each year of sampling, 5 replicate samples were taken (so I would have 20 samples in total for each waterbody). I want to compare the two waterbodies over time to see if the water quality variable in each waterbody follows a similar trend. I feel like this is similar to an ANCOVA problem where a person could test for equal slopes. But this problem is a bit different because instead of having two continuous time series, I have 5 replicates of samples for each waterbody in each year. Anyway I'm wondering if ANCOVA is not reasonable, what is the best method to proceed with? A colleague suggested planned linear contrasts might help but I don't understand how that would be applicable to comparing 2 different waterbodies.


title: Posterior predictive checking in Hierarchical logistic model

text: I'm fitting a Bayesian Hierarchical logistic model with a random intercept, having data at two levels: subjects and hospital, logit$(p_{ij})=a_{j}+\beta_{0}+\beta_{1}X_{1i}+\cdots+\beta_{k}X_{ki}+\epsilon_{ij}$ $a_{j}\sim\text{Normal}(\hat{a}_{j},\sigma_{a})$ $\hat{a}_{j}=\gamma_{0}+\gamma_{1}G_{1j}+\cdots+\gamma_{m}G_{mj}+\epsilon_{j}$ For the "subjects" equation I have 3-5 variables, and for the "group" ones, 5-8 variables. I've fitted 6 models including or not some variables, adding or not random effect, and adding or not a model for the intercept. My question is about the posterior predictive checking: ![](http://i.stack.imgur.com/zw3hR.jpg) For each model (in the picture, only the 6th),in black there is the proportion of 1's in the response variable y (is a status variable with only 2 possible values 0 or 1) for each centre; the red points come from:   1. Sample from the posterior $p(\beta's,a's,\sigma_{a}|y)$   2. Simulate from the posterior predictive $p(y^{rep}|\beta's,a's,\sigma_{a},y)$.   3. From the posterior predictive I get $y^{rep}_{ij}\in[0,1]$, but my data $y$ consists only of 0's & 1`s.   4. If $y^{rep}_{ij}\geq 0.5 \rightarrow y^{rep}_{ij}=1$ else $y^{rep}_{ij}=0$. Then for each centre $j$ I have a vector of 0's and 1's (as in original data $y$) and I compute the proportion of 1's and I get 1 red point, then I repeat 30 times to get 30 red points for each centre. Is the 4th step statistically correct? If not, there is a way to go from the $p_{ij}$ to the data $y$ in order to compare? Thank you in advance.


title: Interpreting SPSS mixed linear model output

text: I have a problem interpreting the output of the mixed model procedure in `SPSS`. I found a couple of threads dealing with similar problems, but none helped me solve it. I have a 2x2 repeated measures crossover design with two fixed factors (medication (A/B) and genotype (A/B)) and a random factor (timepoint (1/2)). As I understand it, the output for the fixed effects is the general influence of the factors over all subgroups and the parameter estimates test the differences between levels. As I only have two levels for all my factors of interest the output should provide information about the difference between all levels from one that is used as 'baseline'. This is the output:               Fixed effects Typ IIIa     source      df  df      F       p     Const       1   2,798   1753    0,000     gen*med     1   30,945  0,988   0,328     genotype    1   31,329  4,564   0,041     medication  1   30,953  12,886  0,001     a. Dependent variable: RT.          Parameter estimates     Parameter           est             sd      df      t       p       conf. int. 95 %     constant            385,695463  10,256203   4,450   37,606  0,000   358,320557  413,070370     [med=B] * [gen=A]   -6,220872   8,462102    40,021  -0,735  0,467   -23,323140  10,881395     [med=A] * [gen=A]   -19,875937  8,460116    39,989  -2,349  0,024   -36,974613  -2,777261     [med=B] * [gen=B]   7,742406    4,340784    30,949  1,784   0,084   -1,111266   16,596079     [med=A] * [gen=B]   0b  0   .   .   .   .   .     [gen=A] 0b  0   .   .   .   .   .     [gen=B] 0b  0   .   .   .   .   .     [med=B] 0b  0   .   .   .   .   .     [med=A] 0b  0   .   .   .   .   .      **The questions are** : is it correct to interpret this output as that there is a difference between medication(A) in genotype(A) compared to medication(A) in genotype(B)? If that is the case then I do not understand how to compare the effect of medication within genotype(A). Am I doing something wrong? I assume this output is correct for multiple comparisons. Is that correct? The code is:               MIXED RT BY genotype medication timepoint     /FIXED=genotype*medication genotype medication | SSTYPE(3)     /METHOD=ML     /PRINT=SOLUTION     /RANDOM=timepoint | COVTYPE(VC)     /REPEATED=medication*timepoint | SUBJECT(id) COVTYPE(CS).      I would really appreciate if someone could help here :)


title: Split plot design with repeated measures in package (nlme) - model errors

text: Forgive me if this answer is already somewhere on this site, but after extensive searching I have found nothing to help me. Here is a brief overview of the experiment:   * I have 18 tanks, each of which are split into 4 sections, each section has a different specie of wetland plant. Tank is my whole plot factor. Variety, or species (spp), is my split plot factor.    * I have 2 treatments, each at three levels. Nutrient addition (nut) and clipping (clip)    * I am measuring methane (ch4) emission as my response variable. I measured methane on 2 occasions, represented by variable t for time.  **Data Structure**               >str(ch4data)     'data.frame':   144 obs. of  6 variables:     $ tank: Factor w/ 18 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...     $ spp : Factor w/ 4 levels "Emac","Ewal",..: 4 2 3 1 3 4 1 2 1 4 ...     $ nut : Factor w/ 3 levels "1","2","3": 2 2 2 2 1 1 1 1 3 3 ...     $ clip: Factor w/ 3 levels "a","b","c": 3 3 3 3 2 2 2 2 2 2 ...     $ t: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...     $ ch4 : num  0.382 1.642 1.529 0.245 11.482 ...      **Model**               ch4.model=lme(fixed = ch4 ~ spp + nut + clip + t + spp*nut + spp*clip + spp*t +                    nut*clip + nut*t + clip*t + spp*nut*clip + spp*nut*clip*t +                    (nut*clip)%in%tank, random= ~1|tank, data=ch4data)      this gives the following error               Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at            level 0, block 1      Specifically, my questions are is: **What is wrong with the fit of my model and what does this error mean?** Additional help with fitting this model correctly would be much appreciated. I am new to R, so your response should take that into account.


title: Repeated measures ANOVA within and between factors

text: I have conducted repeated measures ANOVA, with 2 factors: One within factor (`time`, pre-post intervention) and one between factor (`treatment`, two groups). I have measured different dependent variables. I find in all the cases an effect of the within factor, an effect of the between factor, and only in some cases an interaction between both. I want to know how to interpret these results.   * What does the interaction mean?    * Is the significance of the factor `treatment` (between) sufficient to be able to say that one group is superior to other one, that one treatment is better?  EDIT: The measures I took in three moments are: time spent using technologies, frecuency of use, abuse.. The significant interaction is for within factor (moment) X between factor (group), but it's only significant for the time spent using technologies. What means the sugnificance of between factor? Why it's important the intetaction? Thanks * * * The measures I took in three moments are: time spent using technologies, frecuency of use, abuse.. The significant interaction is for within factor (moment) X between factor (group), but it's only significant for the time spent using technologies. What means the sugnificance of between factor? Why it's important the intetaction? Thanks


title: detecting plagiarism on multiple-choice test

text: Suppose an invigilator suspects one student of copying answers off another student's paper during a multiple-choice exam. She later checks their answers and finds some similarities—but on the other hand, there are bound to be similarities given the nature of the exam. How should she go about determining whether her suspicions were founded? In other words, she will surely have to compare the exams to those of other students (who, let us assume, were not cheating). But if the class size is very large, is it reasonable to take a random sampling for comparison? How many would she then take? If there were many questions on the exam, would it also be reasonable to take a sampling of the questions for comparison? Does it make a significant difference whether each question had 2 possible answers (true/false) or, say, 4? I don't have any specific numbers because I'm wondering about how this would work in general. I have a background in mathematics but little training in statistics. How would you describe this analysis in statistical terms? Thank you.


title: Non-linear model robustness

text: I am building a non-linear model aiming to describe the mechanistic process of resource allocation. There several terms, and what makes the model non-linear is competition between lines that are facing each other pair by pair. I tested the significance of terms through AIC comparisons, deleting non significant terms (it is not really hierarchical, I explored a big set of pathways). I have the feeling that sometime, removing one term can change a lot the value of parameters for other terms, so that the values of parameters seem to be completely inter-dependent. I have 2 questions:   1. How can I formally assess whether this is true (that parameter values change a lot when removing terms)? Please note I have a huge number of parameters, so comparing them one by one is very difficult.    2. If it is true, is it due to the intrinsic nature of non-linear model, or is it due to some kind of colinearity in my terms? Standard error are very small, but when comparing confidence intervals of parameters for two models differing by only one term, these CI do not overlap at all. Thanks!


title: Is there a multiple-sample version or alternative to the Kolmogorov-Smirnov Test?

text: I am comparing the size distribution of trees in six pairs of plots where one plot received a treatment and the other a control. Using a Kolmogorov-Smirnov test on each pair of plots I find that $p$ ranges from $0.0003707$ to $0.75$. Are there any appropriate methods for dealing with all of the replicates together, such as a multi-sample extension of the KS test, or is there an appropriate follow up test? Or should I just conclude something like "The size distribution differ significantly $(p < 0.05$) in 2 pairs of plots and marginally ($p = 0.59$) in one pair of plots."


title: How to compare models from different but related datasets?

text: I'm building regression models on four the different but related data set and at the end, I want to test the significance of models. Since my models are built in a different data set, it's not comparable. But there are some hierarchy in my dataset. 1) The output in my main dataset is the sum of the outputs of data set A, B and C 2) The value of feature is the same in all data sets, but all features are NOT present in all data set and some data sets only having a subset of the features of the main data set. 3) The union of the features of the dataset A, B, C are the features of the main data set. So, I want to build a regression model on these data sets seprately and compare the performance and significance of the models togethers.Becuase the dataset are not the exactly the same for all four models, I can not use standard statstics to test the sigificance of models like AIC, CP- statstics,..... Is there anyway to compare the performance of the model ? Here is my dataset(I just included first two data points) for more clarification: ![enter image description here](http://i.stack.imgur.com/kAN5q.png)


title: Constructing hypothesis test $\mu=0$

text: My question is as below: Random sample $Y_1,\dots, Y_n$ from $N(\mu,\sigma^2)$. Let $\bar{y}=\frac{1}{n}\sum Y_i$.   (1) $(\bar{y}-\,\mu)/\sigma=a'x$ where $x\in R^n$, $x\sim N(0,1)$.   (2) $\sum(y_i-\bar{y})^2=\sigma^2 X'AX$ where $X \sim N(0,1)$.   (3) $\frac{n(\bar{y}-\mu)^2}{\sigma^2}=X'BX$ where $X \sim N(0,1)$.   (4) $Aa=0$,$AB=0$. I already solved questions (1)-(4) and know what $A,B$ and $a$ are.   The last question is given   "How can we use these results to find a test for the hypothesis $\mu=0$?" It is really hard to understand what the question wants and how to solve it.


title: Drawing pair of cards. Did my brother play fair or unfair?

text: **The Game** There is a pool of $n$ cards that are marked either by a `A` or by a `B`. There is a proportion $p$ of the cards that are marked by a `A` and a proportion $1-p$ of the cards that are marked by a `B`. The game is the following. I randomly draw a card and put it back in the pool. I then randomly draw a second card. These two first cards form a pair of cards. A pair can be either `AA`, `AB` or `BB`. I take note of this pair on my spreadsheet and replace the second card in the pool as well. I repeat the draws up to the time I drew $n$ pairs of cards to get $n$ entries on my spreadsheet. The Probabilities for each pair to be `AA`, `AB` or `BB` are $p^2$, $2p(1-p)$ and $(1-p)^2$ respectively. Therefore, the probability that at the end of the game I drew $x$ `AA` is $p^{2x} + (1-p)^{n-x}$. The probability that at the end of the game I drew $y$ `AB` is $2p(1-p)^y + (1-2p(1-p))^{n-y}$. The probability that at the end of the game I drew $z$ `BB` is $(1-p)^{2x} + p^{n-x}$. Everything seems correct up to this point? **Question** I ask my little brother to play this game and he can basically do two things. He either plays fair and chose to draw the cards randomly or he can play unfair and chose to draw the first card randomly and then draw the second card non-randomly so that he may end up with lots of `AA` and `BB` or with lots of `AB` for examples. If he plays unfair he may, for example, chose that every time he draws an `A` at first hand, his probability to draw another `A` to complete the pair becomes $p+\Delta p$ resulting to an excess of `AA` and a lack of `AB`. So he plays and he observes $x$ `AA`, $y$ `AB` and $z$ `BB`, where $x+y+z=n$ by definition. The values $x$, $y$, $z$ differ from the expectation of $p^2$, $2p(1-p)$ and $(1-p)^2$. **How can I tell with a 95% confidence interval (using the distributions of my second paragraph) whether my brother played fair or unfair?**


title: Statistical significance in yes/no poll question

text: I have no statistics background, but an trying to complete my first quantitative research survey writeup for my researcher in education class and hope someone can direct me. I created a poll of 62 teachers and am trying to analyze the following 4 questions   1. Do they use social media personally (yes/no)   2. Do they use SM in the classroom (yes/no)   3. If they plan to use in future in classroom (yes/no)   4. What type of secondary school do you work in (public/private) I am running into a challenge with how to complete a "test of statistical significance for this data" The professor suggested using either a t-test or ANOVA and provided links, but the data I have do not seem to fit into the formulas. Is anyone able to provide some guidance about how to test this type of survey? Is there a different way to analyze the data?


title: Representation within a RKHS framework

text: Given a p.s.d kernel $Q$, can minimization/maximization of $Tr(X^TQX)$ over X be represented within a reproducing kernel Hilbert space (RKHS) framework? If there is a primary concern with the trace function being unbounded(for maximization) or trivial zero-solution(minimization); you may consider constraints over $X$; for example-say $X$ being orthogonal. Am trying to see if I can fit this in an RKHS framework. Also, a few thoughts around this function using hilbert-schmidt norms are as follows: I do see that $TrX^TQX$ can be represented as $Tr[(SX)^T(SX)]=||SX||^2_{HS}$, using the hilbert schmidt norm where $S$ is the p.s.d square root of Q. (Ex: $S=U\lambda^{1/2}$, where $Q=U\lambda U^T$ is the eigen-decomposition of $Q$).


title: Unpaired t-test with only summary means

text: Is it appropriate to conduct an unpaired t-test with only summary means from both groups? For example, let's say you want to compare: a. Professor A's scores received on a evaluation by **his** students from one semester; with   b. The average evaluation scores among all professors in the university completed by **all** students from the same semester. The teacher evaluation has 15 Likert questions (6-point scale). Constraints:   1. You don't know how many students completed the evaluations in either group.    2. All you have is 15 mean scores for Professor A, and 15 mean scores for the overall university sample.  It's safe to assume the evaluation is psychometrically sound. An unpaired t-test is conducted to compare means from group 1 (15 means from Professor A's students) and group 2 (15 means for all professors in the University). So $N = 30$. Is this appropriate?


title: Constant hazard in survival analysis and random event occurence

text: I read that in survival analysis, when the hazard function $$\lambda(t)=\underset { dt\rightarrow 0 }{ lim } \frac { Pr(t\le T<t+dt|T\ge t) }{ dt } =c$$ then the event occurrence is **random** and independent (memoryless). I can understand the concept and implications of monotonically increasing hazard function. However is there any theory or proof that when the $\lambda(t)=c$ then the event occurrence is random and independent?


title: Bayes' rule for conference talk attendance

text: This seems like a typical Bayes' Rule question, but I'm having trouble formulating it properly. Say I'm a speaker on a conference tour which moves from city to city. I'm trying to estimate how many people will attend my talk in the next city, given how many attended in previous cities. For example, say the total number of attendees in each city is 100, and that there are five simultaneous, equally likely tracks, so I start in the first city assuming that 20 people will attend (i.e., my first prior is 0.2). In actuality, 30 people attend. How do I use that evidence to predict how many will attend in the next city? How about the one after that? It feels like I ought to be able to write a simple Bayes' Rule expression for this, but everything I try seems wrong. Any help would be appreciated.


title: Good statistical software system for teaching response surface design?

text: I am looking for a good software system to teach students on response surface designs with. Few things I can think of are SAS Proc RSReg or SAS ADX. Is there any other software?