title: How to test/prove data is zero inflated?

text: I've got a problem that I think should be simple but can't quite figure it out. I'm looking at seed pollination, I have plants (n=36) that flower in clusters, I sample 3 flower clusters from each plant, and 6 seed pods from each cluster (18 seed pods in total from each plant). A pod can have between 0 to at most 4 seeds pollinated. So, the data is count, with an upper bound. I'm finding an average of ~10% of seeds are pollinated, but anywhere between 1 - 30% on a given plant, so over dispersed data, and of course, there are 4 missing cluster replicates on 3 plants, so not perfectly symmetrical. The question I'm asking is if this data supports the idea this plant requires pollinators for seed set. I'm finding that the distribution for the number of seeds in a pod looks like there are more 0 pollinated seed pods (6-9 pods out of 16) and more 3 and 4 pollinated seed pods (2-4 for each) than would be expected if seeds in the population were just randomly pollinated. Basically, I think this is classic example for zero inflated data, first a insect either does or does not visit the flower at all (one zero generator) and if it does, then pollinates 0-4 of the seeds in another distribution. The alternative hypothesis is the plant is partially selfing, and it would then be expected that every seed would have the same probability of being pollinated (this data suggests a roughly 0.1 chance, which means 0.01 chance for two seeds in the same pod, etc). But I simply want to demonstrate the data best fits one or the other distribution, not actually DO an ZIP or ZINB on the data. I think whatever method I use should take into account the actual number of pollinated seeds and the number of pods sampled on each plant. The best thing I have come up with is to do some sort of boot strap thing where I just randomly assign the number of pollinated seeds for a given plant into the number of seed pods I sampled, do that 10,000 times and see how likely it is the experimental data for the given plant came out of that random distribution. I just feel there is something about this that should be a lot easier than brute force bootstrapping, but after days of thinking and searching I’m giving up. I can’t just compare to a Poisson distribution because it’s upper bound, it’s not binomial because I need to generate the expected distribution somehow 1st. Any thoughts? And I’m using R so advice there (especially how to most elegantly generate 10,000 random distributions of n balls into 16 boxes that can each contain at most 4 balls) would be most welcome. ADDED 9/07/2012 First, thanks to you all for all the interest and help. Reading over the answers has made me think to reword my question a bit. What I’m saying is that I have one hypothesis (which for now I am thinking of as the null) that seeds are pollinated randomly across pods, and my alternative hypothesis is that a seed pod with at least 1 pollinated seed is more likely to have multiple pollinated seeds than would be expected by a random process. I’ve provided real data from three plants as examples to illustrate what I’m talking about. First column is the # of pollinated seeds in a pod, second column is the frequency of pods with that seed count. plant 1 (total 3 seeds: 4% pollination) num.seeds ::pod.freq 0::16 1::1 2::1 3::0 4::0 plant 2 (total 19 seeds: 26% pollination) num.seeds::pod.freq 0::12 1::1 2::1 3::0 4::4 plant 3 (total 16 seeds: 22% pollination) num.seeds::pod.freq 0::9 1::4 2::3 3::2 4::0 In plant #1, only 3 seeds were pollinated in 18 pods, one pod had one seed, and one pod had two seeds. Thinking about a process of adding one seed to the pods at random, the first two seeds each go into their own pod, but for the 3rd seed, there are 6 spots available in pods that already have one seed but 64 spots in the 16 pods with no seeds, so the highest probability of a pod with 2 seeds here is 6/64= 0.094. That’s a bit low, but not really extreme, so I’d say that this plant fits the hypothesis of random pollination across all seeds with a ~4% chance of pollination occurring. But plant 2 looks much more extreme to me, with 4 pods completely pollinated, yet 12 pods with nothing. I’m not quite sure how to calculate the odds of this distribution directly (hence my bootstrap idea) but I’d guess the odds of this distribution occurring at random if each seed has a ~25% chance of pollination are quite low. Plant #3 I really have no idea, I think there are more 0’s and 3’s than one should expect for a random distribution but my gut feeling is that this distribution for this number of seeds is much more likely than the distribution for plant #2, and may not be that unlikely. But obviously I want to know for sure, and across all plants. In the end I’m looking to write a statement like “The distribution of pollinated seeds in seed pods fits (or does not fit) the hypothesis that plants are not simply partially self compatible, but require visitation of a pollinator to effect seed set. (results of statistical test).” This is really just part of my forward looking section, where I’m talking about what experiments to conduct next, so I’m not desperate for this to be one thing or the other, but I want to know for myself, if possible. If I can’t do what I’m trying to do with this data, I’d like to know that too! I did ask a rather broad question at first, since I’m curious as to whether or not there are any good tests to show if data should go into a zero inflated model in the first place. All of the examples I’ve seen seem to say –“look, there are a lot of zeros here, and there is a reasonable explanation for that, so let’s use a zero inflated model”. That’s what I’m doing right now on this forum, but I had an experience on my last chapter where I used a Poisson glm for count data and my one of my supervisors said “No, glms are too complex and unnecessary, this data should go into a contingency table” and then sent me a data dump of the massive contingency table generated by their expensive stats package that gave the same p values for all my factors + interactions to three significant digits!! So, I’m trying to keep the stats clear and simple, and make sure I understand them well enough to robustly defend my choices, which I don’t feel I can do for a zero inflated model right now. I’ve used both a quasibinomial (for whole plants to get rid of pesudoreplicaiton) and a mixed model for the above data to compare treatments and answer my main experimental questions, either seems to do the same job, but I’m going to also play around with ZINB’s tonight, to see how well that performs. I’m thinking if I can explicitly demonstrate that this data is strongly clustered (or zero inflated) at first, then provide a good biological reason for that occurring, I’ll be much better set up to subsequently pull out a ZINB, than to just compare one to a quasibinomial/mixed model and argue since it gives better results, that’s what I should use. But I don’t want to distract too much from my primary question, how can I determine if my data really is more zero inflated than expected from a random distribution? In my case the answer to that is what is of real interest to me, with the possible benefit for model justification being a bonus. Thanks again for all your time and help! Cheers, BWGIA


title: Finding a comparable Control group for a treatment group?

text: I have a treatment group of size 30 (30 schools in California) that used a math supplemental software. In a simple analysis, I'd like to compare students' average Math growth between our treatment group and a comparable control group. There are many schools in CA that didn't use the software. I'd like the control group to include similar performing schools (their baseline scores be similar to treatment schools with a reasonable margin of error). Also, I'd like the control group sample size be 3 times of my treatment (here 90 schools). There are many many choices of 90 schools out of more than 1000 schools in CA. How would you choose your control group?


title: How to know if a data follows a Poisson Distribution in R?

text: I am an undergrad student and I have a project for my probability class. Basically, I have a dataset about the hurricanes that impacted my country for a series of years. In my probability Book, (Probability ans Statistics with R) there is an (not complete) example of how to check if the data follows a poisson distribution, they begin trying to prove that this 3 criteria are followed: (From my book, page 120 (criteria) page 122-123 example) 1- The number of outcomes in non-overlapping intervals are independent. In other words, the number of outcomes in the interval of time (0,t] are independent from the number of outcomes in the interval of time (t, t+h] , h > 0 2- The probability of two or more outcomes in a sufficiently short interval is virtually zero. In other words, provided h is sufficiently small, the probability of obtaining two or more outcomes in the interval (t,t+h] is negligible compared to the probability of obtaining one or zero outcomes in the same interval of time. 3- The probability of exactly one outcome in a sufficiently short interval or small region is proportional to the length of the interval or region. In other words, the probability of one outcome in an interval of length h is lambda*h. But criterion 3 is left "as an exercise". A- Can someone tell me if there is a more "easy" way to see if my dataset follows a Poisson distribution? B- Can someone explain me criterion 1 and 3 with some type of example (if it is with R, fantastic)? Thanks! Note: Sorry for the long post. Also, I have to convert the data so that I have a table like:                 number of hurricanes       | 0 | 1 | 2  etc.       -----------------------------------------       total years that have      |   |   |       that number of hurricanes  |   |   |


title: How to normalize time series?

text: This is a general question on normalization of data so that all the variables are within the same range.   1. Why do we normalize data in pattern classification?   2. How to normalize time series which is multivariate? I have collected motion time series from Kinect skeleton tracking. In general, what is the process of normalization when the features are the joint angles, x,y,z coordinate of the pixel and in cases of features say flower width, size, color ?


title: How to visualize phylogeographical data?

text: I am working on phylogeography on one model species and I am a beginner. You can imagine one species that came on locality more than 5 mya (million years ago). It was a good environment, so the species expanded very quickly. After climatic changes this model species split in two refugiums (isolated locations for a once widespread species) and evolved into two species. This happened several times and now we have from that one ancestral species six valid species. I have data from MrBayes. Can you recommend some phylogeographical program which can consider also expansion and contraction of model species? **My input data:**   * coordinates (lat, lon) for each species   * geographical data of a continent (altitude, precipitation etc.) **My goal:** Create a map where will be polygons (areas) of each species. These areas will expand and contract (step-by-step not like a do not want create video). **My problem:** I am able to make polygons by MCP, or kernel density, but unable to incorporate barriers into model. For example, polygons of species which lives in mountain region created by MCP include areas of lowland inside...


title: Statistical Methods for Calculating Vending Machine Refill

text: Am looking into statics to help support a project I am undertaking. The project scope concerns **intelligent replenishment / refill of vending machines**. During an onsite service, a technician must make decisions regarding machine refill to optimise sales for the period up to the next refill visit. The vending machines have no call back to base ability. The technician must…   * maximise sales by specifically stocking and selling items that demonstrate superior sell through (this may vary over time and seasonally)   * limit the restocking of inventory for high value items like cigarette packets or perishable items involving use by dates so over supply does not occur, but balance this by...   * ensuring sales are not lost through insufficient restocking of inventory items The statistical model or programmatic code needs to consider the following real world factors…   * Factors impacting short term sales (seasonal factors, gazetted holiday, localised events, machine inventory levels)   * Factors impacting long term sales trends (long term seasonal factors, machine location)   * Sales performance per bay location within a vending machine   * Item properties including use by dates, purchase price, sales price Was wondering whether a definitive highly accurate statistical method existed to achieve these outcomes? I was planning to construct an array of intelligent algorithms to overcome these parameters with the intent of operation within PDA styled smart client devices. A colleague suggested modelling data using Poisson regression but I gather that is just a starting point to overcome all the interwoven variables?


title: What statistical techniques to use to distinguish between two groups?

text: I have a dataset of about 300 people. 200 test positive for a disease, and the rest test negative. I have data on different test scores and imaging results for these 300 participants. So my dataset would look something like this                   status    test1   test2   test3   imaging1         pos        10      10      5        98         neg        8       7       5        77         pos        8       9       5        98         pos        10      10      5        99         pos        10      10      5        100         neg        6       8       4        78      And so forth. Is there a technique that tells me whether the two groups are different, and if so, how different? I know of t-test, but is there something else that can tell me that the positive group is different from the negative group (if so, what %)?


title: Calculating probability

text: If $f(x,y)=2x , 0\leq x\leq 1 ,0\leq y\leq 1 $, find $ P(Y < e^{-X} \cap X > Y)$ Given X and Y have joint distribution. Here is my approach: $$ P(Y < e^{-X} \cap X > Y) = 1- P(Y > e^{-X} \cap X > Y)$$ = $1- P(X > Y > e^{-X})$ = $P(e^{-X} \leq Y \leq X )$ =$ \int_0^1\int_{e^{-x}}^x2x dydx $ = $2e^{-1} - \frac{1}{3}$ Is my solution for this problem correct? Thanks!


title: Is a sample of 1961 people out of a population of 225,000,000 a sample size that is too small? If so, why?

text: I am working by a set of problems on a single question in my textbook. The question explained that there was a recent survey of 1961 adults out of a population of 225,000,000 and asked me to construct a 95% confidence interval estimate of the percentage of all adults on the "Do not call registry". I completed this problem quite easily using the methods prescribed previously in the text. The following problem asked me if the sample size of 1961 was unreasonably small for the population size, which I do not know how to determine. I am starting to think that the question simply wants me to use my own opinionated reasoning.


title: Search for similar documents on growing term-base

text: Could you please recommend any approach (if it exists) for similar document search considering the following:   1. It's necessary to provide term base explicitly. This is a dictionary of specific named entities on which searching are performed.    2. This term base is unstable - new terms can be added to the dictionary. So what I need to do is to build a some feature vector for each document which allows fast comparing documents using provided term base. The critical thing - search must working even provided (search) document's feature vector constructed on newer version of term base than a feature vectors was built.   I understand that question is not only about feature vectors, it's mostly about comparison mechanism. Is it exists any approach using LSA/PLSA to build such feature vector? Any suggestions?   Thank you.


title: Deep learning ... notion of depth

text: I'm reading from http://www.iro.umontreal.ca/~bengioy/papers/ftml.pdf around page 8. Here is a quote that I'm not understanding:               "If we include affine operations and their possible composition with sigmoids       in the set of computational elements, linear regression and logistic regression      have depth 1, i.e., have a single level."      Depth or Depth of Architecture refers "to the depth of that graph, i.e. the longest path from an input node to an output node."![Here is a simple example demonstrating a depth of 4](http://i.stack.imgur.com/p9dWq.png) I can see why the graph represented in the image has depth of 4. I'm not sure how to think about, say, linear regression in these terms.


title: How to calculate standard error of sample quantile from normal distribution with known mean and standard deviation?

text: I know that the standard error of the mean for an iid sample is calculated as $$\frac{\sigma}{\sqrt{n}}$$ However, assuming a normal distribution with known mean and standard deviation, how do you calculate the standard error of an arbitrary quantile? For example, assume   * normal distribution   * population mean = 0   * population standard deviation = 1   * n=100   * quantile = .95 What would be the standard error of this quantile? * * * I ran this little simulation to explore the properties, but I'm still interested in the closed form solution:               set.seed(1234)     generate_x <- function(n) x <- rnorm(n)     k <- 10000     n <- 100          results <- lapply(seq(k), function(X) generate_x(100))          Z <- seq(.01, .99, .01)     qresults <- sapply(results, function(X) quantile(X, Z))          sd_quresults <- apply(qresults, 1, sd)     var_quresults <- apply(qresults, 1, var)          plot(Z, sd_quresults, type='l')      ![enter image description here](http://i.stack.imgur.com/mpXTD.png)


title: Good online resource with tips on graphing association between two numeric variables under various conditions

text: ### Context: Over the while I've acquired a set of heuristics on how to effectively plot the association between two numeric variables. I imagine most people who work with data would have a similar set of rules. Examples of such rules might be:   * If one of the variables is positively skewed, consider plotting that axis on a log scale.   * If there are a lot of data points (e.g., n > 1000), adopt a different strategy such as using some form of partial transparency, or sampling the data;   * If one of the variables takes on a limited number of discrete categories, consider using a jitter or a sunflower plot;   * If there are three or more variables, consider using a scatterplot matrix;   * Fitting some form of trend line is often useful;   * Adjust the size of the plotting character to the sample size (for bigger n, use a smaller plotting character);   * and so on. ### Question: I'd like to be able to refer students to a web page or site that explains these and other tricks for effectively plotting associations between two numeric variables, perhaps with examples.   * Are there any pages or sites on the internet that do a good job of this?


title: How to do one-class text classification?

text: I have to deal with a text classification problem. A web crawler crawls webpages of a certain domain and for each webpage I want to find out whether it belongs to only one specific class or not. That is, if I call this class _Positive_ , each crawled webpage belongs either to class _Positive_ or to class _Non-Positive_. I already have a large training set of webpages for class _Positive_. But how to create a training set for class _Non-Positive_ which is as representative as possible? I mean, I could basically use each and everything for that class. Can I just collect some arbitrary pages that definitely do not belong to class _Positive_? I'm sure the performance of a text classification algorithm (I prefer to make use of a Naive Bayes algorithm) highly depends on which webpages I choose for class _Non-Positive_. So what shall I do? Can somebody please give me an advice? Thank you very much!


title: How to quantitatively analyse a semi-structured interview?

text: I am looking to design a study that involves interviewing children by asking them questions that relate to their values, practices, and identities. Children will be migrants who could identify to varying extents with both their heritage culture and adopted (receiving) culture values, practices, and identities. I'm a bit stuck on thinking of a good way to organize and analyze the data quantitatively. Variables:   * Values (influenced by heritage and receiving culture, + interaction)   * Practices (influenced by heritage and receiving culture, + interaction)   * Identifications (aligned with heritage or receiving culture, or both). Any suggestions would be most welcome!


title: Estimating how quickly people submit their answers

text: In a contest for solving a problem, a person submits his/her answer for the problem at any time during the contest. The contest system provides the person with feedback on the submitted answer. The person can resubmit a modified answer based on the feedback at any time during the contest. This process is repeated until the end of the contest. I have a dataset for the contest as follows:               Table 1         person feedback submission_time            A     50             20            A     80             50            A     90             60            B     10             35            B     11             80         ...      My interest is to estimate how quickly people submitted their answers based on feedback. In other words, how does feedback at the previous time affect the probability of submitting a revised answer? For this purpose, I first considered a survival model and modified the above model as follows:               Table 2         person feedback time0 time1 submission            A      .        0    20      1                 A     50        20   50      1            A     80        50   60      1            A     90        60           0         ...      As seen in the above table, every submission is present except for the final submission. In this case, does a survival model such as the Cox proportional hazard model provide the proper estimates for my research objectivw? Secondly, are there more suitable statistical models to address my objective than a survival model?


title: Performing a Sobel test with huge discrepancies between unstandardized coefficients?

text: I have read pretty much everything I can find on performing a Sobel test on the internet however the huge discrepancy in my unstandardized coefficients makes me feel like I am doing something completely wrong (a = 0.000306, SEa = 0.000067, b = 430037.87. SEb = 258407.87). Did I miss a step somewhere? Thank you for any insights. ![http://i.imgur.com/7ktbT.png](http://i.stack.imgur.com/o1v1r.png)


title: Logistic regression: the standard deviation used in: GLMPOWER

text: I attended a training class from SAS about experimental design in marketing. They advocated the use of their GLMPOWER proc for power analysis for designing experiments. GLMPOWER is a power analysis procedure for ascertaining the required sample size for a general linear model (main effects, interactions and/or specific contrasts between factor levels). For this class, interest rests on designing experiments with regards to binary outcomes (response or no response) analyzed using logistic regression. Given the large sample sizes typically witnessed in marketing, they said an approximation to this problem could be handled by GLMPOWER, which assumes the response is a Gaussian distributed continuous variable. One of the parameters in the GLMPOWER procedure is a standard deviation STDDEV defined as "the error standard deviation, or root MSE" of the model being postulated. **Questions:**   1. This error standard deviation is not the same as the standard deviation of "Y" is it? They state that it is.   2. Further, after stating this, they use SQRT(p(1-p)) as an estimate of this value, where p is the pooled response rate (number of responses / number of attempts) of the experiment. Where does this come from and does it sound like the right value to use? It was well stated that these are approximations only given the assumptions of GLMPOWER for a binary outcome, but can anyone help with the reasoning?


title: Texts on Various Topics in Statistics (GLMs, MCMC, Decision Trees, etc.)

text: I am currently looking for texts (or preferably a specific text) which have a good balance between theory and application and are as comprehensive as possible and are at an _introductory level_ , covering these topics: > Markov Chains, MCMC, GLMs and extended linear models, Decision Trees, > Spatial Statistics, Visualization & Fit, Neural Nets, Ensemble Methods From my research, it seems that these are topics covered in a master's level statistics program, but I've only had an undergraduate background and don't know what are the standard texts for such topics. Edit: Needless to say, I would like more detailed answers/feedback on this topic. Edit2: The answers below are nice, but I would like answers which address all topics if possible. As mentioned in the bounty description, this list of topics will be covered in a future actuarial exam. No source material has been announced for this exam yet.


title: Repeated measures ANOVA design/structural issue

text: As part of my PhD work, I've conducted an inoculation experiment concerned with marine phytoplankton community productivity (dependant variable, as 'no. of cells') vs nutrient availability. I have two factors:   1. Nutrient Treatment (4 levels - including control)   2. Time (sampled at 5 points) The puzzling aspect to this design is that within the phytoplankton community,there are 4 main taxonomic groups (each with 'no. of cell' data), which I need to account for (i.e. I would like to know which of these groups significantly differed both within and between treatments over the entire sampling regime). So, I understand a post-hoc test (e.g. TukeyHSD) will be needed to run multiple comparisons test in order to elucidate group differences, but first need to know which test to initially run (i.e. Two-way repeated measures ANOVA, or multiple ANOVA), and in what structure? Or, is it simply not possible to simultaneously analyse for between phytoplankton group differences?


title: Estimate median from mean, std dev, and/or range

text: I have easy access to the the mean, standard deviation and min/max values for a variable. I also have the number of elements used to compute those values. The data is real numbers with an absolute zero. Given this information, is it possible to estimate the median value?


title: Compared to the ρ in a standard Dickey Fuller test, is the calculation for ρ in an augmented Dickey Fuller Test the same?

text: As the question suggests, we have a standard AR model. Normally, the way to get ρ is through http://www.stata.com/manuals13/tsdfuller.pdf (page 4). You'll notice in the page after that, it talks about ADF where Beta is our new ρ. Are the methods to get the ρ different now?


title: Implementation of Total Variation Regularization Algorithm (Lagged Diffusivity Algorithm)

text: I am trying to compute the derivative of an experimentally-measured quantity as a function of time. The data are fairly noisy, which causes problems. For instance, using finite differences (central differences) to compute the derivative results in a derivative data set which is even noisier than the original data set. The idea ultimately is to plot the derivative of the quantity versus the original quantity to identify whether or not a power law relationship exists (we expect it to show up). Therefore, I would like to have a fairly smooth derivative if I am to have any hope of extracting something close to the actual power law relationship from the data. I have researched this subject fairly extensively over the past couple of days. I have found several useful references (see, e.g., here, here (see Chapter 8), and here) on using Total Variation Regularization (which I have decided to try using), but I'm running into trouble when I try to actually implement an algorithm. (I'm using Mathematica 9.0, and as far as I can tell, I will need to write my own function to do this.) From the first reference, I'm trying to solve for the minimizer $u^*$, which is the derivative of a function $f$ on some closed interval domain, of the following functional: $$F(u)=\alpha R(u)+DF(Au-f),$$ where $R(u)$ is a term that penalizes irregularity (or roughness, or noisiness) in $u$, $A$ is an antidifferentiation operator, $DF(Au-f)$ is a data fidelity term that penalizes differences between $Au$ and $f$, and $\alpha$ is an adjustable parameter to control the relative importance of the two terms. I'm attempting to follow the discrete implementation I in the first reference, page 3. The author says that "we assume $u$ is defined on a uniform grid $\{x_i\}_0^L=\{0,\Delta x,2\Delta x,...,L\}$. Derivatives of $u$ are computed halfway between grid points as $Du(x_i+\Delta x/2)=u(x_{i+1})-u(x_i)$. This defines our $L\times(L+1)$ differentiation matrix $D$...Integrals of $u$ are likewise computed halfway between grid points, using the trapezoidal rule to obtain an $L\times(L+1)$ matrix $A$. Let $E_n$ be the diagonal matrix whose $i$th entry is $((u_n(x_i)-u_n(x_{i-1}))^2+\epsilon)^{-1/2}$ [where $\epsilon$ is some small positive number used to avoid division by zero]. Let $L_n=\Delta x D^TE_nD$, $H_n=K^tK+\alpha L_n$. The matrix $H_n$ is an approximation of the Hessian of $F$ at $u_n$. The update $s_n=u_{n+1}-u_n$ is the solution to $H_ns_n=-g_n$, where $g_n=K^T(Ku_n-f)+\alpha L_nu_n$..." I think there are some typos in there, but this reference seems to be about the most detailed I've found on how to actually perform the algorithm. My biggest problem with this is that the dimensions don't work out, which makes me wonder if I'm not understanding what the variables are standing for. First, shouldn't $K$ be replaced by $A$ (or $A$ with $K$)? Second, the derivative definition is missing a factor of $1/\Delta x$, right? Now, on to the dimensions. For fun, let's say I have just 5 data points. (I actually have ~1500 pts.) Then, $f$ will be a vector with dimension 5x1. Similarly, I can get a $u$ vector with dimension 5x1 by using 2-pt central differencing for the interior points, and 3-pt forward and backward differencing at the two endpoints (to maintain the same order of error ($\Delta x^2$)). The $D$ matrix will have dimension 4x5, since it gives derivatives of the derivative of $f$ (or derivatives of $u$) at the midpoints of the 4 segments by considering the 5 computed values of $u$. Similarly, the $K$ (or $A$) matrix will have dimension 4x5. Matrix $E_n$ has dimension 4x4. Thus, $L_n$ has dimension 5x5, based on its definition above and consideration of the dimensions of the factors used to computed it. Likewise, $H_n$ will have dimension 5x5. The problem comes up during computation of $g_n$. $u_n$ has dimension 5x1. $K$ has dimension 4x5, so $Ku_n$ has dimension 4x1. But we are then supposed to subtract from this $f$, which has dimension 5x1, so the subtraction is undefined. Obviously I have done something wrong somewhere. I would really appreciate it if someone could help with this problem. Thanks. Note: I originally posted this on the Mathematics StackExchange site, but I'm thinking maybe this is a more appropriate place to ask the question?


title: How to calculate degrees of freedom for chi squared test

text: I am very new to stats, and if the chi squared test is not the best way to handle this, please let me know. To simplify what I'm trying to do, I have many separate lists of genes. Each list is a different size in terms of number of entries (one gene being one entry). For each list, I have counted the number of times each specific gene appears (that is my observed). So for example, in one specific list, the specific gene `KRAS` occurs three times (that is entered into a table as my observed). In a different list, it occurs 5× (another observed value). For my expected, I've pooled all lists together into one master list, and counted the number of times each word appears in the master list, and divided by the total number of entries in the master list to get its overall probability based on the whole data set, then multiplied it back to the original size of any one list to get the expected. I know how to then calculate the chi squared test statistic, but I don't know how to get the degrees of freedom from this. I can't even begin to figure out how to put this into a tabular form, but I suppose one way of representing it would be:                     List1   List2   List3   List4   List5 etc.      GeneA  O/E    O/E     O/E     O/E      O/E      GeneB  O/E    O/E     O/E     O/E      O/E      etc.        * where O/E represents two values, one being observed, the other being expected. ULTIMATELY, I just want to figure out if individual genes occurring within individual lists deviate from their occurrences within the whole data set (all lists combined). Could someone suggest how I could get the degrees of freedom for chi squared, or if there might be an easier test to use?


title: Using two one variable linear regressions on a single response variable to compare explanatory variables

text: I ran two single variable linear regressions $A$ and $B$. $A$ had a relatively large effect size $R^2$ = $.68$ while $B$ had $R^2$ = $.10$. Regressing the explanatory variable from $A$ (= $a$) on the explanatory variable from $B$ (=$b$) shows an $R^2$ = $.44$. The correlation between these two variables is $.66$. Is there a relationship between their correlation in the Pearson Correlation Matrix and in the regression? Is it at all strange that despite the variables being correlated their regressions (on the same response variable) yield significantly different results in terms of effect size?


title: How to apply feature selection based on tf-idf threshold

text: Let's say we have the following matrix (typical VSM example):                        word1    word2    word3    word4     doc1     0.838    0.918    0.001    0.482     doc2     0.311    0.948    0.001    0.183     doc3     0.001    0.928    0.592    0.381      I want to initially get rid of some words with a low tf-idf weighting, but I am not sure what to do with words as word1 in my example. word1 is worthless for doc3 but looks important for doc1. Should I take an average of these weights and put a threshold there?


title: Markov Chain State Transition Probability in R

text: I have a dataset which shows the states (3 states) across 11 time points for each participant. I wanted to estimate the Markov Chain state transition probability matrices for time points 2-11 using R. I read both `markovchain` and `HiddenMarkov` package documentations, but did not see any examples similar to what I planned to do. The most relevant example is the one listed under `MarkovchainFit` (in `markovchian` package), but it seems to me the data entry is a sequence of states (not a sequence of states for an each individual participant), and there is only one transition matrix calculated. I appreciate any suggestions. That is also great if any other statistical software can support the calculation. Thank you!


title: Shouldn't a weighted t-test and a weighted linear regression yield the same standard error on the estimated mean?

text: I have a vector of values and a vector of the reciprocal of their standard error. I'd like to estimate their mean and the standard error on the mean. At first I used linear regression using the `lm()` function. Later I found out that the `weights` library provided a `wtd.t.test()` function for doing t-tests on weighted data. One side-effect of the test is an estimate of the mean of the data, with a standard error. I expected both approaches to yield the same result, but they don't. Here is a short but complete example:               library(weights)     test <- c(1,1,1,1,1,1,2,2,2,3,3,3,4,4)     weight <- c(.5,.5,.5,.5,.5,1,1,1,1,2,2,2,2,2)     wtd.t.test(test, weight=weight)     summary(lm(test~1, weights=weight))      The weighted t-test yields 2.64 +/- 0.27, whereas the linear regression yields 2.64 +/- 0.30. Which is correct (if any)?


title: Naive Bayes with Density estimation

text: The way I understand the Naive Bayes estimators is that the characteristic (or naive) assumption is that all Features are conditionally independent. Now on top of this assumption, for continuos features it is also assumed that they are Gaussian distributed. My question is whether this is something inherent in the Naive Bayes estimator or why do we not just use a non- parametric density estimation to determine the one dimensional probability distributions of the features.


title: How to test for parallelism for two linear models?

text: I'm taking a graduate course in regression analysis and I'm suck on a particular homework question that _should_ be very simple to me! I have the following model:               E(y) = B0 + B1x1 + B2x2 + B3x3 + B4x1x3 + B5x2x3      `x3` is coded as `1` if "smoker" and `0` if "non-smoker". Therefore the regression equations are:               x3 = 1: E(y) = (B3 + B0) + (B1 + B4)x1 + (B2 + B5)x2     x3 = 0: E(y*) = B0 + B1x1 + B2x2      Now I know how to test for parallelism if `x2` is absent in the models:               H0: B4 = 0     H1: B4 != 0      But I'm lost as to what to do with the inclusion of the `x2` variable. Parallelism is obviously testing for slope, but I'm not sure where to find the "slope" coefficient. I was thinking about using an `F-Test` but then I realized I don't actually want to test the _whole_ model, just the parallelism. Could someone please point me into the right direction? Even hints would be sufficient.


title: R code for time series forecasting using Kalman filter

text: Does anybody have a good example for Time Series Forecasting/smoothing using Kalman Filter in R?


title: How to test if change is significant across multiple categories?

text: Here is a look at my data. We asked the same respondents (n=~400) to provide us with their current and future consumption as a proportion of total expenditure. Plotted here are the mean proportions of total expenditure for each category for "Now" and "Later", respectively current and future. What I'm looking for is a statistical method that I can use to test if the other categories are increasing while the largest category is decreasing. The graph shows this, to some degree, but I would like to make sure that it is statistically valid. ![Example](http://i.stack.imgur.com/o3g4r.jpg)


title: Distance between nominal vectors: a review

text: I have a seemingly easy question which however is troubling me a bit. I have couples of vectors made up of nominal attributes. They can be of different length and sometimes some of the attributes in one might not be included in the other. See `a` and `b` as two potential examples.                              a     1  mathematician     2       engineer     3  mathematician     4  mathematician     5  mathematician     6       engineer     7  mathematician     8  mathematician     9  mathematician     10 mathematician     11 mathematician     12      engineer     13 mathematician     14 mathematician     15      engineer                         b     1      physicist     2        surgeon     3      physicist     4        surgeon     5      physicist     6      physicist     7        surgeon     8        surgeon     9      physicist     10     physicist     11 mathematician      Do you have in mind a measure that could summarize the dissimilarity between them? The type of measure I am looking for is something like the euclidean distance, but for qualitative vectors. One option I have in ming is to actually compute the euclidean distance among the categorical vectors transformed into frequence vectors. In this way, they would become quantitative and would be of the same length. But my question is, do you find this a sound approach? If someone has more ideas, we could do a review of distance measures for nominal vectors!


title: Gene location clustering

text: A friend is looking at how closely related genes cluster together on a section of DNA. In particular, he is looking at genes that code for mitochondrial phenotypes (what I will call "mitochondrial genes") vs ones that don't. He has picked a distance, say D = 500,000 base pairs, and for each branch of each chromosome, computed -   1. The average density of genes within D base pairs of another gene.   2. The average density of genes on the entire branch.   3. The average density of mitochondrial genes with D base pairs of another mitochondrial gene.   4. The average density of mitochondrial genes on the entire branch. The hypothesis is that mitochondrial genes cluster together more than non- mitochondrial genes do. What statistical test is appropriate for this?


title: Hierarchical model: question on frequentist estimation

text: I am interested in understanding the differences between Bayesian and Frequentist estimation in the context of hierarchical models. Consider $n$ subjects, where for subject $i$ there are $k_i$ successes out of $m_i$ attempts, adequately modeled with a binomial distribution with parameter $\theta_i$. In a hierarchical model one could consider (if adequate for the problem at hand) $\theta_i$ coming from a Beta distribution with parameters $a$ and $b$, or a re-parametrization of them, such as mean $\mu$ and standard deviation $\sigma$. To estimate the group mean $\mu$ (random effects) in a Bayesian setting, one would consider a reasonable prior for $\mu$ and $\sigma$ and then estimate the posterior $p(\mu|k_1,\ldots,k_n,m_1,\ldots,m_n)$ by sampling or using a 2D grid for the joint posterior of $\mu$ and $\sigma$ (making use of products of beta-binomials to express the likelihood term) and then marginalizing with respect to $\sigma$. Would it be correct to state that the frequentist estimation of $\mu$ and $\sigma$ in this hierarchical model is the couple of parameters ($\mu_{ML}$,$\sigma_{ML}$) that maximizes the likelihood term (which I could obtain again on a 2D grid) or should I apply some corrections to regularize if $n$ is small? Is there an R package that I could use to get the maximum likelihood estimate for the beta-binomial model, or is it common to first apply a logit transformation and then use normal models for the population mean?


title: Performing a time series ARIMA model on natural gas power demand using the forecast package from R

text: I've been attempting to forecast natural gas power demand and how it is affected by temperature and price. I'm not sure if I have done everything correctly (relatively new to R), but I do seem to get relevant data other than I can't seem to change my forecast period, nor am I sure this is an appropriate model for this data. Hopefully someone can provide me with some guidance. Data: demand.csv               library(forecast)     data = read.csv("demand.csv")          # Create matrix of numeric predictors     xreg <- cbind(weather=data$Weather,price=data$Price,m1=data$M1,         m2=data$M2,m3=data$M3,m4=data$M4,m5=data$M5,m6=data$M6,     m7=data$M7,m8=data$M8,m9=data$M9,m10=data$M10,m11=data$M11)          # Rename columns     colnames(xreg) <- c("Weather","Price","Jan","Feb","Mar","Apr",     "May","Jun","Jul","Aug","Sep","Oct","Nov")          # Variable to be modelled     demandTS <- ts(data$Demand, frequency=12)          # Find ARIMAX model     demandArima <- auto.arima(demandTS, xreg=xreg)     demand.fcast <- forecast(demandArima, xreg=xreg)     plot(demand.fcast)      Thank you for any help. References: How to setup xreg argument in auto ARIMA in R From auto ARIMA to forecast in R


title: How do you apply fractal models to time series data?

text: Let's say I have access to both R and Python. How could I use these tools to build predictive (or at least descriptive) models of a time series? I am looking not only for the tools but helpful tutorials or lectures about how to discover the Hurst exponent, scaling exponents, and other characteristics. Thanks!


title: How to detect discernible structure in scatter plots?

text: If I have a large number of relationships I want to check that may or may not be linear (as in the bottom row of image below), is there a method that will reject the hypothesis the scatter plot is showing white noise? I want to check for the presence of any relationship without knowing what the form of that may be. What method/algorithm of measuring deviations from white noise captures the most different types of deviation? ![enter image description here](http://i.stack.imgur.com/AjMDQ.png) image from: https://en.wikipedia.org/wiki/Correlations


title: What test should I use to determine if a policy change had a statistically significant impact on website registrations?

text: A client's website was operating under a certain policy for membership sign ups for over a year. At the start of October 2012 the client implemented a new policy for sign ups that was supposed to encourage more registrations. This policy was applied site wide so there is no control group to compare against (except for the before period). I have access to the entire history of sign ups and I am wondering if there is a test that I could use to determine if there is a statistically significant difference in the number of registrations before and since the policy change? Could I do something as simple as average sign ups per day before and after the policy change and do a t-test?


title: Goodness of fit for possibly dependent data

text: My data consists of about a hundred animal visits to three different feeding zones, and I'm trying to determine if there was a preference for one zone over the others. Each animal visits between one and about a hundred zones (think hummingbirds and feeders). The problem is that 1) the sequence of zones a given animal visits are likely not independent because each individual animal may have its own preferences, and 2) the visits of different animals are likely not independent, because it is possible (likely) that, say, the first animal and the fifth are the same (they are not identifiable), so the two sequences of visits to zones by the first and fifth may be highly correlated. If everything was independent, then I could use a chi square test. Any suggestions? Edits based on the comments: The data is like:               bird1<-c(1,2,1,1,3,2,1,2)     bird2<-c(3,2,1,2,1,2,3,2,1,3,2,3,2,1)     bird3<-c(1)      ...               birdn<-c(2,3,1)      where 1,2 & 3 are labels for the zones visited before the bird left the area. Now any bird may have its own preferences. And bird1 may be bird3, because once they leave the area, if they come back we have no way of knowing it was the same bird (once it is out of eyesight). The null hypothesis would be that the three zones are equally preferred; the alternative hypothesis would be that some zone is preferred. If everything was independent, I could do a chi-square test comparing the observed frequencies to `c(1/3,1/3,1/3)`, but there may be dependence here, and simulations may prove difficulty because there is no obvious model for the dependence.


title: Uncertainty on zero counts for binned result

text: I do a counting experiment where I count observations as a function of two float parameters $x_1$ and $x_2$. This leads to a two-dimensional histogram where each bin corresponds to the number of observations with $x_1$ and $x_2$ in some range. I now see a lot of bins with zero counts in them (even though more counts are in principle possible) and wonder how to best describe the uncertainty in these bins. Knowing that I had $N$ observations which can fall into $k$ bins, I could assign an "expected rate" $n=\frac{N}{k}$ per bin and calculate the uncertainty with Poisson statistic. I however also know that my underlying process does not distribute events flat in $x_1$ and $x_2$ (but with a distribution I would like to extract from the measurement), how can I consistently assign an uncertainty in zero-count bins? I do not want to rebin my data or drop these zero bins to not loose information.


title: Basic question about gradient descent

text: I am trying to understand how gradient descent works. My understanding is that the update rule works as follows: $$ x_{n+1} = x_n - \alpha f'(x_n) $$ This seems magical to me.   1. How did this update rule come about?    2. How do I know $x_{n+1}$ lies on $f(x)$?  A simple explanation will be appreciated. Simple proof and/or pictorial description would be nice too!


title: How to determine which category is most likely given an observation with a set of characteristics

text: I have a question concerning the data analysis methods to use for a specific situation. Here is the situation: There is a dataset from an e-commerce site about its customers' purchases with 875 observations total. Each observation consists of 5 values. Scales of value measurement available for each observation (for each customer) are summarised in the following:   1. Package Type (Nominal): Type 1, Type 2, Type 3   2. Sex (Nominal): Male, Female   3. Age Group (Ordinal): NULL, Younger than 20, 20-25, Older than 25   4. Location (Nominal Scale): NULL, Region 1, Region 2   5. Order Count (Ratio Scale): Integers  NULL represents missing data. The task is to identify which Package Type is preferred by which client type, composed of Sex, Age Group, Location and Order Count. Put another way: what is the Package Type is most likely given the set of characteristics, consisting of Sex, Age Group, Location and Order Count? What am I asking for is not a ready solution for this problem - this just wouldn't be so interesting :). I want you to head me towards **the methodology** that would be used in answering this question. What **branch of Statistics** might handle this problem? Maybe you could advise me some good **classic book** covering the subject or the forum thread?


title: Mann-Kendall trend test of a huge time-series in R

text: I'm using the Mann-Kendall function of the Kendall package in R to compute the statistics of the Mann-Kendall trend test of a huge time-series (19 millions elements). It has been running for 22 hours and it hasn't yet finished. Can you suggest a faster approach?


title: Marginal effect of Probit and Logit model

text: Can anyone explain how to compute the marginal effect of Probit and Logit model in layman's terms? I am new to statistics and I am confused about these two models.


title: Cross validation with two parameters: elastic net case

text: I want to know the cross validation procedure to find the two parameters of elastic net presented by Zou and Hastie on the prostate dataset as example. I can't improve the error rate lasso with k-fold when I use elastic net.


title: How does number of observations supporting alternate hypothesis on a test of a variance have to scale so that null is rejected?

text: **Informal explanation** : In the course of my research I've run into the following problem: I am observing a machine that outputs random numbers. Most (if not all) of these random numbers come from the Gaussian distribution $\mathcal{N}(0,\sigma_0^2)$. I know the value of the variance $\sigma_0^2$. However, some proportion of the numbers I observe may come from another distribution $A$. I know that the mean of $A$ is zero, and that $A$ is symmetric (not sure if symmetry helps us). I also know that its variance is $\sigma_a^2>\sigma_0^2$, but I do not know the exact value of $\sigma_a^2$ (unlike $\sigma_0^2$.) The presence of random numbers from $A$ indicates a problem with the machine I am watching. Note that even in that case, my machine will _mostly_ output random values from $\mathcal{N}(0,\sigma_0^2)$. I would like to construct a hypothesis test on the set $\mathbf{X}$ of $n$ observations I have collected from my machine to determine whether my machine is malfunctioning. Assume that $n$ is large. Null hypothesis is that machine operates normally (i.e. $x_i\sim \mathcal{N}(0,\sigma_0^2)$ where $x_i\in\mathbf{X}$). Alternate hypothesis is that machine is malfunctioning (that is, $x_i \sim \mathcal{N}(0,\sigma_0^2)$ with $x_i\in \mathbf{X}_N$, and $y_i\sim A$ with $y_i\in \mathbf{X}_A$, where $\mathbf{X}_N\cup \mathbf{X}_A=\mathbf{X}$ and $\mathbf{X}_N\cap \mathbf{X}_A=\emptyset$.) My question is: what is the proportion $\pi$ of observations generated by $A$ in my set of observations $\mathbf{X}$ so that big given enough size of that set $n$, I achieve arbitrarily small type I and type II errors? I.e. how many out of $n$ observations in $\mathbf{X}$ have to come from $A$ for me to determine with negligible error that my machine is broken? My intuition, based on the CLT, suggests that that if more than $\sqrt{n}$ observations come from $A$, then null hypothesis will be rejected with vanishingly small sum of errors for large $n$. Is my intuition correct? If so, how can I back it up? If not, why? Does it depend on the form of $A$ (what if I told you that $A$ is the sum of $\mathcal{N}(0,\sigma_0^2)$ and another, possibly Gaussian, distribution?) **Attempt at a formal statement of the problem and conjecture** : For a zero- mean symmetric mixture distribution $(1-\pi)\mathcal{N}(0,\sigma_0^2)+\pi A(\sigma_a^2)$, I construct a test between two hypotheses: $$\begin{align} H_0:&\pi=0\\ H_1:&\pi\neq 0 \end{align}$$ with $0\leq\pi\leq1$, $\sigma_a^2>\sigma_0^2$, with $\sigma_0^2$ known. Null hypothesis states that the distribution is a known "pure" Gaussian, while the alternate hypothesis is that it's a Gaussian "contaminated" with $A$. For that test I need to characterize the value of $\pi$, in terms of $n$, that allows me to separate the two hypotheses. I conjecture that for $\pi>1/\sqrt{n}$ I can reject $H_0$ with arbitrarily small sum of probabilities of type I and type II errors given arbitrarily large (but finite) number of observations $n$. I am ok with $n$ depending on $\sigma_0^2/\sigma_a^2$ -- I assume that those are given to me and do not change as sampling progresses (I think given large enough $n$, whether I can separate the hypothesis depends only on $\pi$, but please correct me if I am wrong.) I think someone in the statistics community has to have looked at this, and the answer is out there, as this seems like a problem that would come up fairly often. Unfortunately, I've started getting into mathematical statistics very recently: my background is computer science with some graduate-level course work on probability theory, signal theory, and information theory. I would appreciate any help you have to offer as I learn this. **Attempt at solution:** I have tried the following track for solving this: suppose I collect $n$ observations in $\mathbf{X}$ and use chi-square test for testing the following hypothesis: $$\begin{align} H_0:&\sigma^2=\sigma_0^2\\ H_1:&\sigma^2>\sigma_0^2\\ \end{align}$$ If $\sqrt{n}$ observations from $A$ separate the hypotheses for this test, then I've got my proof. I can compute the statistic for the test (sample variance times $(n-1)/\sigma_0^2$). Given a very large $n$ that I have, it becomes normally- distributed by CLT. Now, let's assume that out of $n$ observations, $\sqrt{n}$ came from $A$ (rest are from $\mathcal{N}(0,\sigma_0^2)$. My test statistic then is a sum of two random variables: the sum of squares of $\sqrt{n}$ values coming from $A$ and the sum of squares of $n-\sqrt{n}$ values coming from $\mathcal{N}(0,\sigma_0^2)$. These two random variables are also (approximately) normally-distributed (very large $n$). Then I think I can look at the performance (i.e. errors) of the test, given the distribution of the test statistic -- I think I need to bound it asymptotically. Unfortunately, the Chernoff and Chebyshev bounds seem to be in the "wrong" direction. That's where I'm stuck. Am I on the right path? Did I make mistake(s) in my assumptions? Any ideas on how to proceed from here?


title: How to use stata for linear regression

text: Source |       SS       df       MS              Number of obs =    1797     -------------+------------------------------           F(  1,  1795) =  815.71            Model |  273.788275     1  273.788275           Prob > F      =  0.0000         Residual |  602.478503  1795   .33564262           R-squared     =  0.3124     -------------+------------------------------           Adj R-squared =  0.3121            Total |  876.266777  1796  .487899097           Root MSE      =  .57935          ------------------------------------------------------------------------------                y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------        no_claims |  -.1914365   .0067028   -28.56   0.000    -.2045826   -.1782904            _cons |  -1.960553   .0306965   -63.87   0.000    -2.020758   -1.900348     ------------------------------------------------------------------------------      I got the above data by typing `regress y no_claims`. Y is model by Y = $\alpha$ + $\beta x$ + $R$. But I am not sure where $\alpha$ and $\beta$ is represented in the above table.


title: How to compare Harrell C-index from different models in survival analysis?

text: In a dataset with survival event, I calculated Harrell C-index from three different models. Furthermore, I calculated the 95% C.I. for the three different models. So the next question is to compare the discrimination ability of the three models. Could anyone kindly suggest the method to compare the Harrell C-index?


title: So how would you include Bayesian estimates in a meta-analysis?

text: Inspired by this question, and particular "Problem 3": > Posterior distributions are somewhat more difficult to incorporate into a > meta-analysis, unless a frequentist, parametric description of the > distribution has been provided. I've been thinking a great deal recently about incorporating meta-analysis into a Bayesian model - primarily as a source of priors - but how to go about it the other direction? If Bayesian analysis does indeed become more popular, and becomes very easy to incorporate into existing code (the BAYES statement in SAS 9.2 and above comes to mind), we should more frequently get Bayesian estimates of effect in the literature. Let's pretend for a moment that we have an applied researcher who has decided to run a Bayesian analysis. Using the same simulation code I used for this question, if they went with a frequentist framework, they'd the following frequentist estimates: `log relative risk = 1.1009, standard error = 0.0319, log 95% CI = 1.0384, 1.1633` Using a standard, all-default and uninformative priors BAYES statement analysis, there's no reason to have a nice, symmetric confidence intervals or standard errors. In this case the posterior is pretty easily described by a normal distribution, so one could just describe it as such and be "close enough", but what happens if someone reports a Bayesian effect estimate and an asymmetrical credible interval? Is there a straightforward way to include that in a standard meta-analysis, or does the estimate need to be shoehorned back into a parametrically described distribution that's as close as possible? Or something else?


title: Caret re-sampling methods

text: I am using the library "caret" in R to test various modelling procedures. The trainControl object allows one to specify a re-sampling method. The methods are described in the documentation section 2.3 and include: boot, boot632, cv, LOOCV, LGOCV, repeatedcv and oob. Although some of these are easy to infer, not all of these methods are clearly defined. What are the procedures corresponding to these resampling methods?


title: Is it enough to check residuals versus predicted values when assessing linearity assumption in multiple regression?

text: Is it enough to check residuals versus predicted values to see if the linearity assumption in multiple linear regression is satisfied? Because for many predictors, you can't visualize the plot with the outcome variable.


title: How do I develop a metric for comparisons that involves a combination of variables at different scales?

text: I am trying to select the best answer from a set of samples, using two variables at different scales. In the data below, each row indicates the accuracy of a given measurement, the number of times it occurs in the sample set, and the answer that it indicates. Accuracy:0.00886982324861, Frequency:5.0, Answer:1.0   Accuracy:0.0104663914334, Frequency:1.0, Answer:2.0   Accuracy:0.0112727390014, Frequency:1.0, Answer:2.0   Accuracy:0.0143046058573, Frequency:3.0, Answer:1.0   Accuracy:0.0251741710747, Frequency:1.0, Answer:1.0   Accuracy:0.0322055218681, Frequency:1.0, Answer:2.0   Until now, I have been ignoring the accuracy and simply selecting the answer that occurs more frequently i.e., Answer 1 occurs 9 times, Answer 2 occurs 3 times, so select answer 1. A product of frequency and accuracy produces a metric that takes both into account, but the different scales weights frequency disproportionately more than accuracy. How can I calculate a metric that allows me to use both accuracy and frequency, but accounts for their different scales? Your help is greatly appreciated, John *EDIT Thanks for your response. In the data above, accuracy is a distance measurement from the closest category to which I am attempting to quantize this data. For example, suppose I have category 1 and category 2, the sample has a value of 1.4 so I would describe it's accuracy as .4 away from category 1 and .6 away from category 2, indicating that category 1 is the answer.  In the above data, if I multiply the frequency of sample one by it's accuracy ( 5 * .0088) the result is .0443 while the second sample gives me .0104 and the fourth sample yields .0429. From the data, the most accurate sample occurs 5 times (frequency 5) so it should clearly be better than the fourth most accurate sample occurring three times, but without some type of scaling to the results are only .0014 different. What if I had a sample that occurs once but its accuracy is .9. That would have a product of .9 which is larger than the first sample.


title: How to convert abstract value into number in Stata?

text: I have a variable called aa, it has 5 values: very satisfied, satisfied,dissatisfied,very dissatisfied,neither satisfied nor dissatisfied. I need to convert the value of "very satisfied" into "5", "satisfied" into "4", and so on.   I tried this command:               xtile bb=aa,nq(5)      but the outcome is weird, because it only those 5 value into 1, 4, 5. I cannot find the number of 2 and 3 in the result. This means "very dissastified" and "dissatisfied"are converted into the same number "5", "very satisfied" and "satisfied" are converted to the same number, "1". I don't know why, who can help me out?


title: In which setting would you expect model found by LARS to differ most from the model found by exhaustive search?

text: A bit more info; suppose that   1. you know before hand how many variables to select and that you set the complexity penalty in the LARS procedure such as to have exactly that many variables with non 0 coefficients,   2. computation costs are not an issue (the total number of variable is small, say 50),   3. that all the variables (y,x) are continuous. In what setting would the LARS model (i.e. the OLS fit of those variables having non zero coefficients in the LARS fit) be most different from a model with the same number of coefficients but found through exhaustive search (a la regsubsets())? Edit: I'm using 50 variables and 250 observations with the real coefficients drawn from a standard gaussian except for 10 of the variables having 'real' coefficients of 0 (and all the features being strongly correlated with one another). These settings are obviously not good as the differences between the two set of selected variables are minute. This is really a question about what type of data configuration should one simulate to get the most differences.


title: How to model time-varying correlation

text: Suppose I have two time-series variables, $\{x_t\}$ and $\{y_t\}$, where $t\in[1,T]$. I would like to model the correlation $\rho(x_t,y_s)$ as some function of $t$,$s$, and the difference $t-s$. In other words, $\rho_{t,s}$ may take on a different value for any valid combination of $(t,s)$, a total of $\frac {T(T+1)}2$ correlations, but I would like to economize on the number of estimated correlations (as well as possibly improving the output) by applying some sort of model. In the actual application I have in mind, these are macroeconomic and/or financial variables. The derived correlations are used to derive a full pseudo-correlation matrix, which is transformed into the nearest true P.S.D. correlation matrix using Higham's (2002) algorithm.


title: Fitting multivariate, natural cubic spline

text: _note:_ with no correct answers after a month, I have reposted to SO **Background** I have a model, $f$, where $Y=f(\textbf{X})$ $\textbf{X}$ is an $n \times m$ matrix of samples from $m$ parameters and $Y$ is the $n \times 1$ vector of model outputs. $f$ is computationally intensive, so I would like to approximate $f$ using a multivariate cubic spline through $(X,Y)$ points, so that I can evaluate $Y$ at a larger number of points. **Question** Is there an R function that will calculate an arbitrary relationship between X and Y? Specifically, I am looking for a multivariate version of the `splinefun` function, which generates a spline function for the univariate case. e.g. this is how `splinefun` works for the univariate case               x <- 1:10     y <- runif(10)     foo <- splinefun(x,y)     foo(1:10) #returns y, as example     all(y == foo(1:10))     ## TRUE      **What I have tried** I have reviewed the mda package, and it seems that the following should work:               library(mda)     x   <- data.frame(a = 1:10, b = 1:10/2, c = 1:10*2)     y   <- runif(10)     foo <- mars(x,y)     predict(foo, x) #all the same value     all(y == predict(foo,x))     ## FALSE      but I could not find any way to implement a cubic-spline in `mars` **_update_** since offering the bounty, I changed the title - If there is no R function, I would accept, in order of preference: an R function that outputs a gaussian process function, or another multivariate interpolating function that passes through the design points, preferably in R, else Matlab.


title: Fixed Effects Gravity Model for forecasting..with time constant variables Stata

text: I have a large panel data set of global trade flows and I would like to create a model for forecasting between two specific country pairs. I am just having trouble figuring out how to go about it in Stata using fixed effects and yet still accounting for the individual effects. I found this research paper and would essentially like to recreate the process explained on p.299: http://ageconsearch.umn.edu/bitstream/43996/2/martinez.pdf > A problem we faced with FEM is that we cannot directly estimate variables > that do not change over time because the inherent transformation wipes out > such variables. However, these variables can be easily estimated in a second > step, running another regression with the individual effects as the > dependent variable and distance and dummies as explanatory variables I can easily estimate `xtreg` with trade volume regressed on gdp/gdp per capita, but what are the steps I then take to estimate the individual effects (dist, island, etc) in Stata? Note: I assume they are using FE OLS, but because my data contains many 0 trade flows, is it possible to apply this method to an `xtpoisson` instead?


title: Granger causality test

text: How best can I use the Granger causality test in time series data and understand it better because I have never used it. I want to analyze long run relationship and bi-directional relationship between two variables.


title: General truncated Poisson distribution: R routines

text: I face an applicaton with (a possible) five-truncated Poisson distribution - I count periods at least five days long, thus data of the type (5, 8, 12, 5, 5, 10, ...). Much theory and software is devoted to the zero-truncated Poisson distribution. For instance, the R package pscl has the routine zeroinfl. My intention is to perform a regression analysis with the counts as dependent variable. In "Regression Analysis of Count Data" 2nd ed, (Cameron and Trivedi), a log- likelihood function is given in Eq. (4.36) for a general truncated Poisson distribution, r>0. Hence, what remains is the numerical work of optimization. My question: Are there already any packages for R for this general situation? I came across the package aster, is that a good idea? Are there others? Best regards, Carl.


title: Conditional expected value from a regression model using ordinary least squares

text: I have a query regarding part (a) of the following question. I cannot figure out how to calculate the conditional expected value of collections for the month of Easter. Is it not possible to calculate the expected value for (a) with the given information, as no dataset is given with the question, is it possible to calculate the expected value for (a) without any data / only using the info provided below? If so could somebody point me in the right direction of how this would be done? Any assistance or direction would be much appreciated. > 4. Suppose you were recently hired to be the controller for a religious > order. Since your duties include managing short-term cash flows, you need a > forecast of the timing and magnitude of offerings. You expect Easter Sunday > will be an event in which attendance and general offerings are expected to > be above normal. > > To estimate the Easter Sunday effect on order cash flows, you estimate the > following regression model using ordinary least squares:                  COLLt  = b0 + b1 TIME + b2 DEaster + et             Where:                COLLt  = Monthly collections donated for general budget;                Time = Time index;                DEaster = Dummy variable equal to 1 for the month in which Easter occurs                           and zero otherwise.      > (a) What is the conditional expected value of collections for the month in > which Easter occurs?   >  (b) Your expectations are that collections will be above normal for the > month in which Easter occurs. Describe how you would test for this effect. > Be specific in your description of testing procedures. Therefore would (a) be something like: $E(X|Y=y) Coll_{t} = B_{0} + B_{2}(1) + E_{t} X = B_{0} + B_{2}(1) + E_{t}$ or would you keep the B1 slope & time index also so answer is: $X = B_{0} + B_{1} Time + B_{2}(1) + E_{t}$ Also in regards to part (b) the only way i can think of testing this is by comparing each month against easter to determine if the collections were above normal or comparing average collections from every month but easter against the easter collections. Any assistance with this would be great, stats is so difficult.


title: Mutual information decrease with coarse-graining

text: Let $X,A,Y,B,C,D$ be random binary variables. $D$ is independent from $X,A,C$ and $C$ is independent from $Y,B,D$. Is it true that: If $I(Y:B|D=0)\leq \epsilon$ then $I(X\oplus Y:A\oplus B|C=0,D=0)\leq \epsilon \times I(X:A|C=0)$. This result seems intuitive and I managed to show it for $\epsilon=0$ and $\epsilon=1$, but maybe in the general case we should have a function $f(\epsilon)$ on the right hand side and not simply $\epsilon$. If the result is true, how would you show it using the data processing inequality, the chain rule etc.? Thanks for your help.


title: smoothing nodes values on a graph given adjacency matrix

text: I am currently looking for a method to smooth values on a graph (composed of vertices and edges). For example I have a graph with a set of nodes V and I want to be able to smooth it. I could have values on the edges, the problem would be the same using dual graph I guess. I am looking for equations I could use, given a smoothing parameter. I have tried to derived mathematically something analogous to exponential smoothing but I failed and I wonder if there is truly a close form solution.   If there is not does anyone have an idea on how to do that with random walk or diffusion processes ? I fear this would work well on a clique but not on a tree that has long and skinny leaves because it would have trouble updating the bottom of the graph.   Thanks.


title: Interpreting the lasso coefficients

text: I have used lasso logistic regression on some data and I have some non zero coefficients for some of the features. I want to know based upon the coefficient values how do I rank the features?


title: Heckman/Two part model with endogenous binary variables

text: I have a model with the continous variable(Y) that takes zero values. So, I am trying to use the limited dependent variable model (Two part/Heckman) , but I have binary endogeneous variables (two) (X1 and X2) also which I think are jointly determined. I am guessing now I have three probit (including probit part of the heckman(Two part). So, this is the trivaariate probit model. Am I in the right path? or is there any standard modelling procedure to handle this issue.


title: Interpretation of neural network outcome

text: I using R Package neuralnet on 40 input variables and output is predicting single variable Failure. Here is how neuralnet is being used               >>nn <- neuralnet(     Failure ~ MILStatus+EmissionStatus+MisfireMonitor+FuelSystemMonitor+                        CompComponentMonitor+CatalystMonitor+HeatedCatalystMonitor+EvapSystemMonitor+                        SecondAirSystemMonitor+AcRefrigerantMonitor+OxygenSensorMonitor+OxygenHeaterMonitor+                        EgrSystemMonitor+TurnSignalStatus+OilPressureLamp+BrakeIndicatorLight+                        CoolantHotLight+AbsActiveLamp+LowBattery+HighCoolant,data=m,  hidden=2, err.fct="ce",linear.output=FALSE)          >>nn     1 repetition was calculated.            Error Reached Threshold Steps     1 5984.787       0.008882457  5115          >>nn$result.matrix                                        1     error                               5.984787e+03     reached.threshold                   8.882457e-03     steps                               5.115000e+03     Intercept.to.1layhid1              -1.780980e+00     MILStatus.to.1layhid1              -5.022515e+02     EmissionStatus.to.1layhid1         -3.019886e+00     MisfireMonitor.to.1layhid1         -5.053422e+02     FuelSystemMonitor.to.1layhid1      -5.069534e+02     CompComponentMonitor.to.1layhid1   -1.782713e+02     CatalystMonitor.to.1layhid1         4.644851e-01     ....      Error is : 5984.787. How do i interpret this? Also I m not sure what is result.matrix is indicating. Just want to make sure i searched a lot on this forum and found similar questions for nnet package but not with neural net package


title: Bucketising time series

text: I have a time series where the x-axis represents hour of the day and the y-axis represents the avgerage value of some variable, say speed of all the cars in New York. I want to bucketise this data (or say cluster) on their Y-axis value. How can I decide the optimal number of buckets and what parameters should I consider to find out these optimal buckets. One simple way might be to bucket according to the speed range like 0-20 miles per hour is bucket 1, 21-30 is bucket 2, and so on. But I want to decide these ranges using the time-series.


title: Does ensembling (boosting) cause overfitting?

text: I'm using SPSS Statistics Base 20. Using Analyze $\rightarrow$ Regression $\rightarrow$ Automatic Linear Modeling I've input about 50 variables. When using no boosting, the reported accuracy of the model is 21%. The software does not tell me what is meant by accuracy, though. When I then enable boosting and set the amount of ensembles to 400, it takes about half an hour to compute, and finally I get a model where all variables are used with an accuracy of about 70%. Is it likely that this bump in accuracy is because of overfitting, or does this boosting technique really make my model more accurate?


title: How to get last available row in SPSS?

text: I have the following column: ## Growth               0.1     0.5     0.23     .     0.5     0.67     0.23     0.42     .     .     .     0.12     0.20      Where dots indicate that there was no data available at that time. Now I want to compute the difference in growth between the current case and the **last available case**. Thus, in the second last case this should become: `0.12 - 0.42 = -0.30`. `Lag(growth)` takes the last case, and returns nothing if it is not defined. So how do I get the last available case? **Response** : This seems to work:               define !lastAvailable(target !Tokens(1) /level !Tokens(1) /Varys !CMDEND)        compute id= $casenum.        dataset copy tmpset2.        select if not miss(!target).        compute diff= !target - lag(!target, !level).        match files /file= * /file= tmpset2 /by id.        exec.        dataset close tmpset2.     !enddefine.          !lastAvailable target=growth level=3


title: Pattern finding approaches

text: I'm looking for patterns in a signal. I had used successfully a cross- correlation approach. However, I found some references to Gibbs sampling approaches for finding patterns. My question is: Are both approaches complementary one to the other, or are they redundant in some way?


title: Estimation of AR(1) process

text: Suppose the stochastic process ${X_t}$ satisfies the equation $$X_t=\phi X_{t-1} + Z_t \tag{A}$$ where $\phi>1$ and $Z_t$ is a white noise. Then iterating forward we get that the only stationary solution of the equation above can be written as: $$X_t=-\sum_{j=1}^{\infty}Z_{t+j} {\phi}^{-j}$$ Therefore, we can write $$X_t-\phi^{-1} X_{t-1} = -\sum_{j=1}^{\infty}Z_{t+j} {\phi}^{-j} +\sum_{j=1}^{\infty}Z_{t-1+j} {\phi}^{-j-1}=\phi^{-2}Z_t$$ And thus, defining $\tilde{Z}_t=\phi^{-2}Z_t$, we have that $X_t$ also satisfies the following AR equation: $$X_t=\phi^{-1}X_{t-1}+\tilde{Z}_t\tag{B}$$ My question is simple. Suppose I have some data on $X_t$ and I run an OLS regression with one lag and no intercept on it. What solution will I estimate, (A) or (B)? Why?


title: Standard error of proportion that takes into account population size

text: Is there another method to estimate standard errors of proportions other than $\sqrt{p(1-p)/n}$ where an estimated total population size is taken into account? For example if I sample 100 out of population of 100, I should have a standard error of 0. If I sample 100 out of a billion, it should have a higher error. I read in a blog that the formula above should be used when the population is at least 10 times bigger than the sample, but I don't know where this came from or if it is true or what happens if a survey is close to that 1/10 sample.


title: What's MGD sequence in time series?

text: What's MGD sequence in time series theory? what's the full name of MGD, I google it and can't find its definition.


title: How do you specify a binary time-series-cross-section model in lme4?

text: I have a standard binary time-series–cross-section (BTSCS) model that I would like to specify as a mixed effects model using the lme4 package. I've read elsewhere that time-series–cross-section (TSCS) models can be conceptually understood as multilevel models (i.e., years [time] are nested within countries [cross-section]). How would you specify the model if the DV is binary, which coincides with a broader issue of temporal dependence? Assume the standard BTSCS model as it's routinely applied in Stata. For example:               logit y x1 x2 x3 spline1 spline2 spline3, cluster(country)      How would this be specified in `lme4`? Are `year` and `country` separate random effects? Are they nested? What about the cubic splines and the issue of temporal dependence (see: Beck, Katz, and Turner (1998), PDF)? I'm proficient in using `lme4`, but, to date, the random effect had been simple and easy to identify. Any help in tackling a more complicated puzzle would be appreciated.


title: How to find weber fractions from a confusion matrix

text: I have a confusion matrix for a 4 class classification problem which looks this: ![enter image description here](http://i.stack.imgur.com/2jRMo.png) My hypothesis is that group 3 is easiest to classify. Somebody suggested me that I can find weber fractions from this confusion matrix and perform further statistical analysis. I am unable to understand how I can do the same. I am wondering if anybody can suggest me how weber fraction could be estimated or any other statistical test can be applied to this confusion matrix. Thanks in advance.


title: Filtering a dataframe

text: Still learning basic functions in R, The subset function seems to only filter based a condition based on single column with or without multiple conditions? How can I easily filter data from a dataframe?   1. when you are provided with multiple conditions   2. When the condition needs to be applied across the available columns. Example: Given a data frame containing               name    D1      D2     D3      D4     julius  "A"     "A"    "B"     "B"     cate    "D"     "E"     "A"     "C"     karo    "A"     "D"     "C"     "E"      say I want to filter this dataframe so that only names where any of the D1 to D4 is an 'E' then I should have,               name    D1      D2     D3      D4     cate    "D"     "E"     "A"     "C"     karo    "A"     "D"     "C"     "E"      Say that the D1 can be a big list of columns, how or what is the recommended approach to perform this filter? Thank you


title: Questions regarding correlation

text: I'm writing a thesis about cooperating competition between companies. I have a few problems with the last two research questions: RQ7: Are there any differences between the perception of costs and benefits of the cooperating competition between the companies that have previously experienced the concept of cooperating competition and those who had not. RQ8: Are there any differences between the perception of costs and benefits of the cooperating competition between different kinds of industries in which the companies conduct their buisiness. So, for RQ7 I have 3 major perceptions of costs which are scales ranging from 1(I highly don't agree) to 7(I highly agree), the same goes for perceptions of benefits and i have a nominal variable for previously experienced cooperating competition, which is just yes or no. RQ8 is similar to RQ7 except it has a nominal variable with eleven possibilities for different kind of industries. I used custom table, dragged perceptions in columns and previous competition in row and I got the means for every perception for yes and for no (previous competition) and for different kinds of industries. Looking at the results I don't know if I answered my questions, is there maybe a better tool to for answering these questions, I'm thinking about correlation of some sorts, but don't know how to put it all together. please take a look and tell me what you think. Thank you!


title: Predicting number of events with 99.9% probability based on tests of four devices

text: I've never taken a stats course so I really don't know where to begin on this. I am using microcontroller that according to the datasheet, has a Flash memory that can be reprogrammed a minimum of 10,000 times, with 100,000 times being a more typical number. These numbers seem low to me (based on other devices using a similar technology which claim lifetimes of an order of magnitude higher), so I decided to take one of our products and do destructive testing on it. I started testing last Friday night and it is already up to 600,000 events without an error (as determined by calculating a 16-bit CRC over the 4K bytes of data test being Flashed). Based on this device far surpassing the "typical" value in the datasheet, I have started tests on three more devices. So at some point I will have a failure value for all four devices. How can I use that data to predict an expected lifetime x for the device with a 99.9% certainty? (Not sure if I am saying that right, what I mean is that only 0.1% of devices will fail before x events).


title: Can the multiple linear correlation coefficient be negative?

text: I am using IDL regression function to compute the multiple linear correlation coefficient...               x = [transpose(anom1),transpose(anom2),transpose(anom3)]     coef = regress(x,y, const=a0, correlation=corr, mcorrelation=mcorr, sigma=stderr)      `mcorr` is returned with values between 0.0 and 1.0. Clearly, the result of the following IDL code I dug out of the can can't be negative...               mcorrelation = SQRT(mult_covariance)      But, could you help me verify this? I would like a citation or to see 'the math' behind calculating `mcorrelation`.


title: Statement of result for binomial confidence intervals

text: Confidence intervals for binomial proportions have irregular coverage over the range of possible population parameters (e.g. see Brown et al. 2001 http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.ss/1009213286). How can I formally and usefully describe the properties of the confidence intervals? Say I toss a coin ten times and obtain seven heads. Are the following statements accurate? For the Clopper-Pearson method:   The interval 0.3475–0.9333 has been generated by a method that will on at least 95% of occasions, for any true population proportion, contain the true population proportion. The long-run frequency with which this method would yield confidence intervals containing the true population proportion pertaining to this particular experiment is at least 95%. For the Wilson's scores method:   The interval 0.3968–0.8922 has been generated by a method that will on 95% of occasions, averaged over all population proportions, contain the true population proportion. The long-run frequency with which this method would yield confidence intervals containing the true population proportion pertaining to this particular experiment may be more or less than 95%.


title: Number of parameters in a Markov chain

text: I was reading the book Pattern Recognition and Machine Learning by Bishop, which stated that for the first order Markov chain with K states, the number of parameters is K(K-1). Why is that? I think it should have been just K-1, because if we have $p(x_n|x_{n-1})$ for n=2 to K that will only be K-1, correct?


title: How detailed should a data-driven Bayesian prior be?

text: My exact problem is this. I have a number of sources of traffic with different conversion rates. I have good evidence that conversion rates vary based on the source. For each traffic source I have how many leads we got, and how many converted. I want an estimate of each source's conversion rate. I do not need individual sources to be accurately estimated (this would be impossible given that many traffic sources only have 0-5 measured conversions). Instead I would like to be able to say that, across all traffic sources estimated at 1%-1.5% average conversion rates, the average future conversion rate is likely to be somewhere in that general ballpark. The obvious approach of using past performance to predict future performance for each individual source fails badly because of the well-known phenomena of regression to the mean. For instance, among small traffic sources with no conversions there is a reasonable expectation of future performance being better than that. And likewise small traffic sources with lots of conversions are unlikely to maintain their record. My naive idea is to take my data, and use it to produce a reasonable Bayesian prior for the true conversion rate of a random traffic source. Then for each source I can start with that prior and produce a posterior distribution for the true conversion rate of that source. And then my estimate of the average conversion rate for that source will be the average of the posterior. My initial idea for how to fit data to a reasonable prior works like this. My prior will be the sum of piecewise linear functions, which will individually be 0 to a starting point, rise to a midpoint, then fall to the ending point after which it is 0 again. For a given division of the likely range of conversion rates into these intervals and midpoints, the prior I'd produce would be the one to maximize the sum of the logs of the likelihoods that each traffic source would have the conversions that it did. My problem is that the more pieces I divide my interval into, the more closely I can fit my existing data set. But at some point I'm clearly over-fitting. Are there any guidelines that I can use to get a sense when I'm over-fitting my prior to my data? I'm thinking there might be something like a statistic that I can compute to test how closely my measured data fits the prior - if the measured data has a better fit to the prior than a "random" data set should, then I've probably gone too far. I would be grateful if anyone can suggest a statistic, an alternate rule of thumb to avoid over fitting, or a different approach to the original problem. Since I do not have access to a university library, please only suggest specific books or papers behind paywalls if they are guaranteed to be very relevant to my problem.


title: Very high correlation between two variables in SEM as compared to SPSS

text: I am testing a model which involves three variables. I am testing a partial mediation model. In SPSS, when I calculate bivariate correlation among these three variables, the maximum correlation (say between B&C) I get is -.530 . However, using AMOS and testing for discriminant validity, I get the correlation as -.814. This is very high and is creating problems for discriminant validity. Other two correlations between A & B and between C & A remains more or less the same (-.314 and .359 vs -.369 and .429) I am not able to understand why it is happening and what I should be doing to control it.


title: What is the interpretation of interquartile range?

text: I have daily measurements of nitrogen dioxide for one year (365 days) and the **interquartile** (IQR) is 24 microgram per cubic meter. What does "24" mean in this context, apart from the definition of IQR which is the difference between **the 25th and 75th percentile**? How would you explain this figure to a journalist, for example? Thanks


title: Pseudo Inverse Instead of Inverse with LDA?

text: I have implemented the LDA algorithm. However when I had to get the inverse of a matrix Matlab threw an error and I had to use pinv (pseudo inverse) instead of inverse. Did I do something wrong or is it usual sometimes to have to use pinv instead of inv? **PS:** My matrix is a square matrix as usual.


title: Two questions about MLE

text: $$f(x)=\alpha\beta x^{\beta-1}e^{-\alpha x^{\beta}}$$ ![enter image description here](http://i.stack.imgur.com/6UECS.jpg) 1) How do we know that n must be multiplied by both alpha and beta only after we take the log, why don't we simply multiply by n in the original function and then take the log? 2) In the derivative with respect to $\beta$, why do we have a log in the function (after the minus sign)? Where does this come from?


title: Logistic regression without negative samples

text: I have a data set of RNA reaction values of breast cancer. I want to figure out which RNAs are essential genes by Logistic Regression & LASSO. The data set has no negative samples. What should I do to ensure the method will give a reliable result?


title: What type of Design is this?

text: Participants are asked to recall either a happy or sad event in their life. Once a mood of happiness of sadness has been elicited all participants are then presented with 3 stories where the central character was either involved in a happy, sad or neutral event. The results are as follows: Happy Story:   Happy- 87, 75, 83   Sad - 55 , 65, 63 Sad Story:   Happy - 85, 78, 68   Sad - 62, 58, 48 Neutral Story:   Happy - 87, 82, 47   Sad - 48, 63, 47 Could someone please tell me what kind of design this is as I have to put it into SPSS and do appropriate tests on it.


title: What am I supposed to do if Cronbach's alpha is negative?

text: I have gathered data from 130 questionnaires and I have 30 variables. Cronbach's alpha is badly negative. I have checked everything I could think of, but the result did not change. SPSS file is attached http://www.wikiupload.com/9M9YHEV6I4Q1HPR


title: Multistage sampling in R

text: I've got a dataset similar to this:               pat_id   epis Care Type     1       1722650 Acute Care     1       1723120 Rehabilitation care     2       1584309 Acute Care     2       1585705 Rehabilitation care     3       1726487 GEM     3       1664031 Acute Care     3       1726488 Acute Care     3       1726489 Rehabilitation care      Each patient has multiple "episodes/care types". I want to sample 50 patients but also approximately proportional to the "Care Type" from the population (say, 50% acute care, 30% rehabilitation care, 20% GEM). One way I thought of doing is to split the data set into say "acute", "rehab" and "gem" then sample 25 from "acute", 15 from "rehab" and 10 from "gem". But there would be an overlap using this approach "e.g. Patient 1375 would be in gem, acute and rehab. Is there a R package that would handle this sort of sampling easily?


title: Is the distribution of the minimum of two other distributions a mixture distribution? Or is there a better term?

text: This is a terminology question motivated by a review that I got on a paper. In the following I believe that $y$ would be considered to be distributed according to a _mixture distribution_ : $$y \sim \left\{ \begin{array}{ll} f_1(\theta_1) & \mbox{with } p ~ \mbox{probability} \\ f_2(\theta_2) & \mbox{with } (1 - p) ~ \mbox{probability} \end{array} \right. $$ where $f_1$ and $f_2$ are probability distributions. **But would it be correct to say that $x$ is also distributed according to a mixture distribution?** $$ x \sim \text{min}(z_1,z_2) \\ z_1 \sim f_1(\theta_1) \\ z_2 \sim f_2(\theta_2)$$ If not, is there a better term? Or, at least, a general term for a distribution that is _some_ combination of two other distributions? (I mean, $x$ is certainly the result of "mixing together" $f_1$ and $f_2$ in some way, but I don't know if this is enough to call it a mixture... :)


title: Interpret residuals in random and fixed effects models

text: How to interpret residuals in random and fixed effects models? I am a medical doctor. Kindly explain in simple words. **Updated** Here is an example:               Weight (Fixed)      Weight (Random)    Residual (Fixed)      Residual (Random)     Relative weight    Relative weight     Std Residual         Std Residual          32.27                21.56               -4.61                -1.62     25.14                21.02               -0.41                -0.46     18.52                20.21               2.16                  0.52     14.39                19.39               2.10                  0.62      9.68                17.81               2.56                  1.04


title: Using Holt-Winters for forecasting in Python

text: [I first posted this question to Stack Overflow here but didn't get any replies, so I thought I'd try over here. Apologies if reposting isn't allowed.] I've been trying to use this implementation of the Holt-Winters algorithm for time series forecasting in Python but have run into a roadblock... basically, for some series of (positive) inputs, it sometimes forecasts negative numbers, which should clearly not be the case. Even if the forecasts are not negative, they are sometimes wildly inaccurate - orders of magnitude higher/lower than they should be. Giving the algorithm more periods of data to work with does not appear to help, and in fact often makes the forecast worse. The data I'm using has the following characteristics, which might be problems:   * Very frequently sampled (one data point every 15 minutes, as opposed to monthly data as the example uses) - but from what I've read, the Holt-Winters algorithm shouldn't have a problem with that. Perhaps that indicates a problem with the implementation?   * Has multiple periodicities - there are daily peaks (i.e. every 96 data points) as well as a weekly cycle of weekend data being significantly lower than weekday data - for example weekdays can peak around 4000 but weekends peak at 1000 - but even when I only give it weekday data, I run into the negative-number problem.  Is there something I'm missing with either the implementation or my usage of the Holt-Winters algorithm in general? I'm not a statistician so I'm using the 'default' values of alpha, beta, and gamma indicated in the link above - is that likely to be the problem? What is a better way to calculate these values? Or ... is there a better algorithm to use here than Holt-Winters? Ultimately I just want to create sensible forecasts from historical data here. I've tried single- and double-exponential smoothing but (as far as I understand) neither support periodicity in data. I have also looked into using the R forecast package instead through rpy2 - would that give me better results? I imagine I would still have to calculate the parameters and so on, so it would only be a good idea if my current problem lies in the implementation of the algorithm...? Any help/input would be greatly appreciated!


title: What research tool to use when researching a project?

text: I am currently working towards a Statistics BSc and am frequently having to organise my research. I am currently using a mixture of a mediawiki, and a bibliography browser plugin, to make an archive of my online data sources, and references to pdfs with useful information. but its becoming a bit of a mess. Is there any research tool that is orientated towards statistics? would-likes include   1. //TODO, task tracking so I can have a central list of outstanding tasks   2. some cloud storage, so I can access my archive from any machine some suggestions have been   * org-mode   * zotero   * mediawiki (I guess that anything I could store the files as text, I can distribute my work using SVN/git etc)