title: Any suggestions for storing diverse data sets for on-line visualization?

text: We are trying to design a data repository for storing datasets, from many different sources, in such a way that they can easily be used in online data presentations. This data repository will be part of a website that reviews, interprets and discusses scientific reports and data on public health. It is targeted at scientists, medical experts and policy makers. A prominent aspect of the site will be data visualizations: charts and maps. Data used in visualizations must be presented in tabular format and be available for download, as well. We are considering to use a JavaScript library for the charts, possibly high charts. Maps will probably be mainly images with click area's, but some server site tool may be used for specific maps. The data that we use, is collected from other organizations, in the form of database dumps and ASCII files, or sometimes simply copied from scientific studies. We are considering to link to open data sources, in the future. Previously, we have created a generic data model in which each data set was stored in the same generic format: a data set was linked to a variable and for each measure the following dimensions were available: period, area, age, sex and social economic status, all of them optionally present in a data set. Finally, the measures were stored as floating point numbers. As an example, think about _mortality of coronary heart diseases_. A record in the dataset would be the absolute number of deaths in the year 2014 for Dutch males in the age group 65 to 75. Meta data describing the complete set, certain domain values of a domain or individual records were supported. Domain values, like 2011, male and 65-75 were shared between data sets. This allowed for searching across data sets and, to some extent, facilitated combining different data sets in one visualization. This data model, however, required extensive preparation of data before it could be imported, including un-pivoting and mapping of domain values. Another aspect was versioning of data sets. In some cases, the most recent available data should be used, in other cases, a fixed set is to be used, even if more recent data is loaded into the repository. In that case, each visualization is linked to a specific version. This version is only upgraded if the text accompanying the visualization is reviewed and updated as well. Imagine a trend chart in which the latest years show some new trend. This newer data is to be shown only when the editor has done an analysis on that trend. So, this is not so much about large data sets, but rather about large numbers of diverse data sets. Characteristics:   * The presence of various measures. Sometimes, more than one measure was present, for instance an absolute number as well as a percentages or a standardized rate or confidance intervals.   * Various aggreagation levels in the dimensions: years, months, days;   5-year age groups or 0-65/65+, etc.   * Meta data and footnotes on sets, dimensions and/or single values. It must be possible to show these in tabular presentations with footnote symbols linking to descriptions.   * Versioning of data sets, where newer versions remove, correct and/or add values.   * Aggregation may or may not be trivial. For absolute numbers,   aggregation is a summation. This is easily performed and will be   supported. For standardized rates, aggregation is complex and outside the scope of the repository. In such case, pre-aggregated totals must be added to the dataset, to be able to show marginals in cross-table presentations. These totals must be distinguishable from base   numbers, as it must be possible to leave them out of visualizations. Imagine a stacked bar chart over years, where the series are age   groups. A total over all ages must not be included in the stacked   bars.   * We must be able to back-track data used in presentations to the   original source, to be able to explain were they came from. So   ideally, nowhere in the work-flow, data has to be copied and pasted nor edited manually. On top of the data repository, there will be a presentation composer. This composer must be capable of querying the repository. These queries should ideally be generated from choices the editor makes: filtering, selecting dimensions as x-axis, series, etc. Linking to a newer version of a data set should be supported, without completely redefining the presentation. At this point, we consider supporting various of the larger datasets with a dedicated data model, and some sort of adaptors for the presentation composer, mainly to avoid various issues with importing data sets in a generic model. Does anybody recognize this scenario and have some suggestions about best practices or maybe even an out-of-the box solution?


title: Are there differences in Bayesian and frequentist approaches to EDA?

text: Very simply put: Are there any differences in Bayesian and Frequentist approaches to Exploratory Data Analysis? I know of no inherent biases in EDA methods as a histogram is a histogram, a scatterplot is a scatterplot, etc., nor have I found examples of differences in how EDA is taught or presented (ignoring a particularly theoretical paper by A. Gelman). Finally, I looked at CRAN, the arbiter of all things applied: I haven't found packages tailored to a Bayesian approach. However, I thought CV might have a few people who could shed a light on this. Why should there be differences? For starters:   1. When identifying appropriate prior distributions, shouldn't one investigate this visually?   2. When summarizing data and suggesting whether to use a frequentist or Bayesian model, shouldn't the EDA suggest which direction to go?   3. The two approaches have very clear differences on how to handle mixture models. Identifying that a sample likely comes from a mixture of populations is challenging and directly related to the methodology used to estimate the mixture parameters.   4. Both approaches incorporate stochastic models and the selection of model is driven by understanding the data. More complex data or more complex models necessitates more time in EDA. With such distinctions between stochastic models or generating processes, there are differences in EDA activities, so shouldn't there be distinctions arising from different stochastic approaches? Note 1: I'm not concerned with the philosophies of either "camp" - I only want to address any gaps in my EDA toolkit and methods.


title: How to do a paired test when some cases have no pair?

text: I have a question regarding paired and unpaired tests. I know the difference between both tests. I am using R for Wilcoxon test. My test is a paired test, however the X and Y do not have the same length and R is giving an error. I do not want/think I should use unpaired test. For example, I have a subject that makes some cookies every hour. I use some special kind of treatment to increase the stamina. Before the treatment he can only work 4 hours (after day he gets tired) and produces some cookies. After the treatment he produces more cookies per hour and work 7 hours without getting tired. X contains the number or cookies for each hour before the treatment and Y contains the number of cookies after the treatment. X contains 4 values and Y contains 7 values. Now if I want to use paired test, R gives error. What should I do? Is there any solution or explanation for such kind of situations? Can I add just NA NA? **"This is just an example please do not point out mistake in the example, it is just to give you an example."** Thank you. ## EDIT Here is basic R script that I use               someData <- read.csv(file="cookie_data.csv",head=TRUE,sep=",")     wilcox.test(someData$X, someData$Y, paired=TRUE)          Sample Data:     X,Y     2,3     3,2     3,3     2,2     ,3     ,7     ,2      When I use this script, R does not give any error. However, when I print someData$X, it prints 4 values and after that it start writing NA NA NA. I noticed R automatically filled blank values with NA. This script gives me p-value but I do not know if it is correct.


title: How do I cluster data where some observation can not be in the same cluster?

text: Say I have 12 (x,y) positions:               data =         1.0000    0.8500         2.0000    1.1000         3.0000    1.2000         4.0000    0.9000         1.0000    2.0000         2.0000    2.1000         3.0000    1.9000         4.0000    2.0000         1.0000    3.1000         2.0000    3.1000         3.0000    2.9000         4.0000    2.8000      I want to classify these points into horizontal lines. From a scatterplot it is easy for a human to see that the points form 3 horizontal lines; one around y=1, another around y=2 and the last around y=3: ![enter image description here](http://i.stack.imgur.com/Tyr7f.png) **Which classification algorithm should I use to classify the points?** The easiest algorithm seem to be k-means. One problem with this is that I need to specify the number of clusters. How do one typically get around this problem? I was thinking that I may run it with an increasing number of clusters, n, and calculate the total square deviation from cluster mean and once the decrease in total square deviation relatively to the last number of clusters, n-1, is less than a certain threshold (say 1%) then n is a sufficient number of clusters. Another algorithm I was considering using is hierarchical agglomerate clustering. The solution here would form a dendrogram, and like for the k-mean I may choose a level with n clusters for which the total deviation is below some threshold criteria. One important piece of information that I do not understand how to use though in either of these algorithms is this; no cluster may have the same x position more than once! So far I have just ignored the x coordinate and clustered based on y values. **How can I use the criteria that no cluster may have the same x position more than once with a cluster algorithm?** Thanks in advance for any answers!


title: Independence of variables and ANOVA

text: I performed a test to study the joint effects of hardness and alkalinity (independent variables, IV) on reproduction (number of offspring, dependent variable, DV) of water fleas, but I’m not sure what statistical test is appropriate to analyze data. First, I’m not sure whether both IV are independent. They are often correlated in natural waters, thus I kept them autocorrelated in my experiments (they are positively related). However, I manipulated them independently (by adding different salts to water). I may say “they are independent but autocorrelated”?! So are they really independent? Second question: how to analyze data? The design is not factorial: I tested 4 levels of hardness, 4 levels of alkalinity (this would give 16 waters, but I tested only 10 waters). Graphically, data suggest an interaction between hardness and alkalinity on reproduction of the organisms, but I cannot test it. Supposing the assumptions for ANOVA are met, can I use ANOVA? I first thought of performing ANOVA following the GLM procedure and using simultaneously hardness level and alkalinity level as fixed factors (although they vary to a little extent, as a consequence of the measurement technique). But I was told to perform ANOVA (GLM) for hardness level and, separately, for alkalinity level, but I don’t understand why. Anyway, how should I analyze data?


title: how would I compare 2 linear models with the same response variable?

text: I two have linear models with the same response variable: (2) y = x1 + x2 + x3 and (2) y = x4 + x5 + x6 + x7 + x8 what would be a good way to compare them? I know that if (2) would be y = x1 + x2 then an ANOVA would do (in R). I checked all the questions, but could not find a related one. Edit: note that the F-test is about comparing two models that have not only the same response variable, but also the same predictors (see this handout for example ). Does anyone know if the F-test is also suitable for linear models that only have the (exact) same response variable?


title: How can you construct a social network empirically from multiple sources of data

text: Let's say you have databases $A$ and $B$, and a population of users $\left\\{u_1,\ldots,u_n\right\\}$. You're building a social network graph of $\left\\{u_1,\ldots,u_n\right\\}$. From database A, you get an edge weight feature $a(i,j)$ which connects user $i$ to user $j$. From database B, you get a similar edge weight $b(i,j)$. You want to build a single network incorporating all this data. I can think of a couple ways to do this.   1. Treat edges $a$ and $b$ separately, make a multigraph...   2. Combine the two edge weights together in some way. Suppose you define the edge weights $e(i,j) = a(i,j) + \alpha b(i,j)$, and $\alpha$ is a parameter.  I'm mostly interested in method $2$, but allow for the possibility that $1$ could work well. Is there a way to optimize the choice of $\alpha$? One idea would be: remove a small subset of the edges, and then try to predict where they would be (using a clustering/triangle approach, say), or fitting a likelihood value to where the removed edges are.


title: MCMC and degree of freedoms (small number of observations)

text: I want to estimate a simple Regression. I have 20 observations and 10 Regression-parameters. The degrees of freedom are too small to get reliable point estimates and p-values. I found the following Statement: "Another advantage of MCMC estimation is the problem of degree of freedom. In the maximum likelihood estimation, if the parameters number is large compared with the obervation, the model becomes unstable and the obtained result would not be reliable. But we do not need to worry about it if we use MCMC estimation. In this note, for example, the number of parameteres to be estimated is larger than n, which is the observation number." (Estimating Markov Switching model using Gibbs sampling with a statistical computing software R) Would bayesian inferenz (MCMC) of the Regression solve my problem? Or is the Statement only true in the specific case of the markov switching model?


title: Multilevel model and sample size

text: I want to estimate a complex hierarchical multilevel logit model with spatial structure. The model is a two-level-model. I have 20 regions (upper Level) which contain several districts (lower Level). The dependent variable is employment growth at district-level. The multilevel logit model should include at least three variables at the upper Level and eight variables at the lower Level. I am a little bit concerned that my sample size is too small to produce reliable Point estimates and Standard Errors for such a complex model. There are two Options: 1.) Choosing larger district Areas for which more relevant data is avaidable, resulting in: 20 regions with number of districts per region varying from 3 to 40 (Average number of districts per Region: 13) 2.) Choosing smaller district Areas, resulting in: 20 regions with number of districts per region varying from 8 to 400 (Average number of districts per Region: 60) I like to use larger districts, because more data are avaidable. Is the sample size too small for the multilevel model? In literature this topic is discussed controversially. Some authors say you need at least 50 Groups with an average number of 50 memebers per group in a multilevel model with dichotomous data. Ohter claim that at least 20 Groups with minimum 5 members per group are sufficient. Is there a test for minimum sample size for multilevel models (i use R by the way)? Would an MCMC-approach help me to produce more reliable estimates? Off Topic: Has someone try to incorporate spatial effects (spatial lag or spatial error) in multilevel models? Is that possible? Thans in advance


title: Application of likelihood ratio test to test the Markov property

text: Do you know a reference (freely available on the web) where the likelihood ratio test is applied in order to test for the Markov property? The setting is a directly observable discrete Markov-chain with given transition matrices. The concrete application is a model for credit rating transitions. Some statistical tests are applied in this paper Time series properties of a rating system based on financial ratios but there are too little details. Evaluating the Markov property by Frank Bickenbach and Eckhardt Bode is an example for testing the Markov property but it is behind a pay wall. EDIT: The concrete test is: Null hypothesis: rating transitions do not depend on past rating distributions. The alternative hypothesis: rating transitions depend on past rating distributions (in the sense of a 1st order Markov chain).


title: Is it okay to compare fitted distributions with the AIC?

text: suppose I have a data set $x_1, \ldots, x_n$ and I would fit a normal, an exponential and a uniform distribution to them. The fitting function spits out a bunch of goodness-of-fit statistics, e.g. the AIC, BIC, chi-square, Kolmogorov-Smirnov, etc. I am trying to convince someone that the AIC is not appropriate here, because we have different log-likelihoods, and sometimes different number of parameters, depending on the distributions. I would prefer the p-value of the Kolmogorov-Smirnov-Test to compare the fits. Is my approach justified? How can I convince my coworker the AIC is not okay here (he likes to see a cited paper or something equivalent)? Thanks in advance! **edit:** Specifically, I was shown this article: http://www.vosesoftware.com/whitepapers/Fitting%20distributions%20to%20data.pdf I have no idea what to say to this. Page 4 lists the flaws of the chi-squared, Kolmogorov-Smirnov etc, and page 5 and 6 praise the AIC. Is he right?


title: Analyse within-subject data as between-subject data

text: Presence (or better sense of presence, i.e. the extent to which one feels that a mediated experience is not actually mediated) is a notoriously difficult variable to pin down and the question of how to measure it is hotly debated in its respective community. One particular debate concerns the the validity of standard self-report questionnaires used to measure a subject's sense of presence. In a between-subjects experiment, people performed a search task in either a real or a virtual office. They then rated their sense of presence on a set of standard questionnaires. Participants in the virtual condition self-reported a similar sense of presence than participants in the real world condition. While the authors only conclude that presence questionnaires should pass a "reality test" in order to be useful, I personally think that the results might have been different in a within-subjects experiment. Presence is such a subjective and abstract concept to the layman that I think it unlikely that subject A's feeling of presence under condition X can be meaningfully compared to subject B's feeling of presence under condition Y, but that condition X gives a subject a frame of reference on which to evaluate condition Y, if the subject is exposed to both (Note, IMO this would entail that presence could only be regarded a measure for relative differences between conditions rather than an absolute scale). This is somewhat substantiated by the fact that answers on presence questionnaires have a notoriously high variance. I have already designed an experiment (which has four different conditions) as a within-subjects experiment using 48 subjects and balancing the order of conditions across subjects. Then I got to wonder whether I could test my assumption with my data. So I am thinking about analysing my data as within-subjects data (i.e., using a repeated-measures ANOVA) but also taking the very first condition each subject is exposed to and simply treat that subset of the data as if it came from a between-subjects experiment. However, I am wondering whether this is a valid statistical approach and what statistical test would be appropriate, rather than simply stating, for example, that in the within analysis condition X was different from condition Y while in the between it was not. For example, would it be valid to treat _mean for condition X using within data_ and _mean for condition X using between data_ (for all four conditions) as factors for an ANOVA to see if the mean for a condition significantly differs for a within and a between subjects experiment? NB: I could potentially add more subjects which only get exposed to one condition in order to equalise the number of data points available for the within and between subjects analysis but I would prefer to do the analysis on the original subjects only, as (a) the experiment is quite lengthy and cumbersome, and (b) I was actually hoping that keeping the same pool of subjects for both analyses would give me a stronger indication of whether my hypothesis holds, as it eliminates any variability due to some subjects appearing in one set of data but not the other.


title: I want to show some properties for Brownian motion

text: Let $B(t)$ be a Brownian motion. Show that the following processes are Brownian motions on $[0,T]$ 1) $X(t)=-B(t)$; 2) $X(t)=B(T-t)-B(T)$, where $T\lt \infty$; 3) $X(t)=cB\left(\frac{t}{c^2}\right)$; 4) $X(t)=tB(\frac{1}{t})$, $t>0$ ,$X(0)=0$. Thanks for help.


title: What is the meaning of Mean Sq for for quantitative variables in AOV?

text: In an open course, Mean Sq was explained using a factor variable. Suppose your linear model is               y = b0 + b1*I{x == A} + b2*I{x == B}      where A, B are Factor variables. It seems that               Mean Sq of b1 = Mean{(y - Mean(y))^2 | x==A} / Mean{(y - Mean(y))^2}      which means: Inter-class deviation / deviation of all But for quantitative variables, what's the meaning of Mean Sq ? Hope I have described the question clearly ...


title: What does a confidence interval with a negative endpoint mean?

text: I have 10 iid r.v. with Bernoulli distribution with $X_{i} = 1$ for a positive result. I'm given $\sum_{i=1}^{10} X_i= 1$ and need to find a two-sided 99% confidence interval for $\theta$. So $\alpha = 0.01$ and $z_{\alpha/2} = 2.575$ based on the standard normal distribution. The variance is given $\sigma^{2} = \hat{\theta}(1-\hat{\theta})$, with $\theta = 1/10$ from the given information. The confidence interval is: $$ \left(\hat{\theta} - z_{\alpha/2}\frac{\sigma}{\sqrt{n}}, \hat{\theta} + z_{\alpha/2}\frac{\sigma}{\sqrt{n}}\right) $$ so my confidence interval is $(-0.143, 0.343)$. Have I messed up a calculation somewhere? I'm not sure what this confidence interval tells me...perhaps that the sample size is too small to really get any useful information on $\theta$?


title: How to apply PCA with regression ?

text: I'm trying to predict y : from variable x1,x2,x3,.....x1000 say I got equation something like :               y = intercept  + aX1 + bX2 + cX3 + ... zX1000      Now I want to select only meaningful X How can I apply a concept of dimension reduction such PCA to my model ? say once I found that there are only 3 components say ( P1 P2 P3 ) that can explain 95% - how can I actually apply this back to the model ? Can anyone explain this to me ?


title: Summarizing fixed effects variables in Stata

text: Using a data set of 209 countries, I have run a fixed effects regression for non-OECD countries. When I summarize my variables, it gives me all of the variables in data set, but I only need a summary of those variables that have been used in my fixed effects regression. For instance, when I run fixed effects, the number of groups is 138. So, my total observations for a variable should be 138 multiplied by 26 years (1980-2006), which is only 3588 observations. How can I ask Stata to generate this?


title: Test data generation for neural network handwriting recognition

text: Can someone share some Octave/Matlab code or algorithm to pre-process a photo taken from mobile camera of a handwritten digit. After pre-processing, the data should have similar characteristics to the MNIST data set digit images. I have a neural network trained using the MNIST data set. Now I want to test my implementation by taking handwritten digit images using a phone camera and saving it on my computer. I want to give this image as input to test my neural net implementation. Thanks in advance !!


title: Using Zipf's law to select top K ngrams?

text: Lets say we extract all the ngrams (unigrams, bigrams, trigrams) along with their frequency from a text. Now if we want to select say some Top K ngrams, is there a way to estimate a good value of K using Zipf's law ? e.g. if Ngrams are -> w1:10, w2:9, w3:9, w4:8, w5:5, w6:3, w7:1, w8:1, w9:1, w10:1. Here are Ngrams are w1,w2,w3 etc and numbers after them are their frequencies. Ideally we would like to have K here as 5 or 6. Is there a way to estimate K using Zipf's law ?


title: Distribution of a simple random sample

text: Homework Question Hello, I have a simple homework question that asks: _We know that when we sample, with equal probability from a finite population {x1, x2,..., xN} **without replacement** , we obtain a simple random sample instead of an (iid) random sample. Let X1,...,Xn (n < N) be the sample obtained. It is claimed in class that X1,..., Xn_, though no longer independent, are still identically distributed according to the pmf: f(x) = _1/N, if x = xi, for some 1 <= i <= N; 0, for all other x._ _The first choice X1 obviously follows this distribution. Show that the second choice X2 is distributed according to f(x) too._ My work: I did some research on exchangeable random variables. My thinking is that if I can show the sequence of random variables X1,..., Xn is exchangeable, then it follows that any of the Xi's can be the first random variable in the sequence. Therefore, each Xi will follow f(x). My question: I am unsure if I properly understand what being exchangeable means. Also, I am not sure how you would properly prove that this sequence is exchangeable. I would really appreciate any tips on going about it. Thank you.


title: Does a choice-based conjoint study have to be full profile?

text: In a Conjoint study, it's difficult for respondents to evaluate more than 6 product attributes at a time. Often this is solved via the use of Adaptive Conjoint Analysis (ACA), in which the questionnaire is modified for each individual respondent as the survey is being taken. In ACA, it is not necessary to show the full profile -- i.e. all attributes -- of each product. Partial profiles are shown. A study can include up to 30 attributes, but respondents are never asked to consider more than 5 at a time. However, in order to do ACA, as far as I know of at this time, one must have the survey hosted by a very expensive (e.g. $10,000) conjoint-oriented survey host. My question is, is it possible to do partial profile Conjoint, without using ACA?


title: Estimating minimum sample size when independent variable is categorical

text: Let's say the data is faults in a program. Faults are seeded in a method of the program. Assume that I can find all faults in the program. But, because of cost to find all of them, I want to sample the faults. Problem is how to estimate minimum size of the sample. Independent variable is method names in the program. And dependent variable is the number of faults according to each method. I have no idea to determine the sample size because the independent variable is categorical. One possible idea is to assume that the distribution is normal and estimate the sample size on confidence intervals. But I wonder how to get sample mean and variance for the categorical variable. I'd like to know if there is a better idea or typical solution for this problem.


title: Watermarking data for datamining

text: I'm in a work group that analyzes medical data. Unfortunately there's a lot of distrust if measured data gets to a competitor or is manipulated. So I was wondering if there would be a way to "watermark" the measured data before it leaves the house in a way that would not affect some chosen statistics. Searching for it found mostly commercial solutions for marking audio or video - which isn't applicable for us.


title: How to run a multiple nonlinear regression in Excel or R?

text: I am very inexperienced with R and have only a limited background with Excel but have some data that I need to run a multiple non-linear regression with. What is the best way to do this? With R or Excel? What commands would I use?


title: Linear regression dependence on end points?

text: When playing a bit with a live linear regression visualization at http://www.geogebratube.org/student/m28709 I noticed that point at the ends of the X range have a large influence on the slope of the regression line. Is there an accepted way to give weights depending on the X coordinate to reduce the influence of possible outliers at the border of the X range? Or maybe another method?


title: How to find new standard deviation for new data from old one?

text: I have a calculated standard deviation, mean and number of elements in list. For every new numerical entry to list, I should compute new standard deviation via using old one and new entry data. I have stucked how to do in paper so my code implementation is also. How can I formulate new standard deviation after new entry is done? What I have;   * current mean    * current std   * # of elements in list Pseudocode;               for each new data, x1             return new std deviation


title: Statistical distance tests

text: I am using around 200 variables to compare 2 populations, group A and group B. Here's a representation of what the data looks like               Var1    VarA    VarT    Var3    group     2       0.56    10       152    A     3       0.43    11       187    A     5       0.25    15       149    A     7       0.12    12       132    A     2       0.78    8        132    A     1       0.55    7        178    A     0       0.35    7        240    A     0       0.36    7        205    A     0       0.3     6        137    A               Var1    VarA    VarT    Var3    group     5       0.24    10       205    B     8       0.76    8        120    B     7       0.43    9        531    B     4          0    5        203    B     5       0.52    11       215    B     2       0.63    14       224    B     1       0.46    12       211    B     2       0.64    1        212    B     6       0.55    7        134    B     5       0.78    7        222    B      I was trying to find the variable that shows the biggest difference between these 2 groups, by comparing the data column-by-column (e.g. "VarT" from group A with "VarT" from group B). I am working in R. I found the statistical distance tests (such as Kullback–Leibler divergence, Cramér–von Mises, Bhattacharyya distance and Kolmogorov–Smirnov test) useful in this case. The problem I am having is that I cannot decide which test to use, what are the preferences among these tests? Unfortunately I wasn't able to find any comparisons between these tests. My data is numerical only with no NA values, and it is highly variable. Direction towards other methods for testing the "overlap" between the 2 distributions in each variable are welcome. Many thanks,


title: Multivariate logit

text: I have n individuals with k possible financial products where they can put a % of their incomes. I need to predict for each individual the percentage of income that they allocate in each product. Is there a way to model concurrently all the k percentages, exitimating also the correlation matrix? In another work I used a Multivariate probit in Stata (Multivariate probit regression using simulated maximum likelihood, L.Cappellari 2004): it's the same thing that I'd like to have now, but in that case the Y have to be 0 or 1, now I have Y that is a rate between 0 and 1. Suggestions?


title: Comparing slopes in mixed-effect model

text: My data looks like the attached picture. The dependent variable indirectly measured physical activity. I tried to use mixed-effect model rather than RM Anova, because my actual data is imbalanced. ![enter image description here](http://i.stack.imgur.com/rshMh.png) In this data, my research question is to figure out   1. How treatment effects are different   2. How genotype effects are different   3. How physical activity change over time.  What I really want to do is to compare each slope of the physical activity change within day1, day2, and day3 (like..animal a: (4.5-2.5)/4 vs (4.5-2.9)/4 vs (4.8-2.6)/4, not exactly in this way but.. that's what I want to compare). I am currently using R so I wrote a code using lmer function. I formulated Day and Hour variables as continuous variable because I want to see the linear trends. I didn't assume the cross level interaction (between-within var interaction) in the model.               lmer(DV ~ 1 + treatment*genotype + Day + Hour + Day*Hour + (1 | ID), data)      I am assuming that testing Day*Hour interaction term may be one way to look at it, but I am not sure.. Could you give me any advice about it or give me another way to look at it?? Thank you so much for your time in advance.


title: Help me understand $p$-values in Bayesian glm

text: I am trying to run a Bayesian logit on the data here. I am using `bayesglm()` in the `arm` package in R. The coding is straightforward enough:               df = read.csv("http://dl.dropbox.com/u/1791181/bayesglm.csv", header=T)     library(arm)     model = bayesglm(PASS ~ SEX + HIGH, family=binomial(link="logit"), data=df)      `summary(model)` gives the following output:               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  0.10381    0.10240   1.014    0.311         SEXMale      0.02408    0.09363   0.257    0.797         HIGH        -0.27503    0.03562  -7.721 1.15e-14 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 2658.2  on 1999  degrees of freedom     Residual deviance: 2594.3  on 2000  degrees of freedom     AIC: 2600.3      Please walk me through this. I understand that this code uses a very weak prior (since I am not specifying the prior means) so the output is going to be practically the same if I used `glm()` instead of `bayesglm()`. But the output should still be in the Bayesian spirit, right? What are the $p$-values and $z$-values here? Aren't these frequentist inference tools? Are they interpreted differently here?


title: Can you compare regression slopes from a mixed model?

text: I know that we can do a t-test to compare the slopes of 2 independent regressions line using the coefficients and standard error (example is here: https://stat.ethz.ch/pipermail/r-help/2009-April/195079.html) Does this same approach work for mixed models?


title: What exactly is the "wid" argument in the ezANOVA function in the R package "ez"?

text: I'm very new to all this, and I am testing different ways to perform a two-way type III ANOVA on my data.   * I have tried `anova()` from the `stats` package, after fitting a linear regression with `lm()`;   * I have tried `Anova()` from the `car` package, using the same linear regression (and this gives me the same result as `anova()` when I use `type="II"` \- I thought `anova()` used type I SS by default?).   * And I am now trying to use `ezANOVA()` from the `ez` package. With this last one, I can't understand what the `wid=.()` argument is (even reading the help), and as it is not optional, I can't leave it blank. What I am trying to use is as follows, with its result:               > attach(data)     > library("ez")     > ezANOVA(data=data, dv=.(AG.DW), wid=.(), within=.(Genotype, Treatment), type=3)     Warning: Converting "" to factor for ANOVA.     Error in sort.list(y) : 'x' must be atomic for 'sort.list'     Have you called 'sort' on a list?      Is this the right script? What is `wid` and what should I fill it with? Concerning my data, the columns `Genotype` and `Treatment` are my two factors, and I want to see if there is an interaction when looking at the aboveground dry weight of my plants (column `AG.DW`). My data is balanced. I am sorry if information is missing or inaccurate here: this is my first contribution here, and I am only discovering statistics at the moment (and I can't see how to join my data file).


title: Comparing two experimental groups and a treatment group using ANOVA

text: I performed a two-way ANOVA. I have one control group and 2 experimental groups. One experimental group received one strategy, the other group received another kind of strategy. Now, I want to know which strategy had a better effect on the group. Which table in the two-way ANOVA will help me understand which strategy had a better effect? (I used pre-test and post-test)


title: bag of words in an online configuration, for classification / clustering

text: I have a set of image documents. I extract text keywords from this images using OCR to represent each image as a bag of words (a vector where each value is the number of occurrence of a word in the document). Then I can apply a classification or clustering algorithm on the obtained dataset. However, this vector representation as bag of words is possible only if I have the entire set of documents (to be able to have the entire vocabulary, i.e. all words). How can I do that (i.e. extract my bag of words) if I'm in an online configuration (using an online clustering) where the documents are available one by one (for a data stream) and each document should be processed as soon as it is available ?


title: Stationary Process in Plain English

text: How would you describe stationary process in plain English to someone with no mathematical background, using real life examples? The target audience is adults with reasonable intelligence, but most have been out of school for years, if not decades.


title: Calculate extreme tails of distribution

text: This is how my empirical distribution look like: I have calculated a genetics index of fixation in populations, called FST. My values have the maximum value of one. This gives me an empirical distribution of FST values for the markers on each chromosome. Now, I have a list of markers that I would want to know whether they fall in the 5 percent tail of my empirical distribution or not. This sounds very simplistic but I thought as I have all the data values, I can compute the 5 percent by saying I have 5 out of 100 so how many out of the number of my sorted FST values. Do you think this is a correct approach or I am being too simple-minded? Thank you in advance, Best, Lucia


title: Sheather-Jones bandwidth algorithm implementation in java?

text: The Sheather-Jones method for selecting an appropriate bandwidth for kernel density estimation generally produces better results than simpler methods such as Silverman's rule of thumb and Scott's variant of this rule of thumb, both of which give bandwidths which are generally too large. I use the Sheather- Jones method frequently in R but I'm looking for an implementation in Java. Is anyone aware of such an implementation? I have been unable to find one online. The paper describing the method is here.


title: Minimum sample size

text: Probably a basic question, but here goes:   * Suppose we have 1000 people who respond to a survey (500 men, 500 women)   * Respondents can choose to "like" a series of pages.   * Each respondent can "like" as many pages as they want; it is possible to never like anything.   * For page 1, 200 men like it and 150 women like it.   * For page 2, 300 women like it and 50 men like it.   * For page 3, 4 women like it and 1 man likes it. My general question is: for small response scenarios (i.e., page 3), what is the minimum sample size before I can really make a conclusion about the male/female general interest? I can compute a z-score, but I am not sure if it takes this "sub population" into account (I may be confused about this). I know that my definition of significance will depend on what error rate and confidence I can tolerate, but let us suppose it's the standard confidence interval (95%)


title: Questions on paramatric and non-parametric bootstrap

text: I am reading the chapter on Frequent Statistics from Kevin Murphy's book " _Machine Learning - A Probabilistic Perspective_ ". The section on bootstrap reads: > The bootstrap is a simple Monte Carlo technique to approximate the sampling > distribution. This is particularly useful in cases where the estimator is a > complex function of the true parameters. > > The idea is simple. If we knew the true parameters $θ^∗$ , we could generate > many (say $S$) fake datasets, each of size $N$, from the true distribution, > $x_i^s \sim p (·| θ^∗ )$, for $s = 1 : S, i = 1 : N$. We could then compute > our estimator from each sample, $\hat{\theta^s}=f (x^s_{1:N})$ and use the > empirical distribution of the resulting samples as our estimate of the > sampling distribution. Since $\theta$ is unknown, the idea of the > **parametric bootstrap** is to generate the samples using $\hat{\theta}(D)$ > instead. > > An alternative, called the **non-parametric bootstrap** , is to sample the > $x^s_i$ (with replacement) from the original data $D$ , and then compute the > induced distribution as before. Some methods for speeding up the bootstrap > when applied to massive data sets are discussed in (Kleiner et al. 2011).   * **1**. The text says:  > If we knew the true parameters $\theta^*$ ... we could compute our estimator > from each sample, $\hat{\theta^s}$...         but why would I use the estimator of each sample if I **already** know the true parameters $\theta^*$?    * **2**. Also, what is the difference here between the empirical distribution and the sampling distribution?   * **3**. Finally, I don't quite understand the difference between **parametric** and **non-parametric** bootstrap from this text. They both infer $\theta$ from the set of observations $D$, but what is exactly the difference?


title: Fast procedurally built inverse model

text: Suppose I have a DB into which the following key values are entered (in that "random" sequence)               {16, 32, 256, 2, 8, 64, 4, 1, 128, 512}      Assuming the values get sorted on insertion ($B^x$ tree for instance), what's the fastest way to maintain a **running model** of the data? ### What I want: A running model would be a function that delivers record at which a certain element would be found. Since the DB $A$ would look something like               record:  |000|001|002|003|004|005|006|007|008|009|     value:   |  1|  2|  4|  8| 16| 32| 64|128|256|512|      since the data can be plotted optimally with $2^x$, a model would be $\frac{\log(x)}{\log(2)}$. If we wanted to look for 334, we would expect it to be at $\frac{\log(334)}{\log(2)}=8.38$. ### What I have so far: The best solution I have so far is to maintain a smaller superDB $S$ that keeps every $n$th element and linearly interpolates between them. assume that $n=3$ then $S$ would look like               record:  |000|001|002|003|     position:|000|003|006|009|     value:   |  1|  8| 64|512|      then if we wanted to find the position of $i$=334 we'd first seek in $S$ until $ \mathbb{value}> i$, and mark the record position $k$. then the expected position would be $$(i-[(k-1)\rightarrow\mathbb{value}])* \frac{[(k)\rightarrow\mathbb{value}]-[(k-1)\rightarrow\mathbb{value}]}{[(k)\rightarrow\mathbb{position}]-[(k-1)\rightarrow\mathbb{position}]} + (k-1)\rightarrow \mathbb{position} $$ so k=3. (334-65)*((9-6)/(64-512))+6 = 7.80133928571, but this is clearly wrong because record `008` is obviously greater than 334, so it can't be below that.   * $+$ computationally nonintensive (1 multiplication, 1 division, 4 dereferentials)   * $+$ error is never greater than $n$.    * $-$ error can be almost $n$ in skewed data.   * $-$ requires pre-seeking   * $-$ worst case scenario gives no gains. An option might be to insert new reference values into the SuperDB based on error, and to remove some when the error is slim, if for instance if $1-\frac{f(k, k-1)}{f(k+1, k-1)}< \mathbb{allowed \, error}$, then record $k$ could be removed from $S$, and a record inserted if $\frac{\mathbb{seeks \, from}( f(k, k-1) )}{n}> \mathbb{allowed \, error}$ to build an increasingly better model of the system. In combination with the above, parametric interpolation (bezier, quadratic, circular, etc...) could perhaps be used, which however come with the drawbacks of increased overhead. What's a better way to do this?


title: samples from forecasts of VAR time series model in R

text: I'm trying to do a power analysis for a future experiment with time series financial data. We'll be splitting the data by random (actually, stratified) geographies, so we have a control and experimental group. As part of this process, I'm building a time-series model of historical data in the two groups. They're highly correlated, as you'd hope and expect, and so I'm using a VAR model from the `vars` package. That seems to work fine. What I'd like to be able to do is to pull correlated samples from the predictions of this model. If I'm modeling univariate time series with the `forecast` package, I can use `simulate` to pull random samples from the forecasted series, which works great. What's a good way to generate samples from the model's forecasts that incorporates the covariance of the time series? The `vars::predict` method only gives confidence intervals for each series independently, and presumably I don't want to independently draw from those distributions, or I'd lose the additional power from the correlation. Is there a relatively easy way to do this?


title: How to call glm when response variable is categorical in R?

text: I have following data stored in a file. I am applying 'glm' in R to find linear regression equation to best predict the 'output'.               > tmpData        logOfOutput randomSample multiplied part1 part2 randNormalMean100Std20 output     1    0.0000000           33         11     1    19               89.65387      1     2    0.6931472           76         24     2    18              128.23471      2     3    1.0986123           12         39     3    17              103.70930      3     4    1.3862944           68         56     4    16               99.12617      4     5    1.6094379           50         75     5    15               95.68173      5     6    1.7917595            7         96     6    14              129.27551      6     7    1.9459101           70        119     7    13              104.59333      7     8    2.0794415           55        144     8    12              102.15247      8     9    2.1972246           20        171     9    11               72.43795      9     10   2.3025851           24        200    10    10               80.63634     10     11   2.3978953           32        231     9    11              105.03423     11     12   2.4849067           97        264     8    12               78.10613     12     13   2.5649494           28        299     7    13              107.95286     13     14   2.6390573           99        336     6    14               80.07396     14     15   2.7080502           66        375     5    15              102.01156     15     16   2.7725887           95        416     4    16              119.07361     16     17   2.8332133           42        459     3    17               64.19354     17     18   2.8903718           53        504     2    18              106.23402     18     19   2.9444390           85        551     1    19              151.07976     19     20   2.9957323           48        600     0    20               82.78324     20      I am using the following code to perform the same               fn = "delnowSample.txt"     tmpData = read.table(fn, header = TRUE,  sep= "\t" , blank.lines.skip = TRUE)     cnames = colnames(tmpData)     (fmla <- as.formula(paste(cnames[length(cnames)], " ~ ", paste(cnames[1:(length(cnames)-1)],collapse= "+")))  )     model <- try(glm(formula = fmla, family=binomial(), na.action=na.omit, data=tmpData));     summary(model)       The output that I get is as follow:               > summary(model)          Call:     glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse = "+"))),          family = gaussian(), na.action = na.omit)          Deviance Residuals:           Min        1Q    Median        3Q       Max       -0.37926  -0.11242  -0.03441   0.16087   0.28200            Coefficients: (1 not defined because of singularities)                                                 Estimate Std. Error t value Pr(>|t|)         (Intercept)                                0.2638036  0.3078536   0.857  0.40592         unlist(tmpData["logOfOutput"])             0.9202273  0.2727884   3.373  0.00455 **      unlist(tmpData["randomSample"])            0.0026201  0.0018177   1.441  0.17145         unlist(tmpData["multiplied"])              0.0288073  0.0012359  23.308 1.34e-12 ***     unlist(tmpData["part1"])                   0.2106002  0.0403442   5.220  0.00013 ***     unlist(tmpData["part2"])                          NA         NA      NA       NA         unlist(tmpData["randNormalMean100Std20"]) -0.0006214  0.0024922  -0.249  0.80673         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for gaussian family taken to be 0.04403284)              Null deviance: 665.00000  on 19  degrees of freedom     Residual deviance:   0.61646  on 14  degrees of freedom     AIC: 1.1676          Number of Fisher Scoring iterations: 2      To a large extent it is predicting the Pr(z) correctly as we can see the probabilities of random variable are not significant. The R-square is also high (1-residual.deviance/null.deviance), close to 1. Question 1: In the above data 'part1+part2' is equal to output variable. Is 'glm' not able to identify such type of relations? Question 2: Why the degree of freedom of null and residual deviance are different? Question 3: I need to convert the output variable into categorical variable (i.e. Everything <=10 is 'no' and more than this is 'yes'). What is the best way to call 'glm', when the response variable is 'categorical'. I tried converting 'no' to '0' and 'yes' to 1, and called glm as follows:                model <- try(glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse= "+"))), family=binomial(), na.action=na.omit));      I am getting warning message with this code. Also, I am not sure if this is the correct way to call categorical variable. Edit: I have the following categorical data:               > tmpData                x1  x2 x3  y1     1  0.16294456   1  1  no     2  0.80494934   2  2  no     3  0.28962222   1  3  no     4  0.07177347   2  4  no     5  0.54830544   1  5  no     6  0.67655327   2  6  no     7  0.45189608   1  7  no     8  0.82412502   2  8  no     9  0.09076793   1  9  no     10 0.12221227   2 10  no     11 0.56751754 111 11 yes     12 0.04970992 222 12 yes     13 0.56162037 111 13 yes     14 0.96617891 222 14 yes     15 0.50994534 112 15 yes     16 0.70093692 212 16 yes     17 0.02034940 212 17 yes     18 0.78356903 121 18 yes     19 0.58439662 213 19 yes     20 0.31729282 212 20 yes      And the following code:                 fn = "delnowSample.txt"       tmpData = read.table(fn, header = TRUE,  sep= "\t" , blank.lines.skip = TRUE)       tmpData       model <- glm(formula = 'y1~x1+x2+x3', family=binomial(), na.action=na.omit, data=tmpData)       summary(model)       This one doesn't seem to be working??


title: Specific type of dummy regression

text: I would like to regress the values on the y-axis by the values of the x-axis. As you can see, the relationship is not linear. Values of y are only positive can not exceed 4, values of x are always positive. Is it possible (and do you think it is a good idea) to do the following:   * Use dummy variables for bins of x-values (bins of 0.5 for example)   * As an extra restriction, the value of the dummy variable coefficient in bin i can not be lower than the value of the dummy variable of bin i-1 That way I hope to get a nice smooth monotonic rising regression function through the data. (I want to avoid the decrease of the at the end of the X range of the LOESS fit as it is not meaningful) Do you think this is a good idea and is this possible in R? Thank you beforehand. ![enter image description here](http://i.stack.imgur.com/w1FjH.jpg)


title: Constructing the approximate confidence set for parameter vector beta in least square regression

text: From Elements of Statistical Learning, on page 49 (in the context of least square regression), we are given the approximate confidence set (shown in image) for the parameter vector $\beta$. ![enter image description here](http://i.stack.imgur.com/Xx4Ak.png) I want to implement this in R. However, I have no idea how to make $\beta$ as the subject,i.e. I could not construct the confidence set in the following form $$\beta \leq \hat{\sigma}^2 \chi^2_{p+1} \times something$$ I don't really know how to rearrange the terms (of the confidence set from the image). Can someone kindly help me out here please ?


title: What sample size do I need?

text: I have a population of 709 and the ability to draw a truly random stratified sample. I understand how to do the stratified part once I determine the sample size. One calculation indicated that I need 196 items for 95% confidence. Given that I am working with content analysis, that is a lot of records to review. Can someone please confirm the calculation that 196 items ($n=196$) are needed from $N=709$ for 95% confidence? Is there another way (other than altering the confidence interval) to reduce the required sample size?


title: Machine Learning with Skewed Classes in R

text: I am looking for some suggestions on what methods are appropriate for training a dataset with a high skew in the outcome classes. The ratio of Class 0: Class 1 is about 20:1 and I am looking to maximize the accuracy for identifying Class 1 outcomes. This is similar to oft discussed topics such as cancer detection. I have used some methods before but am trying to find if there is any comprehensive resource / suggestions that talks to the different methods for these cases. Examples of how they are applied in R (packages, etc) or with caret would be useful. It is a sparse dataset with about 100K examples of which 5000 belong to Class 1 and the rest to Class 0. Each example has about 20 features, and includes null values, Thanks.


title: Kernel density estimation incorporating uncertainties

text: When visualising one-dimensional data it's common to use the Kernel Density Estimation technique to account for improperly chosen bin widths. When my one-dimensional dataset has measurement uncertainties, is there a standard way to incorporate this information? For example (and forgive me if my understanding is naïve) KDE convolves a Gaussian profile with the delta functions of the observations. This Gaussian kernel is shared between each location, but the Gaussian $\sigma$ parameter could be varied to match the measurement uncertainties. Is there a standard way of performing this? I am hoping to reflect uncertain values with wide kernels. I've implemented this simply in Python, but I do not know of a standard method or function to perform this. Are there any problems in this technique? I do note that it gives some strange looking graphs! For example ![KDE comparison](http://i.stack.imgur.com/kuaHG.png) In this case the low values have larger uncertainties so tend to provide wide flat kernels, whereas the KDE over-weights the low (and uncertain) values.


title: Clarification on interpreting Wald's test and Likelihood ratio tests

text: I am running multinomial logistic regression analysis on my data. The response variable is the number of calves produced each year (0,1, or 2). I am trying to evaluate the influence of the _X_ variables on the odds of producing a calf. My _X_ variables are predation risk (WR; continuous), age of mother (age; categorical or continuous), time (wolf; categorical). First, I have 4 different age classification schemes (i only show 2) - I want to know which one of the age of mother would be "best" to use. I could use it as continuous variable - or as categories based on biological reasoning for senescence in older moose (old ladies don't invest in reproduction as much). So, I thought I would use a likelihood ratio test.               library(mlogit)          modata.model1 <- mlogit(no.C ~ 1 | 1, data=modata, reflevel="1", na.action = na.omit)      modata.model2 <- mlogit(no.C ~ 1 | age, data=modata, reflevel="1", na.action = na.omit)      modata.model3 <- mlogit(no.C ~ 1 | age2, data=modata, reflevel="1", na.action = na.omit)        * * *                lrtest(modata.model2,modata.model3)          Likelihood ratio test          Model 1: no.C ~ 1 | age     Model 2: no.C ~ 1 | age2       #Df  LogLik Df  Chisq Pr(>Chisq)         1   4 -213.22                              2   4 -207.57  0 11.309  < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      QUESTION: to interpret this output - there was a significant difference in the loglikelihood when we comparing the continuous age to a categorical age with 2 classes. The loglik is smaller for model 1 and therefore it would be better to use? Or do I have that backwards? Next I was going to use the Walds test to evaluate nested models. To see if the addition of a variable was worth it.               modata.model7 <- mlogit(no.C ~ 1 | age+WR, data=modata, reflevel="1", na.action = na.omit)           modata.model8 <- mlogit(no.C ~ 1 | age+WR+wolf, data=modata, reflevel="1", na.action = na.omit)       * * *               Wald test          Model 1: no.C ~ 1 | age + WR     Model 2: no.C ~ 1 | age + WR + wolf       Res.Df Df  Chisq Pr(>Chisq)     1    244                          2    242  2 0.5828     0.7472      QUESTION: this tells me that there is no significant improvement when there is an additional variable of wolf added?? So, then I can use the smaller model or do I use the one with the smaller Res.DF? In addition to confirming my interpretations of the results I have 2 side questions... 1)to get the null model for `mlogit` library - is my `modata.model1` correct? I want the intercept only model to compare against. 2) Hosmer and Lemshow suggest by getting Wald values to get significance levels for each coefficient - in mlogit, thats the same as using `summary(model)` and there they provide the t-values with p instead of needing to do an additional Walds test? (NOTE in the below model i use 2 category age class instead of continuous)               summary(modata.model8)          Call:     mlogit(formula = no.C ~ 1 | age2 + WR + wolf, data = modata,          na.action = na.omit, reflevel = "1", method = "nr", print.level = 0)          Frequencies of alternatives:         1     0     2      0.652 0.244 0.104           nr method     6 iterations, 0h:0m:0s      g'(-H)^-1g = 7.4E-06      successive function values within tolerance limits           Coefficients :                   Estimate Std. Error t-value Pr(>|t|)        0:(intercept)  0.38730    0.40144  0.9648 0.334657        2:(intercept) -2.40317    1.04546 -2.2987 0.021523 *      0:age21       -1.39607    0.43989 -3.1737 0.001505 **     2:age21        0.53952    1.07275  0.5029 0.615012        0:WR          -1.46584    0.64797 -2.2622 0.023686 *      2:WR          -0.19214    0.60856 -0.3157 0.752206        0:wolf1        0.56055    0.63292  0.8857 0.375797        2:wolf1        0.42642    0.72375  0.5892 0.555744        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Log-Likelihood: -203.11     McFadden R^2:  0.053606      Likelihood ratio test : chisq = 23.009 (p.value = 0.00079359)


title: Fitting a Mixed Model with Random and Repeated effects in SAS

text: I have want to fit a linear regression with repeated measures and random effects. The data come from clinical observations. In CT images The dependent variable is the diameter of a lymph node lesion which is determined by a CT- image. Each picture is rated by two readers. Each lymph node is observed in a follow-up up to 7 times. Within one patient there could be several lesion. So I have: Repeated measures: 2 raters at each timepoint, up to 7 timepoints Random effect: Several lymph nodes within on patient. The aim is to estimate the difference between the two raters and the change of the diameter in time. How should I write the model equation and how can I fit the model in SAS? Thanks for any help. Variables: LesionID PatientID Reader time diameter               PROC MIXED DATA=Dataname EMPIRICAL;       CLASS LesionID Zeitpunkt Patid Rater;       MODEL DIAMETER= Zeitpunkt Rater ;       RANDOM  ???;       REPEATED ???     RUN;


title: How to normalize bimodal (or multimodal) distributions?

text: If I have multiple data series, a = [a1, a2, ... a100] ~ bimodal with mu_a1, mu_a2, sigma_a1, sigma_a2, b = [b1, b2, ... b100] ~ bimodal with mu_b1, mu_b2, sigma_b1, sigma_b2, c = [c1, c2, ... c100] ~ bimodal with mu_c1, mu_c2, sigma_c1, sigma_c2, how do I normalize a, b, and c to follow the same distribution? Say, for the purposes of combining the data? If a, b, and c followed a Gaussian distribution, for example, all I'd have to do is subtract the means and scale to the same standard deviation. But what about in bimodal or multimodal distributions with multiple means and standard deviations? Thanks in advance.


title: which regression model: nested, censored data

text: I have the following data structures and would like to know your opinion which regression models are most suitable (and available in R) model 1: \- 20,000 cases \- nested in 500 spatial units (only intercept variations are needed) \- the response variable is censored (only values >1) are possible \- 10 independent variables (only linear main effects are needed) model 2: \- 20,000 cases \- nested in 500 spatial units (only intercept variations are needed) \- the response variable is categorized in 3 classes where class 1 < class 2 < class 3 \- 10 independent variables (only linear main effects are needed)


title: How to compare difference between two time series?

text: I am working on my thesis where I'm examining how strong emotion people show to different events. My problem is (1) that I have VERY little experience with statistics and math, so I'm kind of lost with all different methods and would be really happy if a 'easy' answer could be given (without lots of algorithms and stuff). So the main question: In the graphs below I want to check if the yellow line is significant different in the second graph compared to the first, and also how different. How is this done in the easiest way? ![graphs](http://i.imgur.com/y8Y5lMU.jpg) I have searched for an answer for two days, but not being able to find something that could help me. So I hope that someone can help me on the right way!


title: Forecasting agricultural commodity prices with R

text: I would like to create a predictive model in order to forecast the price of an agricultural raw material. I got time series for the prices and the production of this raw material, and also for the prices of oil. I think I'm going to model it with a VECM, but I'm just starting to use R. 1) How do I use R to create such a model ? I started to have a look at the vars package but I'm not sure how to use it... 2) The data I have for the production are given annualy while the one for the prices are given monthly. How can I use both types of data, given the fact that I want to have the more information possible (and thus keeping monthly data) ? Thanks in advance


title: Time Series Analysis Forecasting with Asset Bubbles

text: I'm working on a final project for my class and I'm doing it on an asset's price over time. The scientific question I want to ask is "Is it possible to reasonably forecast the price of asset x with small uncertainty to predict an upcoming crash using a previous crash of the same asset?" Is this feasible? This is for my college time series class. I'll be using R.


title: F-test for Lack-of-Fit in SPSS

text: Some googling revealed that doing the F-test for Lack-of-Fit in SPSS is not so trivial. It seems one has to “trick” SPSS to do that. See for example this. Can anybody suggest a better source of information on how this can be done? I have SPSS 16. Of course I know it can be easily done using R but I am interested in the SPSS way. Thanks


title: What is the effect of having correlated predictors in a multiple regression model?

text: I learned in my linear models class that if two predictors are correlated and both are included in a model, one will be insignificant. For example, assume the size of a house and the number of bedrooms are correlated. When predicting the cost of a house using these two predictors, one of them can be dropped because they are both providing a lot of the same information. Intuitively, this makes sense, but I have a some more technical questions:   1. How does this effect manifest itself in p-values of the regression coefficients when including only one or including both predictors in the model?   2. How does the variance of the regression coefficients get affected by including both predictors in the model or just having one?   3. How do I know which predictor the model will choose to be less significant?   4. How does including only one or including both predictors change the value/variance of my forecasted cost?


title: Randomized SVD and singular values

text: Randomized SVD decomposes a matrix by extracting the first k singular values/vectors using k+p random projections. This works surprisingly well for large matrices. My question concerns the singular values that are output from the algorithm. Why aren't the values equal to the first k-singular values if you do the full SVD? Below I have a simple implementation in R. Any suggestions on improving the performance would be appreciated.                rsvd = function(A, k=10, p=5){            n = nrow(A)            y = A %*% matrix(rnorm(n * (k+p)), nrow=n)            q = qr.Q(qr(y))            b = t(q) %*% A            svd = svd(b)            list(u=q %*% svd$u, d=svd$d, v=svd$v)         }              > set.seed(10)              > A <- matrix(rnorm(500*500),500,500)              > svd(A)$d[1:15]          [1] 44.94307 44.48235 43.78984 43.44626 43.27146 43.15066 42.79720 42.54440 42.27439 42.21873 41.79763 41.51349 41.48338 41.35024 41.18068              > rsvd.o(A,10,5)$d          [1] 34.83741 33.83411 33.09522 32.65761 32.34326 31.80868 31.38253 30.96395 30.79063 30.34387 30.04538 29.56061 29.24128 29.12612 27.61804              B <- matrix(rnorm(500*50),500,500)  # rank 50          > rsvd(B,10,5)$d      [1] 86.48035 83.02114 81.03988 80.04358 77.24979 76.10945 74.47357 74.08382      [9] 72.85898 72.06897 69.59526 67.70750 66.53867 62.96446 61.50838          > svd(B)$d[1:15]      [1] 92.44779 91.47689 88.71948 88.08170 87.24533 85.13312 84.14741 83.71757      [9] 82.80832 81.43005 80.73903 79.92959 78.87421 78.33509 77.38431      As Joris pointed out, I have this posted on stackoverflow as well. you can find the revelant conversation here http://stackoverflow.com/questions/4224031/randomized-svd-singular-values Also see the relevant paper by Martinsson et al: A randomized algorithm for the decomposition of matrices


title: MLE on a restricted parameter space

text: Suppose $X_1, \ldots, X_n$ are a random sample from a normal distribution with mean $\theta$ and variance $1$. Find the maximum likelihood estimator of $\theta$, under the restriction that $\theta \geq 0$. I have found the MLE when there is no restriction on $\theta$. $$L(\theta|\mathbf{x})=(2\pi)^{-n/2} \text{exp}\Big(-\frac{1}{2}\sum_{i=1}^{n}(x_i-\theta)^2 \Big)$$ Differentiating $L(\theta|\mathbf{x})$ with respect to theta, and setting it equal to zero yields $$\sum_{i=1}^{n}(x_i-\theta)=0 \implies \theta=\frac{1}{n}\sum_{i=1}^{n}x_i =\bar{x}$$ The second derivative is negative, hence the MLE of $\theta$ is $\hat{\theta}=\bar{X}$. How would I incorporate the restriction $\theta \geq 0$.


title: Amoeba Interview Question

text: I was asked this question during an interview for a trading position with a proprietary trading firm. I would very much like to know the answer to this question and the intuition behind it. Thank you! **Amoeba Question: A population of amoebas starts with 1. After 1 period that amoeba can divide into 1, 2, 3, or 0 (it can die) with equal probability. What is the probability that the entire population dies out eventually?**


title: Probability in mafia/werewolf party game

text: Okay, first post on this Stack Exchange board, and I'm no mathematician so I kind of lack the proper jargon. But here goes. I'm currently playing a game known as werewolves, which is based on the better-known mafia party game, on a forum. There's a number of players (11 in our case), among which a number of werewolves (let's assume 3; actual number unknown but can be a variable in the calculations). Among the non-werewolves there could be a witch which has a life potion, and a guardian angel who can choose someone (but not herself) to protect during a night. At night, the werewolves choose someone to attack. If that person was protected by the guardian angel, they fail. If not, then the witch can save the dying victim with her potion. After the first night, we observed no-one had died. This means either the angel had made a lucky choice, or the witch saved the victim. Note that the angel doesn't have any information about anyone and might accidentally protect a werewolf (who wasn't in any danger). Nor do the werewolves know anything about anyone's roles but their comrades'. It's easy enough to calculate the odds that the angel would make a lucky choice. It's the chance that she chooses a non-wolf multiplied by the chance that the attacked person was chosen. The chance she chooses a non-wolf is (11 - 1 - 3)/10. The chance that the person protected happens to be the one that gets attacked is 1/8. So this ends up with 7/80 or 0.0875. I ran a little program that simulated this to verify this outcome. Now we're trying to make out what the odds are that the result after the first night, namely no-one died, are due to the guardian angel. I believe that it is simply the chance that the angel happened to protect the right person. But another player in the discussion believes the chance to be much greater, because we observed that no-one died and thus the scenario of a successful murder is excluded. According to me that is falsely using the outcome as a condition. Such as stating that a coin flip had a 100% chance of being heads after doing the flip and seeing it land as heads. However, this is somewhat more complex and I might be missing something. What is the correct way here to determine the odds that given the current situation (no-one died last night) this was brought on due to a lucky choice of the angel? EDIT: added some more info. The witch is told at night who was attacked if that person did not happen to be the angel's protégé. It is assumed that if the witch is attacked and not protected, she will always use the life potion on herself. Otherwise, the choice to use the potion is hers. Werewolves can't target each other. They know each other's identities and can communicate. They collectively agree upon a single target.


title: Reference in Genetic Algorithm Neural network

text: I am planning to work on my thesis with genetic algorithm neural network (GANN), since that is a modern statistical technique instead of by using discriminant analysis which is a traditional statistical technique. My idea is by using this GANN, this method could be used to classify the disease in human which is in categorical. Anybody can give me some literature/link/references actually i can download it also, so i can understand how GANN is work, I am still a basic for these, so I need it completely include the literature for the formula. Thanks!


title: Uncertainty analysis

text: Here is my situation. I am trying to predict the 'entire' distribution of the dependent variable, not just the mean( or conditional mean). Does it then make sense to seprateley predict quantiles of this variable to learn about the new predicted CDF? I intend to use this CDF as one input (say sample from it) in a model that also has other random variables. If my approach isnt advisable, please recommend me an alternative.


title: Multiple testing correction when confirming a correlation with an extended dataset

text: Consider a given dataset on which I performed a set of independent tests on different (non-independent) variables, and obtained several p-values. I estimated the number of independent tests as being $n$ and corrected the p-values obtained using Sidak or Bonferroni correction. Now let us say I obtained for one of these a marginaly significant result and acquired an extended dataset to conclude if the correlation is significant or not. What is tne number of independent tests I should use for the correction on this larger dataset, $n$ or $n+1$?


title: Discrete uniform random variable(?) taking all rational values in a closed interval

text: I just had an (intellectual) panic attack.   * A continuous random variable that follows a uniform in a closed interval $U(a,b)$: a comfortably familiar statistical concept.    * A continuous uniform r.v. having support over the extended reals (half or whole): not an r.v. proper, but a basic Bayesian concept for an improper prior, useful and applicable.   * A discrete uniform taking a finite number of values: let's throw a geodesic dome, no big deal.  But what about a function that has as its domain all the rationals that are included in a closed interval with integer bounds (start with the $[0,1]$ if you wish)? And we want to use it in a probabilistic framework, requiring that each possible value has equal probability with all the others? The number of possible values are countably infinite (which characterizes a lot of discrete distributions), but how to express the probability of a single value given that we want probabilities equal? Can we say-show-prove that such an entity is (is not) a random variable? If not, is this another incarnation (perhaps already well-known) of an "improper prior"? Is it possible that this entity is in some well-defined sense, however special, "equivalent" to a continuous uniform r.v.? Or I just committed a cardinal(ity) sin? It appears that the fact that the domain is a closed interval does not let me let go. Bounded things are usually manageable. The questions are many in order to be indicative of the internal maelstrom- I am not asking to get answers to each one of them. At any time that I may come up with any insights, I will update. **UPDATE:** the present question just acquired a constructivist sequel here.


title: Data requirements for building an HMM

text: I'm trying to describe the requirements that should be placed on data when building a Hidden Markov Model. Should the data fit the requirements for a Markov chain? That is, must it have discrete, sequential states, a stationary probability distribution of each state, and IID properties? What exactly are the requirements for a Markov chain?


title: Markov-Switching model of distribution

text: I was wondering if anyone knows of any work concerning markov-switching models and more precisely that changes the distribution of the error terms? I've found some literature for conditional mean, variance, and degrees of freedom but i can't find anything for different distributions. thanks.


title: Why doesn't Wilks' 1938 proof work for misspecified models?

text: In the famous 1938 paper ("The large-sample distribution of the likelihood ratio for testing composite hypotheses", Annals of Mathematical Statistics, 9:60-62), Samuel Wilks derived the asymptotic distribution of $2 \times LLR$ (log likelihood ratio) for nested hypotheses, under the assumption that the larger hypothesis is correctly specified. The limiting distribution is $\chi^2$ (chi-squared) with $h-m$ degrees of freedom, where $h$ is the number of parameters in the larger hypothesis and $m$ is the number of free parameters in the nested hypothesis. However, it is supposedly well-known that this result does not hold when the hypotheses are misspecified (i.e., when the larger hypothesis is not the true distribution for the sampled data). Can anyone explain why? It seems to me that Wilks' proof should still work with minor modifications. It relies on the asymptotic normality of the maximum likelihood estimate (MLE), which still holds with misspecified models. The only difference is the covariance matrix of the limiting multivariate normal: for correctly specified models, we can approximate the covariance matrix with the inverse Fisher information matrix $J^{-1}$, with misspecification, we can use the sandwich estimate of the covariance matrix ($J^{-1} K J^{-1}$). The latter reduces to the inverse of the Fisher information matrix when the model is correctly specified (since $J = K$). AFAICT, Wilks proof doesn't care where the estimate of the covariance matrix comes from, as long as we have an invertible asymptotic covariance matrix of the multivariate normal for the MLEs ($c^{-1}$ in the Wilks paper).


title: How to calculate the expected zeros in a Poisson distribution?

text: I am modelling the nights spent at hotels (count data) fitting a few predictors in the model. I'd like to know how to calculate the expected number of zeros in this distribution, as I suspect that I should move to a zero-inflated Poisson distribution.


title: Mean pairwise Pearson correlation

text: Consider a set of $M$ multidimensional observations $X=\\{x_i\\}_{i=1..M}$ for $x_i\in\mathbb{R}^N$, and define $c_{ij}$ to be the Pearson correlation coefficient between $x_i$ and $x_j$. Is there a simple expression or useful approximation for the mean pairwise Pearson correlation $C=\frac{1}{M^2}\sum_{ij}c_{ij}$? Alternatively, is there a simple expression for other measure which is similar in spirit to the "mean pairwise Pearson correlation", perhaps defined differently than $C$ above? (This is motivated by the fact that if $c_{ij}$ was square Euclidean distance, then we had $C=2Var(X)$ and thus I wonder if there is a similar result for Pearson correlaion.) Thanks!


title: inequality in bivariate normal variable

text: Let $U_1=(X_1,Y_1)^T,\dots,U=(X_n,Y_n)^T$ are i.i.d. copies of $U=(X,Y)^T\sim N_2(0,\Sigma)$ where $$ \Sigma= \begin{pmatrix} \sigma^2 & \rho\sigma\tau \\\ \rho\sigma\tau & \tau^2 \end{pmatrix} $$ such that \begin{equation} \sum_{i=1}^n\min\\{U_iU_i^T,c\\} \end{equation} is invertible where $\min\\{A,c\\}$ represents the matrix with element $\min\\{a_{ij},c\\}$. Given $c>0$ and define \begin{equation} b=\frac{\sum_{i=1}^n\max\\{-c,\min\\{X_iY_i,c\\}\\}}{\left(\sum_{i=1}^n\min\\{X_i^2,c\\}\sum_{i=1}^n\min\\{Y_i^2,c\\}\right)^{1/2}}. \end{equation} The problem is how to prove that \begin{equation} -E(\min\\{\sigma\tau Z^2,c\\})\leq b\leq E(\min\\{\sigma\tau Z^2,c\\}), \end{equation} where $Z$ is normal standard distribution. I don't know how can I start this problem. Any advice?


title: recursive model

text: My paper includes two equations Log fdi =(Investment profile,exchange rate volatility, interest rate differential,corruption,bureaucratic quality,telephone) Corruption=f(Bureaucratic quality, law and order political regime,GDP growth,taxtogdp,government stability) I am confused as to which model to apply? can someones suggest should i go for two stage least square? Also i believe that my model is recursive, on that basis how should we estimate a recursive model. Is the methodology then be to appy OLS on both equations?


title: Normalization of attributes in test set for a neural network

text: I have built a neural network classification model and each attribute is z-normalized when building the model by the following X=(x-training set sample mean)/training set sample deviation. If I now have new instances that I want to classify, how should I normalize the row data? Doing the same transformation by using the training set sample mean and sample deviation? What is the best practices in handling this type of issue?


title: Ecalp function on R

text: I am doing some model validation and calibration of training data fit to use with glm and gam models, and I am having a problem with the ecalp function. Does anyone had this problem before? Any solution               Model_calibration_table_TRAIN<-data.frame()     DATA2<-data.frame(TRAIN$Latitude, TRAIN$SPECIES, TRAIN$GLM_pred, TRAIN$GAM_pred)     Model_calibration_table_bnd<-data.frame("Model"=c("GLM", "GAM"))     for (m in 3:4){     w<-ecalp(DATA2[,m], DATA2[,2])     w<-as.data.frame(w)     w$bin<-c(0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95)         caltest<-lm(w~bin, data=w, subset=w>-1)         ctest<-cor(DATA2[,m], DATA2[,2])         Model_calibration_table_train$COR[m-2]<-ctest     Model_calibration_table_train$calibration[m-2]<-caltest$coefficients[2]     Model_calibration_table_train$bias[m-2]<-caltest$coefficients[1]     }      > Error: can't find the function "ecalp"


title: Designing a uniformly most powerful test.

text: > Let $X_1, ... , X_n$ be iid normally distributed random variables $N(\mu, > \sigma^2)$, $\mu \in \Bbb{R}$, $ \sigma^2 > 0$. > > a) Design a uniformly most powerful test with significance level $\alpha$ > for testing $H_0: \sigma^2 = \sigma^2_0$ vs $H_1: \sigma^2 > \sigma^2_0$. > > b) Give a formula for the power of the test. Let $\sigma^2_1 > \sigma^2_0$. $\Lambda = \frac{L(\sigma^2_0)}{L(\sigma^2_1)}$ = $\frac{(\frac{1}{\sqrt{2\pi}})^n (1/\sigma^2_0)^{n/2} e^{-\sum (X_i - \mu)^2/\sigma^2_0}}{(\frac{1}{\sqrt{2\pi}})^n (1/\sigma^2_1)^{n/2} e^{-\sum (X_i - \mu)^2/\sigma^2_1}} \leq k$ for some $k < 1$. $\implies e^{\frac{-\sum (X_i - \mu)^2}{2\sigma^2_0} + \frac{\sum (X_i - \mu)^2}{2\sigma^2_1}} \leq (\frac{\sigma^2_0}{\sigma^2_1})^{n/2}k$ Let $k_1 = (\frac{\sigma^2_0}{\sigma^2_1})^{n/2}k$. We have $ \frac{-\sum (X_i - \mu)^2}{2\sigma^2_0} + \frac{\sum (X_i - \mu)^2}{2\sigma^2_1} \leq \ln(k_1)$ $\implies \sum (X_i - \mu)^2[\frac{1}{2\sigma^2_1} - \frac{1}{2\sigma^2_0}] \leq \ln(k_1)$ $\implies \sum (X_i - \mu)^2 \geq \ln(k_1)[\frac{1}{2\sigma^2_1} - \frac{1}{2\sigma^2_0}]^{-1}$ $\implies \frac{\sum (X_i - \mu)^2}{\sigma^2_0} \geq \ln(k_1)[\frac{1}{2\sigma^2_1} - \frac{1}{2\sigma^2_0}]^{-1}(\frac{1}{\sigma^2_0})$. Let $k_2 = \ln(k_1)[\frac{1}{2\sigma^2_1} - \frac{1}{2\sigma^2_0}]^{-1}(\frac{1}{\sigma^2_0})$. We know that $\frac{\sum (X_i - \mu)^2}{\sigma^2_0}$ has chi-square with n degrees of freedom. So we choose $k_2$ such that $P_{H_0}( \frac{\sum (X_i - \mu)^2}{\sigma^2_0} \geq k_2) = \alpha$ b) $P_{H_1}( \frac{\sum (X_i - \mu)^2}{\sigma^2_0} \geq k_2)$ $= 1 - \int^{k_2}_0$ $(\frac{1}{\Gamma(n/2)2^{n/2}})x^{n/2-1}e^{-x/2} dx$. Are my answers correct? Thanks in advance


title: Nearly constant time series

text: I want to analyse temporal interactions of some time series by means of the Box-Jenkins approach to find out which time series are predictors of another one (with the help of prewhitening and computing cross correlations for some lags) and finally set up regression models. Besides I want to analyse the same time series in a vector autoregressive (VAR-) model approach. When I looked at the plotted series in a first step I saw that some series are nearly constant over time. I have for example 140 time points and except for 10 time points the value is always 0. And the 10 values which are not 0 are nearly 0. My question is: does it make sense to include such nearly constant time series in such multivariate time series models? If so how do I have to interpret the results? Would be very happy if anyone can help me! Thank you very much.


title: Is this JAGS model ok, and can it be made faster?

text: I'm very new to Bayesian analysis, and I've come up with the following model. My goal is to get for each individual "test unit" a distribution that describes the lift in success rate under one of several test conditions.               #data     N <- 50000 #Number of test units     T <- 2 #Number of test conditions     succs<-structure(c(...N x T...))     trials<-structure(c(...N x T...))          #model     model {     for (j in 1:T) {         TestCondition[j] ~ dnorm(0, TestCondition.tau)     }     for (i in 1:N) {         Unit[i] ~ dnorm(0, Unit.tau)     }     for (i in 1:N) {         for(j in 1:T){             succs[i,j] ~ dbin(p[i,j],trials[i,j])             eps[i,j] ~ dnorm(0, eps.tau)             logit(p[i,j]) <- mu + Unit[i] + TestCondition[j] + eps[i,j]             for(k in 1:j-1){ #Only compute pairwise when non-redundant                 delta[i,(j-1)*(j-2)/2+k] <- p[i,j]-p[i,k]             }         }     }     mu ~ dlogis(0, 1)     TestCondition.tau ~ dgamma(2,2)     Unit.tau ~ dgamma(2,2)     eps.tau ~ dgamma(2,1)     }      1) Is the above model appropriate for the problem I'm trying to solve? 2) The model takes a long time for large datasets (which is basically all of my datasets). Are there any ways to reformulate the model so it can handle large datasets better?


title: How to interpret these indices/metrics for comparing partitions intuitively out of these images?

text: Two sets of comparisons were performed between original clustering and the new clustering using several indices and metrics of performance. Below are the two initial clusterings or partitions (these should be the truth or original partitions), they are two because samples were taken from two different locations and that's why we have always two sets of partitions: ![enter image description here](http://i.stack.imgur.com/gJxEl.png) The image below shows the new clusterings obtained by utilizing only a subset of data for each location instead of the whole dataset and then were compared with the original ones (obtained using the whole dataset), i.e., left original with left new, right original with right new.   ![enter image description here](http://i.stack.imgur.com/o5yGe.png) The two partitions were compared with their original coutner part (left with left, right with right) as I said and below were the results:                       Name    Part1 (Lt) Part2 (Rt)     1        Purity 0.9633028  0.7431193     2            VI 0.2451685  1.1486369     3           NMI 0.8673525  0.4062956     4    Split/Join 8.0000000 62.0000000     5 Adjusted Rand 0.8750403  0.2131243     6          Rand 0.9374788  0.6083928     7       Jaccard 0.8800522  0.4045466      **Note:**   `VI` short for "Variation of Informaiton"   `NMI` short for "Normalized Mutual Information"   `ARI` short for "Adjusted Rand Index"   `split/join` is also known as `van Dongen S` So it clear that the left new partition is much similar to its left original counterpart partition. The right new partition was way different from its right original counterpart. Samples were numbered from 1 to 34 (y-axis) and were sampled at five different occasions (x-axis). Purity measure was in line with the above and is easy to grasp. The other indices/metrics as far as I know, each measures a different thing. But to be honest as a non-statistician I cannot relate to the extensive math behind each of them. So my request from the community is how to interpret these indices with respect to the partitions in the images below without using _fancy_ statistical terms (I have already came across a bunch of them elsewhere) but rather in simple plain English. For example (pls correct me if I am wrong), I learned the following:   1. Jaccard distance measures how far one partition from the other, it can be obtained by 1-Jaccard index, which makes sense    2. for `VI`; it measures rather difference than similarity, its values is not between 0 and 1 like `NMI`, it is something between 0 and not more than $2log k$. $k$ is the number of clusters so if 4 as in this case the value should not be more than 2.log4 = 2.772589.   3. `NMI` is between 0 and 1 (NMI = 0, two partitions contain no information about one another, NMI = 1, two partitions contain perfect information about one another)   4. `ARI` is the same as Rand Index but corrected for _chance_. This measure is zero when the Rand index takes its expected value, and has a maximum of one. **Question:**   I would be very grateful to have an intuitive interpretation for Split/Join metric (as proposed by Micans, also known as van Dongen S metric), VI what does it measure exactly in relation to the images of the example shown, Jaccard how different is it from Rand? Pls feel free to correct any misconceptions posted in this question. I want to grasp the intuition behind and learn from you the correct way to interpret these measures and how they could help understand partitions from their different angles of view. Having results is fine but commenting on them is an art requires deep understanding of the tools and the philosophy behind them, may I find this with you? I hope.


title: The right distance for the clustering. Maybe Mahalanobis?

text: I have to do a cluster analysis and I'm asking which distance should I used. I know that 99% of the clustering are made using a euclidean distance, but I heard about the Mahalanobis distance and it seems to be better because it takes into account the covariance matrix of the data. **Question** : Why the Mahalanobis distance isn't more used ? For instance with this data (70% of the variance within these 2 Dim) :![enter image description here](http://i.stack.imgur.com/jP7gz.png) The euclidean distance doesn't fit, so does the Mahalanobis distance can better fit ? **Edit :** By the euclidean distance doesn't fit I mean the clusters which become apparent haven't a circle shape ![enter image description here](http://i.stack.imgur.com/2xNCR.png)


title: How to test that observations are related to time?

text: I have a vector of observations (call it $Y$) and a vector of corresponding time (call it $T$). I want to test whether the change of observations are related to time. The errors here are i.i.d and mean is 0. Is there any method or test can do this? I guess I can try polynomial regression and then do the full model F-test. Is this a reasonable way to do this? Can you give me any other suggestions to do this? Please also give me any assumptions about the observations, $H_0$ and test required by your method.


title: Sample size needed to find all error values

text: This is a pretty basic question as I'm new to stats, but hopefully someone can help! I have a sample size of 30,000 records and so far 6500 have been checked to see if an error exists - of these 6500, 1124 have an error. What I would like to be able to do is the following: 1) See how many records I would need to check to find 97% of all error records (this is a deliberately fixed value) 2) Determine if by checking 6500 records I had actually checked enough to have found all of the (97% of the total) error cases. My question is similar to one found here but is different enough for me to be unsure of how to apply it! TIA, and hope it makes sense! Davin **edit to include sample data** The data in question can be found here:


title: Is a logistic regression biased when the outcome variable is split 5% - 95%?

text: I am building a propensity model using logistic regression for a utility client. My concern is that out of the total sample my 'bad' accounts are just 5%, and the rest are all good. I am predicting 'bad'.   * Will the result be biassed?    * What is optimal 'bad to good proportion' to build a good model?


title: Advice on missing value imputation

text: I am working on insurance data in which a customer has a field named `customer_no_dependent` (customer's number of dependent). Its coming out to be a significant variable( just that it has $p<0.0001$). This variable has almost 20% missing values. For imputation, I thought to determine proxy indicators for number of dependents. I tried age (thinking a person more aged could be having more dependents). I correlated it with premium amount as well to think that a person who has more dependents could be having less disposable income. So low premium paying could be meaning more dependents. I do understand that a demographic variable can't be fully taken out from such logic. Now, if somebody goes into detail, he can prove my imputation to be far from perfect. What should I do in such situation? Would deleting those 20% be a correct solution? 20% for my data would be close to 2 lakh rows which is large amount of information.. I know, this question can have many possible answers. I would be grateful for any pointers how to proceed.


title: How to analyse repeated measure ANOVA with three or more conditions presented in randomised order?

text: ### Context: My question concerns a typical design in my area – a researcher takes a group of subjects (say 10) and then applies three different conditions to them to measure the change in a response variable, e.g. vertical jump height performed after drinking a glucose drink, coloured plain water, and fruit juice (say). Every subject has every treatment, but in a random order with enough time between for effects to ‘wash out’. ### Analysis: Kuehl (2000) (Kuehl, R. O. (2009) Design of Experiments: Statistical principles of research design and analysis, Duxbury Press, CA, p497 2nd Ed.) states: > When each of the treatments is administered in a random order to each > subject... then subjects are random blocks in a randomised complete block > design” and then shows the corresponding analysis. In this case, the subject is a random effect, but a nuisance or blocking factor, and although our statistical model will test the significance of the block factor, we aren’t really interested in its significance. However, many researchers (and reviewers!) think that such a design should be analysed as a repeated measures design with a Mauchly test for the Huynh-Feldt condition (with the treatment as the repeated measure). However, this seems more appropriate when a time factor is being analysed – for example when observations are taken at 0 minutes, 10 minutes, 30 minutes and 60 minutes, for example. In this case the covariance between pairs of time points might reasonably be expected to change, particularly when unequal time intervals are used. [In fact, I use SAS to model different covariance structures in this case (e.g. autoregressive) and use the AIC to choose the best structure, though this is not an approach that is well received by many reviewers.] I understood that when the subject is a block factor, and the different treatments are administered in a random order that is different for different subjects, this means that the correlation between observations is different for each subject so compound symmetry can be assumed. ### Question:   * How should repeated measures ANOVAs with 3 or more conditions presented in random order be analysed?   * Is it reasonable to assume compound symmetry?


title: Generating dependent time series from a given distribution?

text: How can I generate dependent time series from a given marginal distribution? I want to be able to adjust the level of dependence, to influence the predictability of the series, which will be given as input to a Monte Carlo simulation. The dependence parameter can be the correlation, the mutual information, or something along those lines. You may assume the distribution is Bernoulli for discussion purposes. MATLAB code is gratefully accepted.


title: How can I interpret what I get out of PCA?

text: As part of a University assignment, I have to conduct data pre-processing on a fairly huge, multivariate (>10) raw data set. I'm not a statistician in any sense of the word, so I'm a little confused as to what's going on. Apologies in advance for what is probably a laughably simple question - my head's spinning after looking at various answers and trying to wade through the stats-speak. I've read that:   * PCA allows me to reduce the dimensionality of my data   * It does so by merging / removing attributes / dimensions that correlate a lot (and thus are a little unnecessary)   * It does so by finding eigenvectors on covariance data (thanks to a nice tutorial I followed through to learn this) Which is great. However, I'm really struggling to see how I can apply this practically to my data. For instance (this **is not** the data set I'll be using, but an attempt at a decent example people can work with), if I were to have a data set with something like...               PersonID     Sex     Age Range    Hours Studied     Hours Spent on TV      Test Score     Coursework Score      1            1       2            5                 7                      60             75     2            1       3            8                 2                      70             85      3            2       2            6                 6                      50             77     ...          ...     ...          ...               ...                    ...            ...      I'm not quite sure how I would interpret any results. Most of the tutorials I've seen online seem to give me a very mathematical view of PCA. I've done some research into it and followed them through - but I'm still not entirely sure what this means for me, who's just trying to extract some form of meaning from this pile of data I have in front of me. Simply performing PCA on my data (using a stats package) spits out an NxN matrix of numbers (where N is the number of original dimensions), which is entirely greek to me. How can I do PCA and take what I get in a way I can then put into plain english in terms of the original dimensions?


title: Is my (SPSS) independent t test correct, and how should I intrepret it?

text: I have a sample size of 158 students belonging to IT(N=77) and non-IT disciplines(N=81). I have carried out the independent sample t-test but I am not sure whether it's correct or not. Moreover, how do I interpret in my writing that the group differences have impact on other factors.   ![enter image description here](http://i.stack.imgur.com/jBSAm.png)   ![enter image description here](http://i.stack.imgur.com/y4buk.png) I'd also appreciate if you're able to share a reference paper explaining the values to focus on.


title: Increasing Exam Expected Mark

text: You are in an exam, and are presented with the following question: > Write down what mark do you expect to take in this exam... If you get it > right in range of +/-10 % then you will take 10% bonus... if wrong (or not > answered) you will lose 5% Assume that you have no idea of how you are going to perform in this exam. How would you choose a mark that maximizes your expected return mark? In other words, I need to help deriving an equation for the expected mark given a decided mark, if that is possible...